Calibrating the Yield Strength of Archean Lithosphere Based on the Volume of Tonalite-Trondhjemite-Granodiorite Crust

Introduction

Continental crust has played a key role in the evolution of the Earth system (e.g., Cawood et al., 2013). While its production is controlled by plate tectonics on the modern Earth (e.g., Rudnick and Gao, 2003), the tectonic processes that formed felsic crust and continents on the early Earth remain debated and range from plate tectonic to non-plate tectonic settings (e.g., Martin and Moyen, 2002; Bedard, 2003; Foley et al., 2003; Rapp et al., 2003; Harrison, 2009; Smithies et al., 2009; Chowdhury et al., 2017; Johnson et al., 2017; Moyen and Laurent, 2018; Rosas and Korenaga, 2018; Capitanio et al., 2019; Smithies et al., 2019; Chowdhury et al., 2020; O’Neill et al., 2020; Turner et al., 2020). This uncertainty partly stems from the significant gap in the preserved rock record and therefore, in our limited knowledge about the proportion and nature of the early continental and oceanic lithospheres (Cawood et al., 2013; Cawood et al., 2018; Lenardic, 2018; Stern, 2018). Numerical modeling of the early Earth processes and ground truthing its results with the limited rock record is, therefore, an important way to understand the early Earth processes (Rozel et al., 2017; Capitanio et al., 2019; Gerya, 2019; Chowdhury et al., 2020).

One approach to understand the Archean continental crust forming processes is by constraining the physical properties of the lithosphere(s) prevalent at that time. Previous research has shown that whether or not a terrestrial planet manifests some form of plate tectonics will depend on the lithospheric yield stress, i.e., the plastic strength of the crustal and the underlying lithospheric mantle rocks (Moresi and Solomatov, 1998; Trompert and Hansen, 1998; Tackley, 2000; O’Neill and Lenardic, 2007; van Heck and Tackley, 2008; Korenaga 2010a; Korenaga, 2010b; Van Heck and Tackley, 2011). Based on the magnitude of this parameter, a planet can exist in stagnant-lid, episodic-lid or mobile-lid tectonic mode (Moresi and Solomatov, 1998), since changing the yield stress value affects the mantle convection pattern via affecting the viscosity contrast between the lithosphere and convecting mantle. For example, present-day continental crust production mainly occurs at subduction and orogenic settings within a plate tectonic (mobile-lid) regime where average lithospheric yield stresses vary within >100–150 MPa (Kohlstead et al., 1995; Gurnis, 1988; Lowman and Jarvis, 1999; Tackley, 2000). However, this parameter remains poorly constrained for the Archean lithosphere.

Here, we perform high resolution 2D numerical modeling of the lithospheric dynamics under hotter mantle conditions appropriate for the early-to-mid Archean (>3 Ga) time-period, for a range of yield stress values varying between 50 and 250 MPa. We evaluate the models by calibrating the tempo of granitoid formation as obtained from them with the preserved geological record, and thereby, seek to constrain the yield stress values of Archean lithosphere and the consequent tectonic settings responsible for the formation of early felsic crust. Changes in the lithospheric processes due to different yield stress values are likely to affect the pressure-temperature (P-T) conditions of tectonic environments where metamorphism and partial melting processes operate. This implies that yield stress values are implicitly linked to the formation conditions and volumes of Archean felsic crust, which can be independently constrained from other approaches, thus, can be used to trace it back for the ancient lithosphere.

Methodology

Numerical Methods

We model 2D-mantle convection at high mantle potential temperatures (T_p; ∼1,475–1,600°C) inferred for the Archean (Herzberg et al., 2007; Herzberg et al., 2010; Korenaga, 2013; Ganne and Feng, 2017; Aulbach and Arndt, 2019) using the numerical code Underworld 2 (Moresi et al., 2007). This code is based on the finite element method and is coupled with the particle-in-cell approach (Evans and Harlow, 1957). It involves a Eulerian grid with the moving Lagrangian particles embedded in a finite element (Moresi et al., 2007). Each particle represents a specific solid or molten rock-type (termed as rock-particle hereafter) and carries the information about its physical properties like, P, T, density, etc., which can be tracked over the entire model duration. The mass, momentum, and energy conservation equations are solved under incompressible conditions to find the pressure, velocity, and thermal evolution of the mantle convection model. Physical properties such as density and viscosity associated with the Earth’s interior are mapped to these equations through particle indexing. The mass, momentum, and energy conservation equations are:

where u, P, $\eta$, T, H, ρ, g, and t represent velocity, pressure, viscosity, temperature, internal heat, density, gravitational acceleration, and time, respectively, (see Table 1 for the parameter values).

TABLE 1. Physical properties used in numerical modeling.

The temperature dependent diffusion creep rheology is applied using the Arrhenius equation (Rozel et al., 2017):

where ${\eta }_o$ is the reference viscosity at T = 1. Here, the parameters E, V, and d represent the non-dimensional activation energy, activation volume, and depth coordinate of the model, respectively. Following some recent studies (e.g., Sizova et al., 2015; Rozel et al., 2017), the activation energy is applied in such a way that temperature dependent rheology creates a stagnant lid mode of convection in the model (Moresi and Solomatov, 1998). A lithospheric weakening mechanism is implemented using the yield stress $\left({\sigma }_Y\right)$ parameter by assuming the ductile (plastic strength) behavior of the lithospheric rocks (Moresi and Solomatov, 1998; Sizova et al., 2015; Rozel et al., 2017). This is the only weakening mechanism that we imposed in the model. The yield stress equation is:

where ${\sigma }_Y^o$ is Earth’s surface yield stress and ${{\sigma }^{\prime }}_Y$ is the depth dependent yield stress gradient. When the convective stresses exceed the yield stress, we choose the yield stress value to calculate the viscosity. The yielding viscosity is ${\eta }_Y= {\sigma }_Y/2\dot{ϵ}$, where $\dot{ϵ}$ is the second invariant of the strain rate tensor. The effective viscosity is then calculated by ${\eta }_e=\left[{\eta }_T^{-1}+{\eta }_Y^{-1}\right]{ }^{-1}$. We have systematically varied the surface yield stress among the values of 50, 100, 150, 200, and 250 MPa, covering a range of weaker to stronger lithospheric strength profiles. These values are proximal to the yield stress compatible with plate tectonics (Moresi and Solomatov, 1998; Tackley, 2000a; Van Heck and Tackley, 2011).

The dimensional values of the above non-dimensional parameters are retrieved using the relationships:

following scaling methods of Moresi and Solomatov (1998), Korenaga (2010b), and Van Heck and Tackley (2011). The primes represent the non-dimensional form of the variable and will be dropped here onwards, as we will only be dealing with non-dimensional symbols.

Partial Melting

To model felsic and mafic crust production, we implemented partial melting of the mantle peridotite and basaltic crust, including melt extraction, residue formation and eclogitization. For any rock-particle, melting initiates whenever its temperature exceeds its solidus. The adiabatic component of the temperature is added to the temperature solution ($T$) of Eq. 3 for the melt calculation. The P-T dependent peridotite melting model of Sizova et al. (2015) is used for the mantle melting and melt extraction processes to account for the realistic scenario of mantle melting over a range of temperature between its wet and dry solidus depending on the extent of hydration. For the basaltic crust, the volumetric degree of melting, i.e., melt-fraction $\left(M_o\right)$ is considered as a linear function of P-T and is calculated using the equation (Gerya, 2010; Sizova et al., 2015; Chowdhury et al., 2017; Chowdhury et al., 2020):

where T_solidus and T_liquidus represent solidus and liquidus temperatures (Figure 1; also see Appendix) and are taken from the previous studies (e.g., Sizova et al., 2015; Rozel et al., 2017; Chowdhury et al., 2020). The melting of felsic crust has not been considered here.

FIGURE 1. (A) The initial compositional setup used to reach a state of stable mantle convection. The geotherm corresponds to a mantle T_p of 1,600°C and mantle adiabat of 0.5°C/km. The thickness of basaltic crust is defined by a Moho-T of 600°C. (B,C) The compositional and thermal profile derived after a stable mantle convection pattern is established. The lithological particles in (B) has been redefined using the thermal profile shown in (C). This stage is considered as the start of model evolution for calculating the volume of tonalite-trondhjemite-granodiorite crust.

To extract melt from partially-molten lithologies, we define four threshold parameters (M1, M2, M3, and M4) following the algorithm of Sizova et al. (2015) for each rock-particle. M1 represents the non-extractable melt volume that always remains in a partially molten rock-particle, while M2 represents the minimum melt volume required for the melt-extraction to take place. The extracted melt volume is assumed to instantaneously leave the melting region and move to the surface (Schmeling, 2006). M3 is the maximum amount of melt that can be sustained in a crystalizing melt particle without melt fractionation. M4 is the maximum amount of cumulative melt that can be extracted from a rock-particle. Once a rock-particle melts beyond M4, it is considered infertile for further melting. All the threshold parameters are assigned distinct values according to the rock-type (mantle peridotite or basalt; see Appendix A1). The melt volume borne by a rock-particle at any instant of time $\left(M\right)$ during the model evolution is given by: $M= M_o-{\displaystyle {\sum }_n{M}_{\text{ext}}}$ where, M_o is given Eq. 6, $M_{\text{ext}}$ is the amount of melt extracted in a single time-step and n is the number of times melt has been extracted from that rock-particle (${\displaystyle {\sum }_n{M}_{\text{ext}}}$ gives the total melt-fraction extracted from that rock-particle). Whenever M exceeds M2 for a rock-particle, the extracted melt fraction $\left(M_{\text{ext}}\right)$ is calculated as M − M1. This value is then added to the extracted cumulative volume of extracted melt for that rock-particle, which is then compared to M4 to check whether that rock-particle has reached its final limit for melt-extraction or not. The prefix melt-bearing is attached to the lithological name of a rock-particle (e.g., melt-bearing peridotite) whenever it carries some melt (M > 0). Once the cumulative melt-volume reaches M4, the melt-bearing rock-particle is converted to an equivalent depleted lithology like, from melt-bearing peridotite to depleted peridotite.

Density Calculations

The temperature dependent density (ρ) for different rock types is calculated using: ${\rho }_o\left[1-\alpha \Delta T\right]$ where the ${\rho }_o$ is the reference density at 1 MPa and 298 K, $\alpha$ is the thermal expansivity, and $\Delta T$ is the temperature drop across the model (Table 1). Partial melting and melt extraction also affect the density of the resultant rock types. The effective density of a melt-bearing rock type is calculated using (Gerya, 2010):

While ρ_melt represents the temperature dependent density of the melt volume (M), ρ_solid is the density of the solid part of the melt-bearing rock-types - called residue - that includes the effect of temperature, melt-extraction and eclogitization. It is calculated using:

where ${\text{coef}}_{\text{res}}$ and ${\text{coef}}_{\text{phase}}$ represent the coefficients of density changes due to melting and eclogitization, respectively.

The density of solid peridotite decreases due to melt depletion (e.g., Poudjom Djomani et al., 2001), which we have implemented using the relation: ${\text{coef}}_{\text{res}}=-0.04{\displaystyle {\sum }_n{M}_{\text{ext}}}$ (Sizova et al., 2015). This corresponds to a decrease of ∼1% in the density of solid residual peridotite for every 25% of extracted melt volume. For depleted peridotite (i.e., M = M4), ${\rho }_{\text{solid}}$ is calculated at a fixed ${\displaystyle {\sum }_n{M}_{\text{ext}}}$ of 40%. No phase transitions are considered for the mantle peridotite (coef_phase = 0). During the melting of a metabasalt at amphibolite facies, the residue becomes rich in garnet, which makes it denser than the unmetamorphosed basalt (Rapp et al., 2003; Zhang et al., 2013). This effect is modeled using the following relation: ${\text{coef}}_{\text{res}}=0.4{\displaystyle {\sum }_n{M}_{\text{ext}}}$ (Sizova et al., 2015). The effect of phase transitions, i.e., eclogitization on the basaltic density is implemented via the term ${\text{coef}}_{\text{ecl}}$. Following the experimental results (Ito and Kennedy, 1971), the values of this coefficient are varied from 0 to 16% within the P-T dependent garnet-in line [P (in kbar) = 0.014 T (in °C) − 5.4] and plagioclase-out line [P (in kbar) = 0.024 T (in °C) + 4] (Sizova et al., 2015), which approximates the conversion of basalt to eclogite through garnet-granulite. It must be noted that eclogitization is not necessarily dependent of melting and therefore, its effect on density is calculated for all the basaltic rock-particles (both solid and melt-bearing).

Initial Setup and Boundary Conditions

Our 2D model domain is created with a 64 × 192 resolution Eulerian mesh and with 20 Lagrangian particles per cell. The model dimension is 660 km × 1,980 km (Figure 1A). A lithosphere with thick, hydrous basaltic crust is used in the models (Figure 1A) in accordance to the recent geological findings (e.g., Herzberg et al., 2010; Herzberg and Rudnick, 2012; Kamber, 2015; Johnson et al., 2017). The crustal thickness is set by the Moho temperature (Moho-T). Whenever a rock-particle exists within the upper 60 km of the model domain and its temperature is less than Moho-T, it is defined as basaltic crust. We varied the crustal thickness using three different values of Moho-T (500, 600, and 700°C). The relevant average basaltic crustal thickness for the model with mantle T_p 1,600°C and yield stress of 50 MPa case is presented in Appendix Figure A1. These different Moho-T values also simulate the effect of relatively colder and hotter crustal geotherms. Following the studies of Rey et al. (2014) and Sizova et al. (2015), the lithospheric thickness is thermally defined by 0°C at the surface and a temperature at a fixed depth of 100 km depending on the mantle T_p values. We further increased the mantle T_p value by 125, 200, and 250°C above its present-day value of 1,350°C in successive steps, which are consistent with the thermal conditions of Archean mantle (Herzberg et al., 2010; Korenaga, 2013). The ambient mantle temperature linearly increases following an adiabatic geotherm of ∼0.5°C/km below 150 km depth. Above 150 km, the mantle adiabat is linearly linked to the lithospheric geotherm over the depth interval of 100–150 km (Figure 1A). The top boundary is kept constant at 0°C while the bottom boundary is kept at constant temperature guided by the mantle T_p and the mantle adiabatic geothermal gradient (see Table 1). Free slip boundary conditions are imposed at the top and bottom of the model, while side walls are given periodic boundary condition. We have run all the models for a finite time at the beginning to reach a steady-state mantle convection (Figure 1C). Once this is achieved, the lithological layers are redefined according to the established steady-state thermal profile (Figure 1B). The stage is considered as the start of the model evolution.

Constraints From the Continental Rock Record

To evaluate our models, we predict the 1) relative abundances of various felsic igneous rocks and 2) the volume of continental crust that would have been formed during the early Archean (>3 Ga), and compare them against the volumes inferred from the rock record.

Comparison with the rock record of Tonalite-Trondhjemite-Granodiorites

To evaluate our models, we need information about the relative abundances of the different types of Tonalite-Trondhjemite-Granodiorite (TTG) - the dominant felsic lithologies that constitute Archean crust (reviewed in Moyen, 2011; Moyen and Laurent, 2018). Extant experimental and petrological studies show that the Archean TTGs mainly formed through the partial melting of hydrated metabasalts (Foley et al., 2003; Rapp et al., 2003; Moyen and Stevens, 2006; Qian and Hermann, 2013; Palin et al., 2016; Johnson et al., 2017). Our implemented algorithm of basalt melting is, thus, appropriate to model TTG formation. The trace element composition of TTGs further suggest that their formation occurred at three different pressure ranges: (Moyen, 2011): 1) at lower pressures (<1 GPa) where garnet-absent amphibolite residues are stable, called low-pressure (LP-) TTGs; 2) at intermediate pressures (∼1–2 GPa) where both garnet and plagioclase are present in the residue, called medium-pressure (MP-) TTGs; and 3) at high pressures (>2 GPa) where garnet and rutile are stable in the residue but not plagioclase, called high-pressure (HP-) TTGs. The bracketing pressure ranges, however, may vary depending on the source rock composition; for example, MP-TTGs can form at P < 1 GPa (∼0.7–0.8 GPa) from low-MgO enriched metabasalts (Johnson et al., 2017).

In this study, we focus on the formation of LP- and MP-TTGs as they dominate the >3 Ga rock record (cf. Moyen and Laurent, 2018), which is the time period of interest. We modeled partial melting of hydrous basaltic crust within a pressure range of 0.7–2.5 GPa, corresponding to a maximum depth of ∼75 km (Figure 2). Notwithstanding the ambiguities of nomenclature, this pressure range overlaps with the overall physical conditions of LP- and MP-TTGs formation during the Archean (Figure 2) (Palin et al., 2016; Johnson et al., 2017; Moyen and Laurent, 2018). Regarding the relative abundances, MP-TTGs constitute nearly 60% of the total TTG volume whereas LP- and HP-TTGs contribute equally toward the remaining 40% (Moyen, 2011; Moyen and Laurent, 2018). Since we focus on the LP- and MP-TTGs, we renormalize their proportions and obtain a ratio of 1:3 (∼25%:75%) between their abundances. Whenever a basaltic rock-particle melts, the resultant melt is classified as LP-/MP-/HP-TTG depending on in whose domain the melting took place (Figure 2). Melt is extracted once the melt-fraction reaches the threshold limit for extraction (M2) as mentioned above. The extracted melt volumes for each of the TTG-types are computed at each time step and added to the volumes of the previous time-steps. It must be noted that we did not consider water and basalt as separate entities and instead, assumed the basaltic crust to be sufficiently hydrated to produce the different TTG-types (e.g., Sizova et al., 2015; Chowdhury et al., 2020).

FIGURE 2. Pressure and temperature conditions of Low-Pressure (LP), Medium-Pressure (MP) and High-Pressure (HP) tonalite-trondhjemite-granodiorite formation (Reproduced from Moyen (2011) and Jain et al. (2019) under the Creative Commons CCBY license). The solidii and liquidii for mantle peridotite and hydrous basaltic crust (Sizova et al., 2015) are also shown. See Appendix for details.

Comparison with the volume of ∼3 Ga old continental crust

Comparing the total continental crustal volume between models and nature is not straightforward, particularly because of the large variations among different growth models and the difficulties in extrapolating the numerical model results at planetary scales. Estimates for the volume of continental crust at ∼3 Ga vary from ca. 10 to 120% of its present-day volume between the different growth models (Figure 3) (cf. Cawood and Hawkesworth, 2019). This is because these models are based on different assumptions and datasets. While some growth models are based on the present-day distribution of the continental crust (e.g., Condie and Aster, 2010), some others are schematic (e.g., Armstrong, 1981). A third group of growth models also exist that are based on datasets such as the proportions of juvenile vs. reworked crust inferred from the zircon data, the evolution of atmospheric argon, and the trace element and isotopic composition of mafic-ultramafic rocks (Campbell, 2003; Belousova et al., 2010; Dhuime et al., 2012; Pujol et al., 2013; Roberts and Spencer, 2015; McCoy-West et al., 2019). Here, we consider the volume estimates of this group of growth models since: 1) they do not rely on the present-day crustal distribution; thus, accounting for the current paucity of the early rock record and 2) they all independently show that at least 65–70% of the present-day continental volume existed by ∼3 Ga. It must be noted that besides granitoids, the Archean continental crust is also comprised of greenstone belt successions made up of mafic-ultramafic volcanic and volcaniclastic rock units (Moyen and Laurent, 2018; Smithies et al., 2018). The estimates of continental growth models are however, meant for the felsic component mainly (Dhuime et al., 2012; Pujol et al., 2013; cf.; Cawood et al., 2013).

FIGURE 3. The growth models of the continental crust relative to its present-day volume (Reproduced from Cawood and Hawkesworth (2019) under the Creative Commons CCBY license). Note that the curves of Belousova et al. (2010), Dhuime et al. (2012), Roberts and Spencer (2015), Campbell (2003), Pujol et al. (2013), all of which used different proxies, predict similar volumes of continental crust at ∼3 Ga.

To measure crustal growth within our models, we introduce a parameter called the Differentiation Index (DI). Defined as the volume ratio of the continental crust to the upper mantle, this parameter quantifies how much continental crust has been extracted from the mantle and did not mix back (at least chemically); thus, accounting for the crustal recycling implicitly. Our approach follows the geochemical hypothesis that effectively continental crust has been extracted from the mantle (e.g., Rudnick and Gao, 2003). The DI value is estimated to be ∼2.4% for the modern Earth. It is calculated using the continental volume of ∼7.2 × 10⁹ km (cf. Cawood et al., 2013) and an upper mantle volume of ∼3.03 × 10¹¹ km (based on a spherical Earth of radius 6,371 km and an upper-lower mantle boundary at the depth of 660 km). The DI value at ∼3 Ga is ∼1.7%, which is calculated assuming the continental crustal volume to be ∼70% of the present-day volume and the same upper mantle volume. During the model evolution, we compute the total volume of TTG crust formed at each time-step (as mentioned above). We then use this accumulated TTG volume to compute DI at each time-step and trace over what duration its inferred value of ∼1.7% is achieved.

Model Results

Figure 4 illustrates the evolutionary stages of our reference numerical model, which has a surface yield stress of 50 MPa, a mantle T_p of 1,600°C (250°C higher than present-day), and a Moho-T of 600°C. The model evolution shows that two tectonic processes featuring local compression, called crustal thickening and dripping, occur under the influence of the convective flow of mantle. Thickening of the basaltic crust (Figure 4A) leads to phase transitions (i.e., eclogitization) within its deeper levels, which in turn increases the crustal density at those depths. Owing to a hotter geotherm set up by the hotter mantle and low yield stress, the crust could not support its gravitationally unstable eclogitized base and leads to the development of a Rayleigh-Taylor type instability within this region (Figure 4B). As this instability grows with time, a portion of the eclogitized crust drips back into the mantle (Figures 4C,D ). Such a process has been suggested to occur during the early Archean by previous studies (Johnson et al., 2014; Sizova et al., 2015; Rozel et al., 2017). Our model however, shows that: 1) individual drips may grow in length, reaching up to few hundreds of km (Figures 4A–D) and 2) individual dripping phases, i.e., from the appearance of the Rayleigh-Taylor-instability to the detachment stage, may last for ∼10 Myr. The cycle of crustal thickening evolving into dripping continue to occur over a periodicity of ∼20–30 Myr.

FIGURE 4. Evolution of the reference model with a surface yield stress of 50 MPa, mantle T_p of 1,600°C and Moho-T of 600°C. Mantle convection drives thickening of the basaltic crust (A), which eventually transform into a Rayleigh-Taylor type instability (B) and drips into the mantle (C–E). This cycle of crustal thickening followed by dripping continues. The residue forms during the production of tonalite-trondhjemite-granodiorite melt is shown in pink color. Within a period of ∼60 Myr, the model shows two major dripping phases: the first phase (B–E) between ∼9 and 34 Myr and the second phase (F–J) between ∼49 and 57 Myr. Note that this second phase is transferred across the (periodic) boundary from left to right in the model (G–H) around the ∼51–52 Myr time-period. The TTG formation mostly occurs during the dripping phase as compared to the thickening phase.

Regarding the formation felsic crust, the model shows that drips are more efficient in producing TTGs as compared to the basal part of a thickened crustal segment (cf. Figures 4A–D). The narrow tail region of a basaltic drip easily reaches supra-solidus temperatures due to the conductive heating from the surrounding hot mantle, resulting in its melting and the formation of TTGs. It must be noted that in our models we have shown the residues (pink colored lithology) formed after the production of TTGs. Therefore, although residues get transported at great mantle depths with the drips (Figures 4C,D), they have necessarily formed at depths <75 km since we did not allow melting of the basaltic crust beyond this depth limit. Figure 5 shows how the volumes of modeled LP-TTG and MP-TTG melts as well as their combined volume increase with time during the model evolution. The relative proportions of LP- to MP- TTG melts remain fairly constant (varying from ∼22%:78% to ∼30%:70%) around the ∼1:3 ratio throughout the model evolution (Figure 5B), which agrees well with the rock record (Moyen, 2011; Moyen and Laurent, 2018). In terms of the absolute volume, prominent spikes implying significant growth in the TTG volume, can be seen at ∼17.5, 56.8, and 100 Myr in the melt volume vs. time plot (Figure 5A). Notably, all these time-instants correspond to the active phases of dripping (Figure 5C). In contrast, melt volumes do not show any major increase when crustal thickening dominates the model evolution (Figure 5A), implying only limited mafic crustal melting in this setting.

FIGURE 5. Evolution of (A) the volume of low-P tonalite-trondhjemite-granodiorite (TTG), medium-P TTG and total TTG crust [represented via Differentiation Index (DI)], and (B) the relative proportion of low-Pressure (LP) to medium-Pressure (MP) TTGs through model-time. The dashed green line in (A) represents the geologically constrained DI value (≥1.7) at 3 Ga (termed as Archean DI), corresponding to a continental volume of 70% of its present-day present-value (see text for details). Note that Stage-I, II, and III represent the major changes in the volume of TTG melts and total TTG crust. These stages correspond to (C) the three major, successive dripping phases in the model. The pink-coloured material represents the residue of the TTG melts. It must be noted that the residues form at shallow depths (P < 2 GPa) where LP- and MP- TTGs are prescribed to form, and then they get transported to deep mantle owing to greater density. The relative proportions of LP-TTGs to MP-TTGs remain around the ∼1:3 ratio throughout the model (except during the initial, few million years), which agrees with the preserved igneous rock record.

To test the effects of higher yield stress and Moho-T, we have simulated a model with a surface yield stress of 200 MPa and a Moho-T of 700°C under the mantle T_p of 1,600°C. This model also feature dripping (Figures 6A–C) and crustal thickening phases (Figures 6C–F) similar to the reference model. But, in this model, crustal thickening phases exhibit melting of the basaltic crust and TTG production in noticeable volumes (Figures 6E,F). A higher Moho-T along with a higher yield stress value leading to sustaining of thickened crustal domains over ∼20–30 Myr, facilitates melting at the base of the thickened basaltic crust. This substantiates our contention that Moho-T and yield stress values together influence the production of TTGs within thickened crustal domains. However, even within this setup, the overall volume of the TTGs is controlled by the melt produced during the dripping phase (Figures 6A–G). However, the relative proportion of modeled LP- to MP-TTG melt do not corroborate to the desired value of ∼1:3, rather they show sub-equal abundances (Figure 6H). Models with other values of yield stress, Moho-T and mantle T_p do not show any major deviation in their evolution from that of the reference model.

FIGURE 6. Evolution of a model with yield stress of 200 MPa, mantle T_p of 1,600°C and Moho-T of 700°C. While the overall model evolution (A–F) resembles that of the reference model, an enhanced production of tonalite-trondhjemite-granodiorites (TTGs) at the sites of thickened crust (E,F) is observed in this model. Evolution of (G) the TTG melt- and total TTG crustal volumes, and (H) relative proportions of (low-pressure) LP- to (medium-pressure) MP-TTGs through time. Crustal growth during the thickening phases (D–F) can be seen but the major growth still occurs during the dripping phases (A–C). Additionally, this model predicts an equal proportion of LP- and MP-TTGs during the Archean, which contrasts with the preserved rock record.

Furthermore, the model results show that TTG production and the resulting DI values increase with the lowering of lithospheric yield strength for any given mantle and Moho temperature conditions (Figure 7). When the mantle T_p range is within 1,500–1,600°C, only the models having a lithospheric yield stress of ∼50 to <100 MPa produces enough TTG volume such that their DI values approach or surpass the geological constrained DI value (∼1.7%) of the Earth approximately at 3 Ga within few hundreds of Myr (Figures 7A,B). However, at low yield stresses, the DI reaches ∼1.5–1.7% within ≤100 Myr when the Moho-T is 700°C irrespective of the mantle T_p (Figures 7A–C), implying an extremely rapid growth of felsic crust. In contrast, the DI values reach ∼0.5–1.2 in 100 Myr when the Moho-T is ∼600°C (Figures 7A–C); thereby, suggesting that the felsic crust grows steadily under these conditions and the DI values likely reach the ∼1.7% value in a few hundreds of Myr. A colder Moho (∼500°C) or mantle (<1,500°C T_p) seems to limit the TTG production such that the DI values remain <0.5 in 100 Myr (Figures 7A–C); thus, favoring a slow crustal growth. However, the estimated DI values likely reach the ∼1.7% in less than 1.5 Gyr for colder Moho (∼500°C) and/moderate T_p (∼1,540°C) and/higher yield stress (250 MPa) suggesting that slow crustal growth rates require excessively high lithospheric yield stress values. Even when the relative proportion of LP- and MP-TTGs are considered (Figures 7D–F), models with a combination of hotter mantle (>1,550°C T_p), moderately high Moho-T (∼500–600°C) and low lithospheric yield stresses (<150 MPa) provide the best agreement to the natural rock-record. A significantly hotter Moho (≥700°C; Figures 7D–F) or colder mantle (<1,500°C T_p; Figure 7F) would strongly favor or disfavour the production of LP-TTGs, respectively, (even for weak lithosphere), which in turn deviates the ratio of LP- to MP-TTGs from their geologically constrained value (∼1:3).

FIGURE 7. Plots showing (A–C) the total tonalite-trondhjemite-granodiorite (TTG) crustal volume cast in terms of the Differentiation Index (DI) value and (D–F) the proportions of low-pressure (LP-) to medium-pressure (MP-)TTGs as a function of lithospheric yield stress, as obtained from our numerical models with different mantle T_p and Moho-T values. The DI values are computed at the end of the first 100 Myr of the model evolution in each case. The dashed blue lines in (A–C) represent the geologically constrained DI value (≥1.7) at 3 Ga, while in (D–F) they represent the relative proportion of the LP-to MP-TTGs as preserved in the rock record (see text for details). Note that the models with low lithospheric yield stresses (at hotter mantle and moderate hotter Moho conditions) produce the DI values that are likely to result a continental volume of 70% of its present-day value (DI = ∼1.7) in few hundreds of Myr while matching the rock-record of the proportions of LP- and MP-TTGs.

Discussion

Yield Stress of the Archean Lithosphere

Our modeling shows that Archean lithosphere has somewhat lower yield strength compared to the present-day lithosphere (>100 MPa). The presence of such weaker lithosphere (together with hotter mantle) is critical for producing the large volumes of the Archean continental crust. This is because lowering the lithospheric yield stress facilitate the dripping-off of basaltic crustal segments, where the most favourable melting conditions (∼0.7–1.8 GPa; > 800°C) ascertained for the production of TTGs (e.g., Palin et al., 2016; Johnson et al., 2017; cf.; Moyen and Laurent, 2018), are easily attained. Furthermore, weaker lithosphere (with moderately high Moho-T) can aptly replicate the relative dominance of MP-TTGs over LP-TTGs similar to what is preserved in the rock-record. Higher values of yield stress (>100 MPa) would impede the formation of drips and therefore, would result in a lower abundance of MP-TTGs, particularly if the crust is not significantly thick and/or thermally pre-conditioned for melting. Consequently, we suggest that the yield strength of Archean lithosphere likely varied between ∼50 and 100 MPa, provided the mantle was significantly hotter than today.

The yield stress of material is temperature dependent in most engineering applications (Liu et al., 2002; Senkova et al., 2006; Heidarzadeh et al., 2014; Wu et al., 2014; Altenberger et al., 2015), as well as geological laboratory experiments related to different rock compositions (Evans and Goetze., 1979; Long et al., 2011; Sygała et al., 2013; Idrissi et al., 2016; Hansen et al., 2019). Therefore, owing to the hotter thermal structure and different crustal composition of Archean lithosphere relative to the present-day crust (Chardon et al., 2009; Dhuime et al., 2015; Tang et al., 2016; Brown et al., 2020), it is expected to have a lower yield strength. In fact, studies of how a lithosphere deforms have shown that it involves three main mechanisms: diffusion creep, dislocation creep and Peierls mechanism (Regenauer-Lieb et al., 2006), and ultimately results in a shear heating controlled rupture of the ductile harder lithospheric layer (Regenauer-Lieb et al., 2001; Thielmann and Kaus, 2012; Stern and Gerya, 2018). The yield stress that we apply is an average of a large number of processes, which are primarily thermally controlled and secondary pressure controlled, and determine the rupture of a lithosphere whose ultimate strength is thermally controlled (Braeck and Podladchikov, 2007). Consequently, the effect on the thermally driven tectonic style during the Archean, when the uppermost mantle boundary layer was hotter, is expected.

The results agree with laboratory experiments of rock deformation showing a lower yield stress with increasing temperature (Evans and Goetze., 1979; Long et al., 2011; Sygała et al., 2013; Idrissi et al., 2016; Hansen et al., 2019). Postulations of weaker lithosphere due to hotter mantle temperatures, higher radioactive heat generation and melt impregnation have been suggested by previous studies as well (e.g., van Hunen and van den Berg, 2008; Sizova et al., 2010; Magni et al., 2014; Rozel et al., 2017). Our results converge with these studies in finding a need for a weaker lithospheric strength. Overall, the lithospheric weakening could have occurred due to both viscous strength variations as well as yield strength variations with higher temperature conditions in the Archean Earth.

Tectonic Settings for the Early TTG Formation

Our modeling also bears implications for tectonic settings where the Archean felsic crust formation may have occurred. Previous studies have invoked both plate tectonic settings (mainly subduction) analogous to the present-day scenario (e.g., Foley et al., 2003; Harrison, 2009; Korenaga, 2018; Turner et al., 2020) or at various tectonic sites hosted within a stagnant-lid type tectonic regime (e.g., Smithies et al., 2009; Sizova et al., 2015; Johnson et al., 2017; Nebel et al., 2018; Smithies et al., 2019) for the formation of TTG-dominated crust during the Archean. Although the debate continues, recent studies are increasingly favoring the stagnant (squishy)-lid to transitional mobile-lid tectonic settings (not like modern plate tectonics) for the Archean crust formation (Johnson et al., 2014; Sizova et al., 2015; Chowdhury et al., 2017; Condie, 2017; Hawkesworth et al., 2017; Rozel et al., 2017; Johnson et al., 2017; Cawood et al., 2018; Moyen and Laurent, 2018; Capitanio et al., 2019; Smithies et al., 2019; Chowdhury et al., 2020). Tectonic processes like crustal thickening, dripping and local crustal overturns are suggested as settings for TTG production, particularly during the pre-3 Ga time-period (e.g., Smithies, 2000; Johnson et al., 2014; Sizova et al., 2015; Johnson et al., 2017). Our study highlights the relative dominance of these different non-plate tectonic mechanisms in producing TTGs. The reference model shows that production of TTG dominantly occurs during the dripping phase instead of during the crustal thickening phase. This partly contradicts the results of Sizova et al. (2015), which showed that crustal thickening promotes TTG production. However, they implemented an unusually high, initial Moho-T of ∼1,200°C, which is markedly higher than 600°C value used in our model. Such lower Moho-T values together with the rapid transition of the thickened crustal domains into drips (within <10 Myr) at low yield stresses, impeded large-scale melting of the mafic crust during crustal thickening in the reference model. This becomes evident when we increase the Moho-T and yield stress compared to the reference model and observed a greater proportion of melting within the thickened crustal domains.

Nonetheless, in all our models, the dripping process dominates TTG formation. We infer that this is because dripping continuously brings newer basaltic material to the melting region, maintaining a steady melting of the source rock (Figures 4 and 6). Secondly, dripping leads to crustal thinning and the resultant, mantle upwelling around the locale of the drip. This heats up the surrounding crust, which facilitates its subsequent melting. In contrast, only a thin layer of the lower most crust reaches melting conditions within a thickened crustal segment (Figures 6E,F). Given that crustal thickening occurs at a local scale and is not laterally pervasive, the resultant melting and TTG production is also limited. We, however, do not suggest that crustal thickening did not contribute toward Archean crustal growth. Instead, we imply that this may have contributed significantly only when an unusually hot (Moho-T > 1,000–1,200°C) and thick (>35–40 km) mafic crust was available for melting conditions that may have existed in the Archean but perhaps not everywhere (cf. Cawood and Hawkesworth, 2019).

Conclusions

The preserved record of early Archean felsic crust is more consistent with the presence of weaker lithosphere during active tectonism at those times. Tectonic activity involving such weaker lithosphere atop hotter mantle would have led to the production of voluminous TTG melt prior to ∼3 Ga, consistent with the inference of numerous crustal growth models. Our study further shows that TTG production may have dominantly occurred within mafic crustal drips, whereas crustal thickening made an only subordinate contribution unless the crust was significantly thicker. The low yield stress facilitated the drips to grow in length up to few hundreds of km and keeping them stable over tens of Myr. Furthermore, weaker lithosphere together with the dominance of dripping can also account for the large abundance of MP-TTGs relative to the LP-variant within the early Archean rock record.

Data Availability Statement

The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Author Contributions

PG—Overall research project is carried out from identifying the research problem, model development, and evaluation of results to drafting the manuscripts. GM—Provided support for the numerical model development, applying Geophysical significance, and Initial manuscript corrections. PC—Provided support in interpreting the model results, comparing them with the rock-record and manuscript development. PAC—The whole research project was guided by PAC from identifying the research problem to final manuscript development.

Funding

MIPRS and MGS PhD scholarship provided by Monash University, and Australian Research Council grant FL160100168.

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.

Acknowledgments

We acknowledge the computational resources provided by the Monash eResearch Center and eSolutions-Research Support Services (for the MonARCH HPC Cluster) and National Computational Infrastructure (NCI), Australia. We thank Jeroen van Hunen (Durham University) and two reviewers for their constructive comments that have immensely improved the manuscript. We also thank Christopher M. Gonzalez, Shayne McGregor, Fabio Capitanio, and John Mansour (Monash University) for insightful discussions.

Appendix

Appendix A1: Pressure-Temperature Conditions for Low-P, Medium-P and High-P TTG Formation

The pressure-temperature (P-T) conditions of LP- and MP-TTGs formation during the Archean (Figure 2) can be written as the following functions (Jain et al., 2019):

$760-60 {\left(P-1\right)}^2<T<\mathrm{1,000}-150{\left(\frac{P-1.2}{1.2}\right)}^2(\mathrm{A1})$

$-0.5 \left(\frac{T-870}{220}\right)<P<1.5+0.7\left(\frac{T-700}{200}\right)(\mathrm{A2})$

where, T is in °C and P is in GPa. Low-P TTGs generally forms at P < 0.8–1 GPa while the medium-P TTGs forms at P > 0.7–1 GPa, depending upon the actual source rock composition (Moyen, 2011; Palin et al., 2016; Johnson et al., 2017). For simplicity, our employed functions put this pressure bracket between low-P and medium-P TTGs at 1 GPa, following the traditional definition scheme put forward by Moyen (2011). The P-T conditions for forming the High-P TTG via partial melting of basalts is given by (Jain et al., 2019):

$\mathrm{1,000}<T<\mathrm{1,100}+50{\left(\frac{P-3.5}{3.5}\right)}^2(\mathrm{A3})$

$2.35+0.15\left(\frac{T-\mathrm{1,000}}{100}\right)<P<2.5(\mathrm{A4})$

FIGURE A1

FIGURE A1. The average basaltic crustal thickness relevant for different Moho-T values for a model with yield stress of 50 MPa and mantle T_p of 1,600°C.

TABLE A1

TABLE A1. Melt extraction threshold values.

Appendix A2: Solidii and Liquidii for Different Lithologies Used in the Modeling

For mantle melting, we used the P-T dependent peridotite melting model parametrized in Sizova et al. (2015).

$T_{\text{sol}}=1120.85+133P-5.1P^2(\mathrm{A5})$

$T_{\text{liq}}=1799.85+114P(\mathrm{A6})$

where P is in GPa and T is in ^°C. For the melting of hydrated basalt, we used the following parametrized solidus and liquidus functions (Sizova et al., 2015):For, P < 1.6 GPa

$T_{\text{sol}}=699.85-\frac{\mathrm{70,400}}{\mathrm{1,000}P+354}+\frac{\mathrm{77,800,000}}{{\left(\mathrm{1,000}P+354\right)}^2}(\mathrm{A7})$

$T_{\text{liq}}=1149.85+105P(\mathrm{A8})$

for, 1.6 GPa < P < 10 GPa

$T_{\text{sol}}=661.85+3.5P+6.2P^2(\mathrm{A9})$

$T_{\text{liq}}=1149.85+105P(\mathrm{A10})$

References