Building archean cratonic roots

Jain, Charitra; Rozel, Antoine B.; van Hunen, Jeroen; Chin, Emily J.; Manjón-Cabeza Córdoba, Antonio

doi:10.3389/feart.2022.966397

ORIGINAL RESEARCH article

Front. Earth Sci., 22 November 2022

Sec. Solid Earth Geophysics

Volume 10 - 2022 | https://doi.org/10.3389/feart.2022.966397

Building archean cratonic roots

Charitra Jain^1,2*

Antoine B. Rozel³

Jeroen van Hunen¹

Emily J. Chin⁴

Antonio Manjón-Cabeza Córdoba^3,5

¹Department of Earth Sciences, Durham University, Durham, United Kingdom
²GFZ German Research Centre for Geosciences, Section of Geodynamic Modelling, Potsdam, Germany
³Department of Earth Sciences, Institute of Geophysics, ETH Zurich, Zurich, Switzerland
⁴Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA, United States
⁵Instituto Andaluz de Ciencias de la Tierra (IACT), UGR-CSIC, Armilla, Spain

Geophysical, geochemical, and geological investigations have attributed the stable behaviour of Earth’s continents to the presence of their Archean cratonic roots. These roots are likely composed of melt-depleted, low density residual peridotite with high magnesium number (Mg#), while devolatilisation from the upper mantle during magmatic events might have made these roots more viscous and intrinsically stronger than the convecting mantle. Several conceptual dynamic and petrological models of craton formation have been proposed. Dynamic models invoke far-field shortening or mantle melting events, e.g., by mantle plumes, to create melt-depleted and thick cratons. Compositional buoyancy and rheological modifications have also been invoked to create long-lived stable cratonic lithosphere. However, these conceptual models have not been tested in a dynamically self-consistent model. In this study, we present global thermochemical models of craton formation with coupled core-mantle-crust evolution driven entirely by gravitational forces. Our results with melting and crustal production (both oceanic and continental) show that formation of cratonic roots can occur through naturally occurring lateral compression and thickening of the lithosphere in a self-consistent manner, without the need to invoke far-field tectonic forces. Plume impingements, and gravitational sliding creates thrusting of lithosphere to form thick, stable, and strong lithosphere that has a strong resemblance to the Archean cratons that we can still observe today at the Earth’s surface. These models also suggest the recycling of denser eclogitic crust by delamination and dripping processes. Within our computed parameter space, a variety of tectonic regimes are observed which also transition with time. Based on these results, we propose that a ridge-only regime or a sluggish-lid regime might have been active on Earth during the Archean Eon as they offer favourable dynamics and conditions for craton formation.

1 Introduction

The relative stability and longevity of Earth’s continents compared to its ephemeral oceanic crust are attributed to the presence of their buoyant continental lithosphere, including the cratonic regions that have resisted deformation for several billion years. Cratons are underlain by cold and thick mantle roots (Pearson et al., 2021), which can extend up to several hundreds of kms deep as evidenced by xenolith P-T data (Lee et al., 2011) and observations of low surface heat flux (Pollack et al., 1993; Nyblade, 1999) and high seismic velocities (Jordan, 1979; Grand and Helmberger, 1984). These roots are generally considered to be composed of residual peridotites, which can be classified into fertile lherzolitic (primitive mantle or pyrolytic) and highly melt-depleted harzburgitic lithologies (Boyd, 1987, 1989). It has been proposed that progressive melt extraction from the igneous protolith increases the magnesium number Mg# in these residues, thereby decreasing their densities (Jordan, 1979; Lee, 2003; Schutt and Lesher, 2006). On average, cratonic mantle is cooler than the ambient mantle, which would result in a density increase owing to thermal contraction. Jordan (1978, 1988) proposed that under perfectly isopycnic conditions, the positive compositional buoyancy emanating from these residual peridotites would exactly offset the negative thermal buoyancy of the colder cratonic mantle. However, having this net neutral buoyancy may not be enough to ensure cratonic stability. Pollack (1986) argued that devolatilisation results in a more refractory and stiff residual mantle, thereby making it less susceptible to subsequent melting and entrainment in the convecting mantle. Geological record suggests that the efficiency of devolatilisation and partial melting have decreased over time in response to the secular cooling of the mantle (Herzberg et al., 2010). Experimental results have also shown a substantial reduction in the viscosity of olivine in the presence of water (Hirth and Kohlstedt, 1996). Volatiles are highly incompatible in solids during melting, and therefore it has been hypothesised that cratonic roots are made of melt-depleted and dehydrated peridotites, which are considerably more viscous and intrinsically stronger than the convecting mantle (Hirth et al., 2000; Lee et al., 2011).

The mechanisms behind the formation and stabilisation of these cratonic roots, however, remain debated (Pearson et al., 2021) and three separate hypotheses have been proposed with each having their strengths and caveats (Pearson and Wittig, 2008; Lee et al., 2011). One model involves high degrees of melting within a large thermal plume (Herzberg, 1993; Griffin et al., 2003; Griffin and O’Reilly, 2007; Arndt et al., 2009) leaving behind a melt-depleted, dehydrated, and low-density residue, which would explain the higher viscosity and the buoyancy observed in the roots. This model predicts an increase in peridotite fertility (low melt-depletion and Mg#) with depth as higher amounts of melt would be extracted at shallower depths during decompression melting. However, such gradual chemical stratification has not been widely observed in cratons (Lee et al., 2011). The plume model also fails to explain the subsolidus history of pressure increase experienced by the cratonic peridotites (Lee and Chin, 2014). Another model involves underthrusting or stacking of subducted oceanic lithosphere (Helmstaedt and Schulze, 1989; Canil, 2004, 2008; Simon et al., 2007; Pearson and Wittig, 2008), which could explain the low-pressure protoliths of cratonic peridotites and the lack of compositional stratification with depth, but it also has several drawbacks. The negatively buoyant oceanic lithosphere is more likely to subduct than stack (Arndt et al., 2009), unless fragments of the downgoing plate are captured as it has been suggested for the Farallon oceanic plate (Saleeby, 2003; Luffi et al., 2009; Lee et al., 2011). Large-scale weak zones or fault planes in the supposedly stronger continental mantle might have been necessary for imbrication of oceanic lithosphere to work. Lee et al. (2008) proposed that serpentinisation of the oceanic mantle during the mid-Archean to the early Proterozoic might have provided the weak zones. Furthermore, a mechanism is needed to remove the eclogite from the oceanic lithosphere as the predicted amount is inconsistent with its rare presence among the mantle xenoliths in kimberlites (Schulze, 1989). The third model involves making cratonic roots by accretion and thickening of already buoyant arc lithosphere (Şengör et al., 1993; Ducea and Saleeby, 1998; Kelemen et al., 1998; Parman, 2004), which might have been active during the Phanerozoic. For this mechanism to work, the arc has to be mature enough to be in a compressional state and arc magmatism has to coincide with lithospheric thickening (Kay et al., 2005; DeCelles et al., 2009). Despite the presence of some similarities between rocks from modern continental arcs and Precambrian xenoliths, this mechanism needs further validation by integrating field data with geochemical and petrological studies (Lee et al., 2011).

Wang et al. (2018) explored the viability of compressive thickening models to make cratonic roots in their numerical experiments. They presented a two-stage growth model, in which a tectonic shortening stage of the residual mantle for 10s of Myr is followed by gravitational self-thickening stage that lasts for 100s of Myr. The authors further showed that the long-term secular cooling of the mantle prevents a Rayleigh-Taylor type collapse and helps with cratonic root preservation. Beall et al. (2018) proposed that a rapid transition from heat-pipe dominated vertical tectonics to plate tectonics during the Archean might have provided stresses high enough to thicken the cratonic lithosphere. The authors initialised their lid-breaking model with a buoyant and strong layer that undergoes compression and stabilises. Recently, models by Perchuk et al. (2020) showed oceanic subduction followed by arc-continent collision under Archean conditions. Using a prescribed convergence rate, they demonstrated protokeel formation by viscous underplating of the low-density, melt-depleted sublithospheric oceanic mantle during subduction, which subsequently matures into keels with conductive cooling. Capitanio et al. (2020) showed the formation of cratonic lithospheric mantle by melt-depletion induced dehydration stiffening beneath a long-lived stretching lid. Over time, these regions stiffen, resist deformation and force strain migration and cooling.

Besides formation, the preservation and stability of cratonic roots until present-day has also been studied. Previous numerical modelling work has shown that the positive chemical buoyancy of the cratonic lithosphere is not enough to keep it from entraining into the underlying convecting mantle (Lenardic and Moresi, 1999; Lenardic et al., 2003; Sleep, 2003). Using scaling analysis and considering a weakly temperature-dependent viscosity, Sleep (2003) showed that a viscosity increase of a factor 20 between the cratonic lithosphere and the ambient mantle would help with the craton survival. Lenardic et al. (2003) proposed that high cratonic root viscosity only helps stabilising the cratonic lithosphere when its thickness is twice larger than that of the oceanic lithosphere. Additionally, having higher brittle yield stress in the cratonic lithosphere or the presence of mobile belts (weak zones) around cratons also helps with cratonic root stability. O’Neill et al. (2008) further showed cratonic stability in their models having a compositional viscosity ratio of 50–150 between the root and the asthenosphere along with a root/oceanic lithospheric yield stress ratio of 5–30. The authors also argued that hotter Archean mantle temperatures would result in lower viscosities and reduced stresses around cratons, thereby promoting cratonic survival.

Using simulations with non-Newtonian rheology, Wang et al. (2014) showed resistance of compositionally buoyant cratonic roots against erosion by small scale convection for billion year timescales when initialised with a strengthening viscosity factor of 3 in the rheological prefactor, leading to an effective viscosity increase of about 50 (3ⁿ, with n = 3.5) when accounting for the non-Newtonian contribution. They also reported that for cratonic roots with no intrinsic buoyancy, a higher strengthening factor on the order of 10 leading to an effective viscosity contrast of over a thousand is required to ensure their stability.

Even though models have been instrumental in testing the validity of different geological and geophysical hypotheses, there are either regional scale or upper-mantle models and thus by design, lack mantle plumes and convective activity offered by global models, as well as a self-regulating thermal evolution. Extending the previous work by Jain et al. (2019), in this study, we present whole mantle convection models that are capable of generating cratons in a self-consistent manner, i.e., without any artificial forcing of surface velocities. Furthermore, we test whether lithospheric compression and tectonic thickening can occur naturally in our models and aid with the formation of cratonic roots.

2 Physical model and numerical model

To investigate how Earth built cratonic roots in a self-consistent way, we model its thermochemical evolution during the Archean Eon (4.0–2.5 Ga) with realistic parameter values (Table 1) suited for early Earth conditions. The model includes pressure-, temperature-, and water-dependent viscosity, plasticity, time-dependent radiogenic and basal heating, phase transitions, and melting leading to basaltic and felsic crust production. We introduce a temperature-, pressure-, and composition-dependent solubility function, which controls the water ingassing and outgassing. Moreover, we improve upon the existing melting parameterisation of Jain et al. (2019) by making the solidi and liquidi depend on water saturation and introduce magnesium number evolution.

TABLE 1

TABLE 1. Non-dimensional and dimensional parameters along with the rheological properties for 3 different layers i used in this study (UM = Upper Mantle (dry olivine); PV = Perovskite; PPV = Post-Perovskite). TTG = tonalite-trondhjemite-granodiorite.

2.1 Conservation equations, boundary conditions, and solution method

We use the code StagYY (Tackley, 2008) to solve the compressible anelastic convection equations in a two-dimensional spherical annulus geometry (Hernlund and Tackley, 2008). StagYY employs a Eulerian mesh with radial refinement on which equations for momentum conservation, continuity and heat diffusion/production are solved using the finite volume approach. The computational domain considered here is quadrant annulus, where the grid is divided into 512 (laterally) times 128 (radially) cells. Additionally, ∼ 1.3 million Lagrangian tracers, carrying various quantities (see below), are advected in the computational domain. Free slip boundary conditions are used at the core-mantle boundary and the surface. Side boundary conditions are periodic.

The Stokes equation is given by

0 = - \nabla P + \nabla \cdot τ + ρ g, (1)

where P is pressure, τ is the stress tensor, ρ is the density, and g is the gravitational acceleration, which has a homogeneous magnitude. Density is computed as a function of pressure, temperature, and composition including solid-solid phase transitions, following a 3rd order Birch-Murnaghan equation of state (Tackley et al., 2013).

The continuity equation is given by:

\nabla \cdot (ρ u) = 0, (2)

where u is the velocity vector.

The advection-diffusion equation is given by:

ρ C_{P} (\frac{\partial T}{\partial t} + u \cdot \nabla T) - α T (u \cdot \nabla P) = \nabla \cdot (k \nabla T) + τ : \nabla u + ρ H, (3)

where C_p is the heat capacity, T is temperature, α is the thermal expansivity, k is the conductivity and H is the radiogenic heating. τ : ∇u denotes tensor contraction, such that: τ : ∇u = ∑_ijτ_ij∂u_i/∂x_j, where x_j is the position. The second term on the left-hand side is the heat production/consumption due to adiabatic (de)compression. On the right-hand side, the first term describes heat diffusion, the second term adds viscous dissipation during non-elastic deformation processes, and the third term contributes radiogenic heating (Ismail-Zadeh and Tackley, 2009). The terms on the right hand side are calculated on the mesh, whereas advection and adiabatic temperature changes are performed on the tracer level. Only the changes in the temperature field (as opposed to the temperature itself) due to the right hand side terms are interpolated to the tracers at each time step to limit numerical diffusion.

Tracers are advected in the flow using velocities interpolated from the cell edges, where the Stokes and continuity equations (Eqs. 1, 2) are solved together with a direct PETSc solver (Balay et al., 2012) and Picard iterations of the viscosity. We found that a second order interpolation and a fourth order Runge-Kutta scheme is necessary to prevent gaps from opening in the tracer field. Tracers carry various quantities such as temperature, composition, water, radiogenic heating, time of melting, etc. Tracers are also transported instantly, when magmatism occurs, through the procedure detailed in Section 2.4.

The surface temperature (300 K) is kept constant whereas the core-mantle boundary (CMB) temperature (4500 K) evolves using the model proposed by Buffett et al. (1992, 1996), allowing for the cooling due to the CMB heat flux and taking into account the energy necessary to grow the inner core (see Nakagawa and Tackley (2004) for parameterisation details).

2.2 Rheology

A visco-plastic rheology is considered where the viscous deformation is accommodated by grain-size independent diffusion creep. The mantle is divided into 3 different layers i: upper mantle (1), lower mantle (2) and post-perovskite layer (3), with each layer having different values for activation energy E_i and activation volume V_i (Karato and Wu, 1993; Yamazaki and Karato, 2001). Following the Arrhenius law, the temperature-, pressure-, and water- dependent viscosity is given as (see Table 1 for constants):

η (T, P, w) = η_{0} Δ η_{i} A_{w} \exp (\frac{E_{i} + P V (P)}{R T} - \frac{E_{i}}{R T_{0}}), (4)

where η₀ is the reference viscosity at zero pressure and reference temperature T₀ (1600 K), Δη_i is the viscosity offset between layer i and the reference viscosity, A_w is the water-dependent viscosity multiplicative factor, P is the pressure, R is the gas constant and T is the absolute temperature. η₀ is valid for the phase system olivine and the reference composition (“pyrolite” as explained in Section 2.3), both of which have viscosity multipliers of 1 (see Table 1). To obtain profiles consistent with Karato and Wu (1993); Yamazaki and Karato (2001), the activation volume decreases exponentially with increasing pressure in each layer i according to the relation:

V (P) = V_{i} \exp (- \frac{P}{P_{i}}) . (5)

where P_i is the pressure scale which is different for each layer i as given in Table 1. Based on the experimental data and assuming that the activation parameters are independent of the water content (Karato, 2015), the presence of water modulates the viscosity by factor A_w as:

A_{w} = {[Δ η_{w} + (1 - Δ η_{w}) {(\frac{C_{w}}{C_{w 0}})}^{r_{w}}]}^{- 1}, (6)

where Δη_w is the viscosity offset due to water (or dehydration strengthening), C_w is the cell water concentration, C_w0 is the reference water concentration (100 ppm), and r_w = 1. In this work, the viscosity of the material drier and wetter than 100 ppm are increased and decreased respectively. However, for water concentrations higher than 100 ppm, the weakening would be the same irrespective of Δη_w. Further viscosity jumps of 10 and 0.1 are applied at the upper-lower mantle transition (660 km) (Čížková et al., 2012) and perovskite-post-perovskite transition (2740 km) (Hunt et al., 2009; Ammann et al., 2010) respectively.

To allow for plastic deformation in the lithosphere, plastic yielding is assumed (Moresi and Solomatov, 1998; Tackley, 2000). The maximum stress that a material can sustain before deforming plastically is given by the yield stress σ_Y, which has both brittle and ductile components:

σ_{Y} = \min (σ_{Y,ductile}, σ_{Y,brittle}) . (7)

The ductile yield stress σ_Y,ductile increases linearly with pressure as:

σ_{Y,ductile} = σ_{Y}^{0} + \overset{́}{σ_{Y}} P, (8)

where $σ_{Y}^{0}$ is the surface ductile yield stress and $σ_{Y}^{'}$ is the pressure gradient of the ductile yield stress. Following Byerlee (1978), the brittle yield stress σ_Y,brittle is calculated as:

σ_{Y,brittle} = μ P, (9)

where μ is the friction coefficient. If the convective stresses exceed the yield stress, the viscosity is reduced to the yielding viscosity $η_{Y} = σ_{Y} / 2 \dot{ϵ}$ , where $\dot{ϵ}$ is the 2nd invariant of the strain-rate tensor. The effective viscosity is then given by:

η_{eff} = {(\frac{1}{η} + \frac{2 \dot{ε}}{σ_{Y}})}^{- 1} . (10)

To mitigate large viscosity variations and ensure code stability, viscosity limiters are applied. On the lower end, viscosity is limited to 10¹⁸ Pa⋅s. For the upper end, different values are used in this study. All simulations considered a reference viscosity η₀ = 10²⁰ Pa⋅s and initial mantle potential temperature T_P0 = 1900 K. Additionally, melt-depletion induced density reduction of the residual harzburgite (peridotite) is considered implicitly by including phase changes (see Section 2.3 for details). The simulations start at 4.0 Ga (or 4 Gyr before present-day) as the evolution during the Hadean Eon (4.5–4.0 Ga) is not considered. Time is running forward from 0–1500 Myr and it represents the Archean Eon (4.0–2.5 Ga). Accordingly, the radioactive elements are initialised with concentrations, which are appropriate for the period.

2.3 Phase changes and composition

The model uses a parameterisation based on mineral physics data (Irifune and Ringwood, 1993; Ono et al., 2001), in which the materials are divided into olivine (ol), pyroxene-garnet (px-gt), TTG (tonalite-trondhjemite-granodiorite), and melt phase systems. Solid-solid phase transitions are assumed within the olivine and pyroxene-garnet phase systems as considered previously in Xie and Tackley (2004); Nakagawa and Tackley (2012). The composition is mapped linearly into the fraction of different phase systems and it can either be in the continuum between end-members, harzburgite (ultramafic and depleted material, 75% ol and 25% px-gt) and basalt (mafic igneous rocks, 100% px-gt), or TTG (the felsic rocks that dominate Archean continental crust (Jahn et al., 1981; Drummond and Defant, 1990; Martin, 1994)). The mantle is initialised with a pyrolytic composition: 80% harzburgite and 20% basalt (Xu et al., 2008), thus being a mixture of 60% olivine and 40% pyroxene-garnet. At a depth of 40 km, a basalt-eclogite phase transition is applied which makes basalt denser than harzburgite by around 190 kg/m³ and has a first-order influence on the lithospheric dynamics (Lourenço et al., 2016). At a depth of 290 km, TTG material undergoes a coesite-stishovite phase transition (Akimoto and Syono, 1969; Akaogi and Navrotsky, 1984; Gerya et al., 2004) and its density increases by 168 kg/m³. Additionally, in the deeper mantle, the phase transition to post-perovskite is considered (e.g. Tackley et al. (2013)). The phase change parameters are given in Table 2. Changes in composition arise from melt-induced differentiation, which is described in the next section.

TABLE 2

TABLE 2. Phase change parameters for olivine, pyroxene-garnet, and TTG systems with surface density at zero pressure ρ_s, density jump across a phase transition Δρ, and Clapeyron slope γ.

2.4 Melting and crustal production

We used the melting model initially developed by Xie and Tackley (2004); Nakagawa et al. (2010), and further extended by Jain et al. (2019) for self-consistent generation of oceanic and continental crust by two-step mantle differentiation. Any newly generated melt is partly intruded at the base of the pre-existing crust as molten material and partly erupted at the surface as solid crust. The mass ratio of erupted to intruded material follows a specified constant value of 30:70 (corresponding to 30% eruption efficiency) based on the previous studies by Rozel et al. (2017); Jain et al. (2019), which have demonstrated the important role of intrusive-dominated magmatism towards producing Archean TTG rocks. In this study, we introduce a solubility function and use solidi that are hydration- and composition- dependent.

2.4.1 Hydration- and composition- dependent solidi

In this study, 8 pressure-dependent solidus temperatures corresponding to 4 compositions (harzburgite, pyrolite, basalt and TTG) and for both hydrated and dry conditions are pre-computed (presented in Supplementary Figure S1A and Supplementary Table S1). The solidus temperature used in the melting routines are interpolated from these pre-computed quantities, instead of being extrapolated from the pyrolite solidus. Specific P-T windows (Moyen, 2011) have been previously used to estimate TTG volumes (Rozel et al., 2017) or to form TTGs (Jain et al., 2019), yet a given range of melt fraction itself was not required. Using interpolated solidi, the TTG windows are matched with a realistic melt fraction as the hydrated solidus for basalt is well controlled and extrapolation-related problems are avoided. Once the end members’ solidi are known, the solidus of material of any composition and hydration is obtained via a three step procedure: (1) calculation of the cell water saturation, and (2) interpolation of end-member solidus in the water saturation space, and (3) interpolation of the solidus in the composition space (if required). The detailed parameterisation is given in Supplementary Appendix S1.

2.5 Water ingassing, partitioning, outgassing and solubility

In our simulations, water is one of the prerequisites for formation of TTG rocks and its concentration has a direct influence on the rheology, when dehydration strengthening is considered. At each time-step, the surface rocks in the top 10 km are saturated with water based on the solubility function (ingassing). Like other quantities such as radionuclides or trace elements, water is partitioned during melting and/or solidification with a partitioning coefficient D_p = 0.01. This ensures that the majority of water present in overheated hydrated rocks (i.e., when the temperature exceeds the solidus) is partitioned into the newly formed melt. When the melt is erupted/intruded, water is also transported along with it, thereby removing water from larger depths (outgassing). Furthermore, to account for the limited water storage capacity of the mantle, a temperature-pressure-composition-dependent saturation limit, i.e., solubility is imposed, above which water is instantaneously transported upwards. Figure 1 depicts the solubility of water (in weight %) of each compositional end-member as a function of temperature and pressure. The bottom row shows the water solubility calculated in hydrous minerals with Perplex (Connolly, 2009). The top row shows the approximated solubility that is employed in our simulations to minimise possible extrapolation artefacts. Note that the water solubility in solids is only non-negligible at low temperature and high pressure. Therefore only cold downwellings can bring water into the transition zone. The analytical functions used to calculate the approximated solubility are described in Supplementary Appendix S2.

FIGURE 1

FIGURE 1. Water solubility maps. Top row: analytical functions used in our models. Bottom row: look-up tables obtained with Perplex (Connolly, 2009).

2.6 Magnesium number

The parameterisation for magnesium number (Mg#) employed here is based on the previously proposed empirical relationships obtained by studying experimental melts coexisting with olivine and orthopyroxene. Lee and Chin (2014) identified protolith P-T of 1–5 GPa and 1400–1750 °C for cratonic peridotites and proposed that building blocks for cratonic material were generated by hot shallow melting. The parameterisation details are given in Supplementary Appendix S3.

3 Results

First, two representative simulations (A13 and B8) are shown (compositional evolution illustrated in Figures 2, 3) to highlight the variety of processes and conditions that can lead to the formation of cratons. Following that, a parameter sensitivity study is presented and the influence of different model parameters elucidated.

FIGURE 2

FIGURE 2. Compositional evolution with time of simulation A13 showing the formation of cratonic roots (grey shaded regions). The individual cells may contain more than one material at a given time, and therefore a threshold of minimum local concentration for each material is employed for their visualisation. Cell-based composition field has the following solid and molten materials: TTG ( ≥60% in tangerine), basalt ( ≥60% in dark purple), harzburgite ( ≥40% in teal, lighter shades represent higher harzburgite content and mantle depletion), TTG-melt ( ≥ 50% in peach-orange), basaltic-melt ( ≥ 50% in sky blue), harzburgitic-melt ( ≥ 50% in white), and TTG-basalt-mix ( ≥ 40% TTG and ≥ 40% basalt in light purple). Superimposed are the contours for water concentration and velocity arrows in the lithosphere.

FIGURE 3

FIGURE 3. Compositional evolution with time of simulation B8 in its later stages showing the formation of cratonic roots (grey shaded regions). Superimposed are the contours for water concentration and velocity arrows in the lithosphere. The reader is referred to Figure 2 caption for explanation of different composition fields.

3.1 Simulation A13

Figure 2 illustrates the evolution of a typical simulation with a relatively high lithospheric strength, hereafter referred to as simulation A13. Early decompression melting events caused by mantle plumes result in widespread crust formation. The residual mantle material becomes depleted and has a higher harzburgite content. The newly generated basaltic melt is partly intruded at the base of the pre-existing crust as molten material and partly erupted at the surface as solid oceanic crust, following an eruption efficiency of 30%. Within each cell, the oceanic crust melts to generate TTG melt, provided that water solubility in basalt exceeds 50% and the necessary P-T conditions for TTG formation Moyen (2011) are met. Following the same eruption efficiency, a fraction of this melt solidifies at the surface as felsic crust, while the rest is intruded at the bottom of felsic crust (if pre-existing) as molten material.

As the primary goal of this work is to understand the processes that lead to craton formation, relevant material properties are continuously tracked and visualised. Since cratonic roots are understood to be viscous, buoyant, depleted, and dehydrated, the cratonic roots in this work are represented with lithosphere older than 750 Myr with a Mg# ≥ 0.905, and water concentration between 10 and 50 ppm (unless otherwise noted). The age of the mantle considered here is based on the last time the material underwent melting. By definition, any material is only considered as cratonic root if it has been around for more than 750 Myr, and therefore cratonic roots only appear during the later stages of the model evolution.

The mantle plumes are the major driver behind the lithospheric dynamics and tectonic regimes observed here. For the first ∼250 Myr (not shown here), the intense plume activity generates crust and recycles it by means of negatively buoyant subducting slabs. The surface rms velocities can reach ≈100–300 cm/yr and the dynamics observed here are similar to a mobile-lid regime (Tackley, 2000; Stein et al., 2004), which is often reported in numerical models as an equivalent of plate tectonics where oceanic lithosphere descends into the mantle.

From 250 Myr onwards, as plumes reach the surface, they push away the lithosphere laterally along both directions (with surface rms velocities ≈0.1–0.2 cm/yr) at a spreading ridge, which can be seen as a discrete zone of divergence in the velocity field in Figure 2. The compression of the pre-existing oceanic and continental crust leads to the tectonic shortening and thickening along horizontal and vertical axes respectively. Around 900 Myr, the cratonic roots start to appear and get thickened and buried deeper with the ongoing compression. The thickened oceanic crust undergoes an eclogite phase transition and is 190 kg/m³ denser than the surrounding residual mantle. This dense lithosphere is gravitationally unstable and its lower layers start to delaminate and drip into the lower mantle. Over time, an accumulation of this recycled oceanic crust can be seen at the core-mantle boundary. During 250–1500 Myr, the ridge continues to spread and no new subduction zones develop. These dynamics persist for a long time and could be considered identical to processes shown by models that exhibit a ridge-only regime (Tackley, 2000; Rozel et al., 2015).

3.2 Simulation B8

Figure 3 illustrates the evolution of a simulation with dehydration strengthening and moderate lithospheric strength, hereafter referred to as simulation B8. Similar to A13, B8 shows two-stage mantle differentiation and the formation of both oceanic and continental crust. The plumes are driving the lithospheric compression and tectonic thickening, however, without forming any subduction zones or ridges. Non-yielding lithosphere is covering the entire surface, akin to a stagnant-lid regime (Nataf and Richter, 1982; Christensen, 1984; Solomatov, 1995). However, the lithosphere is still deforming internally owing to naturally occurring lateral compression with surface rms velocities reaching ≈0.05–0.1 cm/yr. Throughout its evolution, this simulation is classified to exhibit a sluggish-lid regime (Solomatov, 1995; Lenardic, 2018), which has been proposed previously for models where surface velocities are finite and lower than that of the underlying mantle. Around 1000 Myr, the cratonic roots start to appear and by 1300 Myr, the underthrusting of oceanic crust is seen as a result of the ongoing natural lithospheric compression. Furthermore, this model shows that the excess eclogite from the cratonic roots, which is expected as a result of stacking, can be efficiently removed with the delamination and dripping processes. Animated versions of model simulations A13 and B8 are available as Supplementary Material.

3.3 Influence of model parameters

Important, but poorly constrained rheological parameters that significantly affect the mantle and lithosphere dynamics in these models are lithospheric strength and the influence of fluids on the material strength. Therefore, we explored the rheological parameter space that is spanned by three key model parameters: lithospheric strength (denoted by maximum viscosity limiter η_max), plastic yielding (denoted by surface ductile yield stress $σ_{Y}^{0}$ ), and dehydration strengthening (denoted by viscosity offset due to water Δη_w). To independently assess the parameter sensitivity towards cratonic formation and stability, a total of 38 simulations were performed, which can be categorised in three different sets presented in Figure 4. For each simulation considered in the parameter space, Figure 4 highlights whether it forms cratonic roots or not. It also illustrates the shifts in tectonic regimes with time, which have been characterised on the basis of models’ surface rms velocities and the qualitative analysis of the observed dynamics. Initially, the lithospheric strength is varied in conjunction with either plastic yielding (set A, with no dehydration strengthening) or dehydration strengthening (set B, with no plastic yielding). Building upon these results, all three parameters are varied simultaneously (set C). The different tectonic regimes identified in our models are further summarised in Section 4.2 and illustrated in Figure 5. For a selection of cases from the parameter space, Figure 6 shows the temporal evolution of mean mantle temperature, surface root-mean-square (rms) velocities, crustal volume (basaltic and felsic) and TTG mass.

FIGURE 4

FIGURE 4. Temporal evolution of different tectonic regimes for the entire parameter space. Also highlighted are the cases with cratonic roots and regional-scale resurfacing events. A selection of cases (with bold names) are further explored in Figures 5, 6.

FIGURE 5

FIGURE 5. Snapshots of different cases showing composition and viscosity fields to illustrate the different tectonic regimes. Note that, even though different η_max values were used in these cases, only one colorscale for viscosity field is shown.

FIGURE 6

FIGURE 6. (A) Mean mantle temperature, (B) surface rms velocities, (C) total crustal volume (basaltic and felsic), and (D) TTG mass for a selection of cases from the parameter space.

3.3.1 Set A

Simulations with η_max = 10²³ Pa⋅s do not grow a lithosphere strong enough to resist any deformation. There is naturally occurring tectonic shortening and thickening of the crustal material, however, the crust is being continuously generated and recycled into the mantle due to the plume activity. The cold downwellings are efficient in cooling down the mantle and surface rms velocities reach up to ≈100–300 cm/yr for low values of $σ_{Y}^{0}$ (60/120 MPa) (for instance, see curve A1 in Figures 6A,B). Higher values of $σ_{Y}^{0}$ (180/240/300 MPa) push the brittle-ductile transition to greater depths and change the recycling mechanism. This is evident from simulations A1-A5 where the observed tectonic regime shifts from a mobile-lid (Tackley, 2000; Stein et al., 2004) to a plutonic-squishly-lid (Rozel et al., 2017; Jain et al., 2019; Lourenço et al., 2020). With higher $σ_{Y}^{0}$ , the lithosphere can thicken more before it deforms plastically in a ductile manner. While simulations with a mobile-lid regime have negatively buoyant oceanic lithosphere subducting as long, thin, and coherent slabs; the simulations with a plutonic-squishy-lid regime have surface rms velocities ≈10–20 cm/yr and they recycle crust in the form of thick and dense blocks. Once transformed into the denser eclogitic material at a depth of 40 km, the thicker lithosphere becomes gravitationally unstable and starts to drip and delaminate (see Figure 5 for a visual comparison of different tectonic regimes).

Simulations with η_max = 5 × 10²³ Pa⋅s also show crustal shortening and thickening while exhibiting a variety of tectonic regimes depending on the value of $σ_{Y}^{0}$ . However, no stable continents form as most crust is recycled by either global- or regional-scale resurfacing events. For instance, A7 shows a mobility burst with surface rms velocities ≈10–20 cm/yr, during which the entire non-yielding lithosphere resurfaces over a period of 200–300 Myr. Previously, similar dynamics have been reported in models as characteristics of an episodic-lid regime, which has been proposed to be active on Venus (Armann and Tackley, 2012; Turcotte, 1993; Moresi and Solomatov, 1998; Noack et al., 2012). This global resurfacing event can be seen as a decline in mean mantle temperature or a spike in surface rms velocity around 800 Myr in Figures 6A,B. The crustal volume fluctuates highly as the recycled crust melts again to generate new basaltic and felsic crust (Figure 6C).

In cases with η_max = 10²⁴ Pa⋅s, lower $σ_{Y}^{0}$ values (60/120 MPa) result in a thick and a single-plate covering the entire surface with surface rms velocities ≈0.01–0.02 cm/yr (Figure 4). Lithospheric erosion from underneath happens but the lithosphere does not yield or compacts laterally, which are characteristics of a stagnant-lid regime (Nataf and Richter, 1982; Christensen, 1984; Solomatov, 1995). These cases have many downwellings in the initial stages (mobile-lid in the first ∼ 250 Myr) triggered by plumes and the mantle cools down drastically. The colder mantle convects less vigorously and causes less deformation. With higher $σ_{Y}^{0}$ values (180/240/300 MPa), the planet exhibits a ridge-only regime and forms cratonic roots, albeit with occasional global resurfacing events. For instance, with no major downwellings, case A13 shows a steady increase in mean mantle temperature and the surface rms velocity is ≈0.1–0.2 cm/yr (Figures 6A,B).

3.3.2 Set B

Simulations with η_max = 10²³ Pa⋅s have shortening of crustal (TTG+basaltic) material and are in a plutonic-squishy-lid regime. While the denser basaltic-eclogitic crust recycles into the mantle, a high proportion of the buoyant TTG crust remains preserved at the surface (compare B5 curve in panels C) and D) of Figure 6). On the way down, the downwellings release water, which accumulates in the transition zone. Following the imposed saturation limits (see Section 2.5 for details), this water migrates upwards, thereby hydrating the asthenosphere and reducing its viscosity by the offset Δη_w. As a result of tectonic shortening, there is underthrusting of the old and depleted lithosphere beneath the TTG crust. However, these stacked domains have water concentrations reaching up to 5000 ppm and hence are not classified as cratonic roots in this study.

Cases with η_max = 10²⁴ Pa⋅s show tectonic shortening and thickening and are in a sluggish-lid regime. Crustal recycling is minimal and only happens in the form of eclogitic drips (see curve B12 in panels C) and D) of Figure 6). As a result, the overall thicker crust retains most of the water and the underlying upper mantle is drier and stronger when compared to cases with η_max = 10²³ Pa⋅s. There is some stacking of depleted lithosphere, leading to the formation of cratonic roots. However, excess eclogite remains in these depleted domains. The mean mantle temperature continues to rise and surface rms velocities remain low (see curve B12 in panels A) and B) of Figure 6).

Cases with an intermediate value η_max = 5 × 10²³ Pa⋅s show crustal recycling which is somewhat higher than cases with η_max = 10²⁴ Pa⋅s (compare curves B8 and B12 in Figure 6C). These cases (B6-B10) show the successful formation of cratonic roots by underthrusting mechanism while efficiently removing excess eclogite from the cratonic roots.

3.3.3 Set C

Based on the absence of cratonic roots from the relevant simulations of sets A and B, the value of η_max = 10²³ Pa⋅s was dropped from further numerical experiments. All cases start with a regime which is somewhere between plutonic-squishy-lid and mobile-lid: surface rms velocities vary over 2 orders of magnitude and crustal recycling happens both as thick blocks and thin subducting slabs. Post-arrival of the first plumes at the surface, simulations exhibit different tectonic regimes depending on the parameters employed. Most cases with η_max = 5 × 10²³ Pa⋅s (C1-C3) show a stagnant-lid regime at different onset times with a single-lid covering the entire surface. Prior to that, a ridge-only regime operates as long as a high lithosphere-asthenosphere viscosity contrast exists. The higher amount of water in the asthenosphere (brought down with recycled oceanic crust) weakens the material and accentuates this viscosity contrast. For instance, a global resurfacing event shifts the tectonic regime from a ridge-only to a stagnant-lid for case C1 around 380 Myr. A correlation between a sudden drop in mean mantle temperature, a spike in surface rms velocities, and decrease in crustal volume is evident from Figure 6. C4 shows underthrusting and formation of cratonic roots, only if the criterion of maximum water concentration allowed is relaxed from 50 ppm to 200 ppm.

All cases with η_max = 10²⁴ Pa⋅s (C5-C8) show lithospheric compression and thickening. Viscous and depleted cratonic roots also form, once the maximum water concentration allowed is increased to 200 ppm. Irrespective of the $σ_{Y}^{0}$ value employed, the dynamics and tectonics depend more on the value of Δη_w. With a lower value of Δη_w = 0.1, cases C5 and C7 are operating in a ridge-only regime. With a higher value of Δη_w = 0.01, cases C6 and C8 are exhibiting a sluggish-lid regime. What differentiates C6 and C8 from B15 is that with the combination of plasticity and dehydration strengthening, the entire lithosphere is old and viscous as opposed to localised cratonic root formation underneath or around continental margins.

4 Discussion

4.1 Cratonic roots: Defining criteria, formation and stability

As introduced in Section 1, several hypotheses have been proposed to explain the origin of cratons. While the models presented in this study lack the many complexities of a natural system, they are able to form cratonic roots through processes, which could be considered as a combination of different formation mechanisms. In our global models, mantle plumes cause decompression melting and generate continental crust by a two-step mantle differentiation. At the same time, the residual material left in the underlying lithosphere is melt-depleted and dehydrated, but it is not as such considered as cratonic lithosphere in this study. Here, one of the defining criteria for cratonic roots is that the lithosphere has to show resistance against possible erosion or deformation by the convecting mantle for more than 750 Myr. i.e., the material does not go undergo any further melting. In these models, mantle plumes also cause lithospheric compression, which in turn thickens the lithosphere and facilitates the recycling of gravitationally unstable eclogitic material from the deeper layers. Such cases show processes similar to underthrusting of pre-existing oceanic lithosphere beneath the crust. These processes have been proposed as a formation mechanism (Helmstaedt and Schulze, 1989; Canil, 2004, 2008; Simon et al., 2007; Pearson and Wittig, 2008) and might have also been responsible for the pressure increase experienced by cratonic peridotites Lee and Chin (2014).

Cratonic roots are considered to be made of residual peridotites (harzburgite in this study), which are left behind as a result of melting of igneous protoliths (pyrolite in this study). Progressive melt extraction increases the magnesium number in these residues and makes them compositionally buoyant by decreasing their densities (Jordan, 1979; Boyd, 1987, 1989; Lee, 2003; Schutt and Lesher, 2006). Assumed to be reflective of the secular cooling of the Earth’s mantle, Archean peridotites have higher Mg# (0.91–0.94) than Proterozoic and Phanerozoic peridotites (0.90–0.92) (Griffin et al., 2003; Lee et al., 2011). Here, the Mg# evolution during protolith to residual peridotite melting is parameterised using the empirical relationships proposed by Lee and Chin (2014). Accordingly, in this study, another criterion for classifying depleted lithosphere as cratonic roots, is that it must have a Mg# of at least 0.905. However, the density reduction of the residual material is not explicitly incorporated (e.g., Lee, 2003). The density decreases implicitly with melt depletion, i.e., with increase in olivine content as a result of the phase changes considered in Section 2.3. As olivine is less dense than pyroxene-garnet below the eclogitic phase transition depth of 40 km, the residual peridotite or harzburgite (75% olivine, 25% pyroxene-garnet) becomes less dense than pyrolite (60% olivine, 40% pyroxene-garnet) below that depth by 29 kg/m³.

Being intrinsically stronger than the convecting mantle, cratonic roots are understood to be dehydrated with low water concentrations. The amount of water bound in peridotites varies with their age and tectonic environment. The clinopyroxenes from Kaapvaal peridotite xenoliths have water concentrations reaching up to 400 ppm (Grant et al., 2007; Peslier, 2010), corresponding to 2–18 ppm H₂O content in olivine. The H₂O content of olivine in the MORB mantle source is estimated to be around 10–30 ppm (Hirschmann, 2006). In comparison, clinopyroxenes from the Archean North China Archean and Proterozoic Colorado Plateau are drier and wetter respectively (Li et al., 2008; Yang et al., 2008). As the most depleted mantle material (harzburgite) in these simulations only has 75% olivine (the rest being pyroxene-garnet), a water concentration of 10–50 ppm is considered as another defining criteria for cratonic roots. For some cases, this upper limit of 50 ppm is relaxed to 100 or 200 ppm.

As the models presented in this study span only 1500 Myr of evolution under Archean conditions, they are not suited to give insights on the long-term stability of cratons. In all these models, the previously formed cratonic roots (for instance, shown in Figures 2, 3, 5) do not resist any further deformation caused by ongoing lithospheric compression as subsequent mantle plumes spread at the surface. They continue to thicken and eventually recycle by destabilisation. At the same time, there are new regions (grid cells) that satisfy the above-mentioned defining criteria (lithosphere older than 750 Myr with water concentration between 10 and 50 ppm and magnesium number of at least 0.905) and become cratonic roots. Considering how dynamic and short-lived both these processes of formation and destruction are, we did not attempt to quantify the lateral or vertical extents of these roots.

Our models are missing natural processes that could make these cratonic roots even more chemically buoyant and intrinsically strong. Previous models have shown that in addition to the positive chemical buoyancy (Lenardic and Moresi, 1999; Lenardic et al., 2003; Sleep, 2003), an effective viscosity contrast of 20–1000 between the cratonic roots and the surrounding mantle is necessary to keep the cratonic roots from recycling (Sleep, 2003; O’Neill et al., 2008; Wang et al., 2014). In this study, the depleted material within the cratonic roots is less dense than the surrounding non-depleted material by only 29 kg/m³ below the eclogitisation depth of 40 km. In future work, these phase systems (see Section 2.3 for details) could be adapted to allow for higher degrees of mantle depletion, i.e., material with 100% olivine. Then the corresponding density contrast between the depleted cratonic roots and the non-depleted surrounding material would increase to 76 kg/m³.

4.2 Tectonic regimes during the Archean and their implications

For the early Earth, a variety of tectonic regimes (post magma ocean solidification (Abe, 1993a,b, 1997; Elkins-Tanton, 2012)) have been proposed albeit without any consensus. And when it comes to narrowing down the onset of subduction-driven plate tectonics, opinions vary widely with a multitude of studies proposing its inception anytime between the Hadean Eon (4.5–4.0 Ga) and the Neoproterozoic Era (1.0–0.54Ga) (Korenaga, 2013; Cawood et al., 2018; Palin and Santosh, 2020). While the main objective of this study was to model self-consistent craton formation, the simulations presented here also offer us a diversity of tectonic regimes as summarised below:

• Plutonic-squishy-lid: In our models, it is characterised by having a thin lithosphere with small and ephemeral plates made by magmatic intrusion, surface rms velocities ≈10–20 cm/yr, and crustal recycling in the form of thick and dense blocks of delaminated and dripping basaltic-eclogitic material. Generated by plutonism, this regime has been suggested to be applicable for early Archean Earth and present-day Venus using global models (Rozel et al., 2017; Jain et al., 2019; Lourenço et al., 2020). It has similarities to the plume-lid tectonics regime described by Sizova et al. (2010); Gerya et al. (2015); Fischer and Gerya (2016a,b) in their regional models.

• Mobile-lid: Similar to plate tectonics, models with this regime have negatively buoyant oceanic lithosphere going down as thin (when compared to thick blocks in a plutonic-squishy-lid regime) and coherent slabs due to their negative buoyancy (Tackley, 2000; Stein et al., 2004; Bédard, 2018) and surface rms velocities can reach ≈100–300 cm/yr in the first 200 Myr. Contrary to our results, Capitanio et al. (2022) showed the formation of cratonic roots in their models with mobile-lid regime by employing melt-depletion dependent rheological stiffening. They observed that low lithospheric strength favours an increase in mobility with widespread rifting and downwellings. The subsequent higher melting volumes and melt-depletion stabilises the lithosphere.

• Ridge-only: In models with this regime, a discrete zone of divergence forms in the lithosphere without the presence of any subduction zones and surface rms velocities are ≈0.1–0.2 cm/yr. Previously, this regime has been suggested to operate on early Earth by Tackley (2000); Rozel et al. (2015) where they employed a depth-dependent yield stress in their models. At the spreading ridge, the lithosphere is pushed laterally along both directions in a 2D setup, resulting in its compaction and thickening. These models are successful in forming cratonic roots, however, these roots may not achieve long-term stability. In some cases (A7, A8, A14, A15), the lower and denser layers of the thickened lithosphere (several 100 kms), which are gravitationally unstable, delaminate along with the entrained cratonic roots.

• Episodic-lid: Classically, models with this regime are characterised by an unyielding lithosphere that experience mobility bursts and as a result, the entire lithosphere is resurfaced in a very short amount of time (Armann and Tackley, 2012; Turcotte, 1993; Moresi and Solomatov, 1998; Noack et al., 2012). In our models, the observed surface rms velocities reach up to ≈10–20 cm/yr during these bursts and the duration of these global resurfacing events is a magnitude higher than what has been reported before (Lourenço et al., 2016). These overturns happen following the ridge-only regime, as opposed to the classical definition, where stagnant-lid precedes them.

• Stagnant-lid: In this regime, the entire surface is covered by a single-plate as the thick lithosphere does not yield (Nataf and Richter, 1982; Christensen, 1984; Solomatov, 1995). With surface rms velocities ≈0.01–0.02 cm/yr, lithospheric erosion from underneath is possible, but there is no lateral compression or tectonic thickening. As cratons are understood to be regions of crustal basement localised underneath continental crust, this global single-plate lithosphere is not considered cratonic in this study even when it meets all the criteria outlined in Section 4.1.

• Sluggish-lid: Models exhibiting this regime have surface rms velocities reaching ≈0.05–0.1 cm/yr (not stagnant), which are lower than that of the underlying convecting mantle (Solomatov, 1995; Lenardic, 2018). This regime is similar to a stagnant-lid in terms of having a non-yielding lithosphere, however, this lithosphere is being continuously deformed internally owing to naturally occurring lateral compression. These models are also successful in forming cratonic roots and show underthrusting of oceanic crust. Similar stacking of subducted oceanic lithosphere has been suggested to have taken place in the Kimberley block of the Kaapvaal craton during the Neoarchean (Pearson and Wittig, 2008).

Based on these results of coupled core-mantle-crust evolution on global scale, we propose the following shifts in tectonic regimes during Earth’s early history:

• In the Hadean, mantle plumes must have been more vigorous due to a higher core temperature and responsible for the majority of the crust production and recycling. Assuming that most of the magmatism on early Earth was dominated by plutonism (Rozel et al., 2017; Jain et al., 2019; Lourenço et al., 2020) similar to present-day (Crisp, 1984), plutonic-squishy-lid on a global scale or plume-lid tectonics on a regional scale might have been operational.

• During the Archean, the rapid cooling of the core (e.g., Nakagawa and Tackley, 2004) and the resulting less efficient crustal recycling might have ensured that the majority of the water is retained in the crust itself. The underlying mantle would have been melt-depleted, strong, and relatively dry. Mantle plumes might have facilitated the lateral compaction and the subsequent tectonic thickening of the lithosphere. During this compaction, the stacking of oceanic lithosphere might have led to the formation of stable cratonic roots. And the denser eclogitic material might have then delaminated and dripped into the lower mantle as shown in our models. These dynamics would have been possible under a ridge-only or a sluggish-lid regime.

4.3 Felsic crustal volume, recycling and TTG proportions

While the growth rate, composition, and preservation of continental crust during Earth’s history remain topics of active discussion in the scientific community (e.g., Cawood and Hawkesworth, 2018; Korenaga, 2018; Hawkesworth et al., 2020), this study shows the evolution of felsic crust during the Archean Eon, which can be compared against geological observations and previous models. The combined volume of basaltic and felsic crust (Figure 6C) and the mass of produced TTG (Figure 6D) presented here are computed after correcting for the geometry used in these simulations (see Supplementary Appendix S4 for details). While a big proportion of the felsic crust (TTG rocks) is produced in the first ∼250 Myr of the Archean with the arrival of the first mantle plumes, the production of TTGs is continuous for the entire 1500 Myr evolution (Figure 6D), which agrees with the observations of TTGs from throughout the Archean (Moyen and Laurent, 2018). Also visualised is the remaining mass of TTG in the top 100 km, which tells us about the intensity of recycling: very high for cases A1 (mobile-lid regime), A7 (ridge-only and episodic-lid regimes) and minimal for B8, B12 (both sluggish-lid regime).

Recent work by Guo and Korenaga (2020) also advocates for rapid continental crust growth during the early Earth. The authors modelled argon degassing with thermal evolution of Earth and estimated that more than 80% of the mass of present-day continental crust was already reached during the Archean. This would correspond to a production of 5 × 10²² kg of felsic rocks during the Archean (Korenaga, 2021). Supplementary Table S4 lists the cumulative mass of TTG produced, which varies between 1.38 × 10²² - 4.55 × 10²² kg depending on the simulation. The tectonic regime has a direct influence on the amount of TTG produced. In cases with a mobile-lid regime (A1 for example), the recycled oceanic lithosphere partially melts, leading to higher production of TTG rocks. In cases with a sluggish-lid regime (B12 for example), crustal recycling is minimal and it generates a lower amount of TTGs.

Moyen (2011) classified the Archean TTGs into three different types based on the pressure of their formation settings: low-pressure (10–12 kbar), medium-pressure (ca. 15 kbar), and high-pressure (20 kbar or higher) TTGs, which account for 20%, 60%, and 20% of the sodic TTGs, respectively. Using global models, Rozel et al. (2017) estimated TTG volumes and proposed that a plutonism-dominated plutonic-squishy-lid regime would be able to reproduce this observed proportions of different TTG rocks. However, in a follow-up study, Jain et al. (2019) reported relative proportions of 85%, 14.5%, and 0.5% based on their models of TTG production. In their work, most of the TTGs were produced at the tip of the deformation fronts when plumes spread laterally at the surface. Therefore, these models lacked certain tectonic processes that might have been responsible for the production of high-pressure TTGs. In the present study, the best case scenario for these relative proportions is 65%, 20%, and 15%, which occurs for simulations operating with a sluggish-lid regime (see Supplementary Table S4 for all simulations). These simulations show lithospheric compression and tectonic thickening forming cratonic roots without any subduction zones. This allows the basaltic crust to get buried to higher pressures and results in an increased proportion of high-pressure TTGs. Also, models with plutonic-squishy-lid or ridge-only regime shift the relative proportions of produced TTG rocks towards higher pressures as basaltic crust gets to thicken over time. Models operating with mobile-lid regime, however, generate mostly low-pressure TTGs as the basaltic crust gets recycled before it can grow thick. In their regional-scale models, Gunawardana et al. (2020) also observed the important role of lithospheric dripping and a weak lithosphere towards reproducing the relative proportions of low-pressure and medium-pressure TTGs from the Archean rock record.

4.4 Model limitations and possible future improvements

Although these models offer us insights about the formation of stable continents during early Earth, they have certain shortcomings, which should be addressed in future studies. First, these models do not consider grain-size dependent rheology (Bercovici and Ricard, 2005, 2014; Ricard and Bercovici, 2009; Rozel et al., 2011), which might help with the stabilisation of cratonic roots (King, 2005). Hall and Parmentier (2003) showed a viscosity increase in regions of low convective stresses due to grain size evolution. Grains might grow bigger in the low strain environments of cratonic roots and might make them stronger. Second, cratonic roots in our models are not as compositionally buoyant as what is expected from geological observations. Third, Mondal and Korenaga (2018); Korenaga (2021) argue that delamination of the dense, thick, and gravitationally unstable eclogitic crust in our models may not occur if a stronger rheology for eclogite is considered (Jin et al., 2001). Fourth, mantle plumes might have a limited impact on the lithosphere dynamics in 3-D models with opposing effects. On one hand, this might reduce the lithospheric compression and tectonic thickening responsible for the formation of cratonic roots. On the other hand, the stability of already formed cratonic lithosphere might improve. Fifth, increasing the resolution in these global simulations might allow us to reproduce geological features such as dome and keel structures (Hickman, 2004; Van Kranendonk et al., 2004; Van Kranendonk, 2011). Finally, future work is needed to self-consistently model the regime transition from pre-plate tectonics active during early Earth to subduction-driven plate tectonics.

5 Conclusions

We have presented global mantle convection models that generate viscous, melt-depleted (with Mg# ≥ 0.905), dehydrated (with 10–200 ppm water content), and stable (older than 750 Myr) cratonic roots in a self-consistent manner. We considered the protolith P-T conditions of cratonic peridotites previously identified by Lee and Chin (2014) and modelled the Mg# evolution in the mantle. In our models, formation of cratonic roots happens as a result of a combination of different tectonic processes and for a range of rheological parameters employed in this study. Initially, decompression melting caused by mantle plumes generates oceanic crust and Archean continental crust (considered as TTG rocks here). Concurrently, the residual material left behind in the underlying lithosphere is melt-depleted and dehydrated. Over time, this lithosphere gets laterally compacted and thickened and evolves into long-lived and strong cratonic roots. Our simulations show a variety of tectonic regimes and transitions between these regimes over time are also observed. In cases with a ridge-only regime, there is no active subduction and the lithosphere is laterally compressed on both sides of a long-lasting spreading ridge. Cases with a sluggish-lid regime also show naturally occurring lateral compression and thickening of the lithosphere, thereby allowing the cratonic roots to mature over 100s of millions of years. Most simulations also show underthrusting of oceanic lithosphere beneath the pre-existing crust, which has been suggested by Pearson and Wittig (2008) to explain the evolution of the Kaapvaal craton during the Neoarchean. The denser eclogitic material within these underthrusted domains is gravitationally unstable and is recycled into the lower mantle by delamination and dripping processes. These simulations, therefore, offer a mechanism to explain the rare presence of eclogites among the mantle xenoliths in kimberlites (Schulze, 1989). We envisage that early Earth exhibited different tectonic regimes during its evolution. Considering the initial vigorous nature of mantle plumes, a plutonic-squishy-lid regime was active during the Hadean Eon and led to intense production and recycling of both oceanic and continental crust. Sometime during the Archean Eon, the secular cooling of the mantle and the increase in lithospheric rigidity caused a regime transition to a ridge-only or sluggish-lid, which offered conditions suitable for the formation of cratonic roots.

Data availability statement

The convection code StagYY is the property of Paul J. Tackley and ETH Zurich and is available for collaborative studies from Paul J. Tackley (paul.tackley@erdw.ethz.ch). The scripts and raw datasets used to produce Figures 1, 6 for this study can be downloaded from Zenodo with the following DOI: 10.5281/zenodo.6630970.

Author contributions

JvH conceptualised the study and outlined the research direction. JvH, CJ, and AR designed the set of numerical simulations. CJ ran the simulations and did post-processing. CJ and AR produced the figures. AR and AM-C developed the water/melting parameterisation and AR coded it in StagYY. EC helped with the Mg number parameterisation and CJ implemented it. All authors contributed to the article and approved the submitted version.

Funding

CJ and JvH were supported by the Natural Environment Research Council under the grant NE/M000281/1. CJ is also receiving funding from the ERC Synergy Grant 856555 for the project: Monitoring Earth Evolution through Time (MEET). AR is funded by ETH Zürich.

Acknowledgments

We would like to thank James Connolly for his advice and for providing us with the water solubility. We would also like to thank Peter Cawood and Adam Beall for their constructive feedback during the review process that helped improve the original manuscript, and Simone Pilia for his editorial work. This work has made use of the facilities of ARCHER UK National Supercomputing Service (https://www.archer.ac.uk) and the Hamilton HPC Service of Durham University. The scientific colour map lapaz (Crameri, 2021) is used in this study to prevent visual distortion of the data and exclusion of readers with colour-vision deficiencies (Crameri et al., 2020). The StagPy library (Morison et al., 2021) and Paraview were used to process and visualise StagYY output data respectively.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.966397/full#supplementary-material

Supplementary Presentation S1 | Parameterisation details for interpolated solidi, water solubility and magnesium number.

Supplementary Video S1 | Evolution of simulation A13.

Supplementary Video S2 | Evolution of simulation B8.

References

Abe, Y. (1993b). Physical state of the very early Earth. LITHOS 30, 223–235. doi:10.1016/0024-4937(93)90037-d