{'fraction_non_alphanumeric': 0.04221317654153475, 'fraction_numerical': 0.028297904417307403, 'mean_word_length': 4.916599768084916, 'pattern_counts': {'":': 0, '<': 27, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Craton Formation in Early Earth Mantle Convection Regimes

How the geological record of cratons reconciles with the tectonic environments in which they formed has remained debated. We use 2D Cartesian geometry numerical models of mantle convection varying temperatures from present day to Archaean‐inferred values, to address the formation of cratons, accounting for melt depletion‐dependent rheological stiffening. For mantle temperatures comparable to present day, melting is negligible and the convective regime depends on the strength of the thermal lithosphere. For mantle potential temperatures higher than present day, high depletion degree and large depleted mantle volumes are formed at low lithospheric strength and high surface mobility, whereas these are negligible beneath a poorly mobile lithosphere. When compared to the models, the record of tectonics and large volumes of high‐degree depleted mantle in Archaean cratons is best explained by a lithosphere initially prone to yielding and mobility. At high mobility, large depletion favors the progressive differentiation of the thermochemical lithosphere, which stiffens and thickens with increasing mantle temperatures. The ensuing reduced heat flow atop a hotter mantle is in agreement with the inferred Archaean thermal evolution, and may rule out the viability of a stagnant lid for the early Earth. Large‐scale depletion stiffening resists plate margin formation and this wanes as heat production decreases, thus may hold the key for the establishment of plate tectonics during secular cooling.

A viable tectonic model for the early Earth must additionally satisfy the constraints posed by the thermal evolution in the Archaean. The role of the lithosphere in mantle convection is critical to the regulation of thermal loss of the Earth (e.g., Korenaga, 2013;Moresi & Solomatov, 1998). Extrapolating convective regimes to Archean hotter mantle conditions leads to excessive mantle temperatures that are not observed, called the "Archaean paradox" (Christensen, 1985;Lenardic, 2006). Dehydration (Korenaga, 2006) and cooling (Moore & Webb, 2013) of the mantle during melt extraction have been proposed to play a more prominent role in a hotter mantle (Moore & Lenardic, 2015), thickening the lithosphere, critically reducing heat flow and explaining petrological evidence of mantle warming in the Archaean, preserved in cratons (Herzberg et al., 2010). Therefore, while the link between the geological record and the geodynamic regime is acknowledged in the Phanerozoic, understanding what the regime of the early Earth and its thermal evolution were, has remained difficult to unravel.
In this work, we address the link between early Earth tectonic environments and mantle convection regime, and provide a quantification of mantle melt-depleted occurrence to compare with that of cratons. We aim at reproducing the conditions for different mantle convection regimes, from stagnant to ML (e.g., O'Neill et al., 2007;Stein et al., 2013), and test the impact of melt extraction and stiffening of the lithosphere, previously proposed on the basis of scaling analysis (Korenaga, 2003). Previous work highlighted the role of melt-related heat extraction and weakening in crustal differentiation (Chowdhury et al., 2017;Fischer & Gerya, 2016a;Gerya et al., 2014;Piccolo et al., 2019) and cooling histories (Lourenço et al., 2018;Moore & Webb, 2013;Nakagawa & Tackley, 2012;Rozel et al., 2017). Other studies have followed the opposite approach and focused on the residual mantle cooling and stiffening as melt is extracted to illustrate the differentiation of the thermochemical lithosphere, the formation of cratons and the implications for the Earth's thermal evolution (Capitanio et al., 2019(Capitanio et al., , 2020Korenaga, 2006;Liu et al., 2021;Wang et al., 2014Wang et al., , 2018. These different approaches are, in fact, complementary as they are all processes relevant to the formation and differentiation of the crust and mantle (Arndt et al., 2002). Here, we follow the latter approach and address the impact of progressive viscosity increase as melt is extracted on mantle stiffening and convection. The outcomes are temperatures, melt-depletion degrees, and compositions of the lithosphere that are compared with depletion degree, volumes, and structures of cratonic keels (Griffin et al., 1999;Pearson & Wittig, 2014).
We show that depletion degree along with the volumes and distributions of depleted lithospheric mantle preserved in cratons pose unique constraints to the interpretation of the tectonic environments in which they formed. High depletion degrees with shallow emplacement depths and large depleted mantle volumes preserved in Archaean keels are matched only by models with relatively low yield strength, suggesting that cratons formed in an environment favorable to high mobility, not a stagnant lid. We then speculate on the role of the craton-forming processes in the thermal evolution of the Archaean, by comparisons with scaling analysis (e.g., Foley, 2018;Solomatov & Moresi, 2000;Stein & Hansen, 2008;Stein et al., 2013;Weller & Kiefer, 2020), and show how the inferences based on the sparse geological record may be critical to the interpretation of the thermal evolution of the early Earth.

Modeling Approach
We focus on the formation of stable lithosphere as the result of mantle convection through Earth's early history, when mantle temperatures were higher and melting and depletion were more pervasive. We account for the effect of melt extraction implementing a dehydration-and latent heat-dependent rheology (Ito et al., 1999;Korenaga, 2003;Phipps-Morgan, 1997;Wang et al., 2018) during melting of the mantle. While melting is implemented as in other studies (e.g., Fischer & Gerya, 2016a;Johnson et al., 2014;Sizova et al., 2015), we model the stiffening of the residual mantle, as it undergoes melting and melt extraction. In a hotter mantle, melting occurs in a deeper layer, thereby differentiating the lithosphere into a thick thermochemical boundary layer where viscosity increases due to a refractory olivine residual (Korenaga, 2003) and to cooling due to latent heat extraction (Moore & Webb, 2013). Because the melting depth increases with mantle temperature, the thermochemical boundary layer must have been thicker in the Archaean, potentially impacting the convective regime of the early, hotter Earth (Korenaga, 2006;Moore & Lenardic, 2015), which is addressed here.
A key parameter in our modeling is the yield strength of the lithosphere. This parameter is known to control the mobility of the lid and its ability to subduct in the mantle and enhances heat transport, thereby impacting the surface tectonics and mantle melting, as well as regulating the heat budget of convection (e.g., Moresi & Solomatov, 1998).
Melting temperatures of basalts preserved in cratonic crust suggest mantle potential temperatures up to +260°C with respect to present day, in the Archaean (Herzberg et al., 2010), and temperature change rates varying from +50 to +100 K Gyr −1 , in the Archaean, to -50--100 K Gyr −1 in the Phanerozoic (Jaupart et al., 2015). These temperatures and trends are achieved using mixed heating conditions, that is, basal heating and internal radiogenic heat production, in numerical convection models (e.g., Lowman et al., 2001;O'Neill et al., 2007;Stein et al., 2013) in a two-dimensional Cartesian geometry. This approach focuses on the lithospheric processes (Chowdhury et al., 2017;Fischer & Gerya, 2016a;Gerya et al., 2014;Piccolo et al., 2019) reproducing potential temperatures and their change to provide quantities comparable to the petrological record.
No attempt is made here to reproduce the exact thermal budget of the Earth. The heating mode of convection depends strongly on the choice of boundary conditions, model space geometry and dimensions (e.g., Guerrero et al., 2018;O'Farrell & Lowman, 2010;O'Farrell et al., 2013), and mantle temperature may vary largely with model aspect ratio (Lowman et al., 2001). To overcome this limitation, we tested the effect of different thermal boundary conditions and internal heat production, in upper and whole mantle convection models, and assessed the impact on the temperature and melting history of the models.

Governing Equations and Numerical Method
Mantle convection is modeled as the flow of a viscous fluid at very high Prandtl number in a two-dimensional Cartesian geometry. We use the geodynamic framework Underworld  to solve the equations of conservation of mass, momentum, and energy using an Eulerian Finite Element Method with Lagrangian particles embedded in the elements (FEM-PIC). The Lagrangian particles allow for multimaterial properties, tracked throughout the history of the model. This is key to the implementation of history-dependent melting and melt extraction used for the melt depletion-dependent rheology.
The conservation of mass equation, enforcing an incompressibility constraint, is where is the velocity vector. The conservation of momentum equation is with the stress tensor and = the force term, the density, and the gravity. The stress tensor is defined as where is the deviatoric stress, p the pressure, and is the identity tensor.
Using the Boussinesq approximation, the conservation of energy equation reduces to with T′ the nonadiabatic temperature, t the time, the thermal diffusivity, and H r the radiogenic heat generation rate per unit mass (Turcotte & Schubert, 1982): where H is the internal radiogenic heat production and c p is the heat capacity. Latent heat is not considered in the energy equation, as this implies assumptions on melt transport and emplacement that might bias the outcomes (Moore & Lenardic, 2015). In order to respect the energy conservation constraint, we consider the system closed, whereas we include the effect of latent heat in the rheology.
The Boussinesq approximation neglects the effect of density change except its buoyancy effect (Christensen, 1995), therefore the adiabatic temperature is defined as a function of depth y: with α as the thermal expansivity. The equivalent temperature defined at the surface is the potential temperature T P . We do not model secular cooling and H r remains constant to keep the models in the same regime throughout their evolution.
We define two nondimensional numbers to characterize the evolution of the models, the internal Rayleigh number Ra i (Solomatov, 1995): where 0 is the reference density, T i is the internal temperature, the horizontal average of the models' output, D the model thickness. And ( i ) is the viscosity at the internal temperature, and the Nusselt number: where kT/D is the flow in the hypothetical case of heat released by conduction only through the mantle thickness, and is the surface average heat flow.
The equation of state accounting as function of temperature difference ΔT and the depletion fraction F (see below) is where ∆T is the superadiabatic temperature difference and the change in density due to depletion. This depletion-dependent density varies linearly by 0.726 kg m −3 per depletion percent (Schutt & Lesher, 2006), whence = −72.6 kg m −3 . This value follows previous geodynamic studies (e.g., Ito et al., 1999;van Hunen and Moyen, 2012).

Constitutive Laws and Dehydration Stiffening
Key to the constitutive laws used is the melt-depletion degree F. We evaluate the melt fraction F (T, p) as a function of the super-solidus temperature T ss and pressure, following McKenzie and Bickle (1988), and using the parameterized solidus and liquidus, T sol and T liq , respectively, from Katz et al. (2003). The dimensionless super-solidus temperature T ss is defined as with a = 0.5, b = 0.25, c = 0.4256, and d = 2.988 (McKenzie & Bickle, 1988).
The melt is considered extracted, that is, it is not added to the overlying crust, whereas we consider the two major effects of melt extraction in the rheology: Dehydration and latent heat. These two effects occur at the same time and result in viscosity increase of the residue by cooling and dehydration, due to heat advection and water extraction from melting olivine, respectively. These are implemented in a constitutive law, which is a function of the depletion degree F and the latent heat T* (see below). We do not physically extract the melt from the models and instead change their density. This approach, commonly used in modeling at this scale (see Gerya, 2014), aims at minimizing the violation to the mass conservation constraint in the overall modeling space, although preliminary test shows that the outcomes are not different from models with no depletion-dependent density. Because we do not advect heat with the melting from the model, as prescribed by the energy conservation constraint, but only use the effect of cooling on the rheology, the estimates of volumes and melting-degree represent an upper bound, which is estimated in Figure 1. We show the case of a mid-ocean ridge geotherm in a hotter mantle, T P = 1,560°C, where melting is maximized ( Figure 1a). In Figure 1c, melting curves are shown with latent heat, F*, and without (thicker lines). An integration of the volumes in the absence of latent heat yields volumes that are 1.58-1.7 times larger, respectively, for the wet and dry mantle. The choice of dry melting provides minimum overestimates. Due to the large uncertainties in temperature and water content in the early Earth mantle and the experimental nature of available data, the comparison with the observables is best achieved by the depletion degrees, which are well constrained (Lee et al., 2011), allowing avoidance of any influence of the hydration assumption. Further discussion is proposed in Section 2.2 and 4 ( Figure 2).
Relevant to our modeling is the effect of water partitioning into the melt on the rheology of rocks and latent heat, which depend on F. We focus on the partitioning of water during fractional melting. This is an incremental process, then the concentration of water in the solid phase, C S, varies as (Karato, 2008): where C 0 is the initial concentration of water and partition coefficient of water in the melt phase is k = C S /C M , where C M is the concentration in the matrix. Fractional melting, assumed here, is better suited for settings where melt is extracted (Braun et al., 2000;Karato, 1986). Other studies adopted a batch (linear) partitioning law, assuming the difference with fractional melting is negligible (Braun et al., 2000;Phipps-Morgan, 1997). The latent heat is defined as T* = T -F 260°C, following Phipps-Morgan (1997). Mantle adiabats for present-day temperatures beneath a mid-oceanic ridge, depletion, and viscosity for different water contents versus depth. (a) Adiabatic temperature T and T*, with cooling due to latent heat, for T P = 1,350°C, with solidi for dry and wet (water content X H2O in weight fraction) and water saturated mantle. (b) Depletion degrees F and F* for dry (solid) and wet (dashed) as a function of T and T*, respectively, (c) viscosities of dry (solid black line) and wet (black dashed line), compared to Hirth and Kohlstedt (1996)  Finally, we embed the effect of dehydration and latent heat (e.g., Korenaga, 2006) during melting using a temperature-dependent creep law (e.g., Rolf et al., 2018;Rozel et al., 2017). The water content-dependence of strain rate in olivine depends quasilinearly on water content, that is, r ∝ , where r = ∼1 (Hirth & Kohlstedt, 1996). Then, the temperature-and dehydration-dependent creep law for viscosity, that is, ∝̇r , with r = 1, becomes: where A = A 0 C 0 is a prefactor, 0 is the reference viscosity, E a is the activation energy and R the gas constant and T* is defined above (values in Table 1). Noting that (1-F) (k -1)/k ≈ exp (F(1k)/k), for F ≤ 0.5, we chose realistic parameters in Equation 13 that yield negligible difference between the results presented here and in Capitanio et al. (2020). Note that 13 does not include the adiabatic temperature-and pressure-dependence, as customarily done in large-scale models (Rolf et al., 2018;Rozel et al., 2017). For the low activation volume of dry olivine these two effects oppose each other having minor impact on the depth dependence of the viscosity profile in the upper mantle (King, 2016).
The implementation of rheological stiffening during melting follows several previous studies (Braun et al., 2000;Choblet & Parmentier, 2000;Hirth & Kohlstedt, 1996;Ito et al., 1999;Korenaga, 2003;Phipps-Morgan, 1997;Wang et al., 2018). As melt is produced in a two-phase system, the strength of rocks weakens nonlinearly, with the largest bulk viscosity drop reached for melting fractions F ≤ 0.1 (e.g., Rosenberg & Handy, 2005). In contrast, in a single-phase system, that is, in systems where the melt is removed (Karato, 1986), dehydration induces an increase in viscosity of the residual mantle. However, when the effects of melt weakening and dehydration stiffening are combined, the latter overcomes the effects of melt weakening and an increase of effective viscosity is the result (Braun et al., 2000;Choblet & Parmentier, 2000).
The stiffening resulting from dehydration and cooling can be broken down to three domains (Hirth & Kohlstedt, 1996;Mei & Kohlstedt, 2000): (a) During wet melting, viscosity increase is negligible (∼10-fold), (b) across the wet-to-dry transition, the viscosity has a sharp jump, of a factor ∼100, (c) during dry melting, viscosity further increases ∼3-4 times per 10% dry depletion (Figure 1c, HK96, dotted grey line). The resulting viscosity is reproduced here using k = 2.1 × 10 −2 for F ≤ 0.1 (domain ii) and k = 1.57 × 10 −1 for F > 0.1 (domain iii; Figure 1c, solid line). The partitioning coefficients chosen are compatible with water partitioning k ∼ 10 −2 -10 −3 , ensuing rapid dehydration of the solid phase, for small F. Following published studies (Ito et al., 1999; Phipps-Morgan, 1997), we neglect the initial stiffening achieved in the wet melting (domain (a) (Figures 1b  Plasticity is implemented using different laws for the mantle and the crust. For the mantle, including dehydrated mantle, we implement a pseudo-plastic flow law using a Drucker-Prager-type yield criterion: with 0 the cohesion stress at surface conditions and the internal friction angle, using the standard value of 34°. To simulate a hydrated crust, we use a weaker rheology weaker than the lithosphere. The crust is defined by the inflow (toward lower temperatures) across the nonadiabatic basal (Moho) temperature of 320°C. There is no attempt to address the nature of the crust, whether basaltic or felsic. The rheology of the crust follows a Byerlee's law for near-surface conditions, where the cohesion vanishes, 0 = 0, where the superscript C indicates the crust, which is defined as The different depth-dependent coefficients in Equations 14 and 15 allow the crust to be weaker than the mantle at near-surface conditions.
The plastic flow law is implemented as following: where Y is either the mantle or the crust's value, that is, Equation 14 or 15, and II = √ ( ∶ ) ∕2 is the square root of the second invariant of the strain rate tensor. The latter is defined aṡ The maximum viscosity in the models is capped self-consistently by the plasticity, which is depth dependent and yields maximum values of ∼10 25 Pa s.

Model Setup and Parameters
We have run different suites of models, for a total of 106, with varying values of internal heat production H, bottom temperatures boundary conditions T b , depth of the mantle, lateral boundary velocity conditions, and cohesion stress 0 ( Table 2 and Supplementary Table ST1 in Supporting Information S1). Using mixed heating-mode convection, that is internal heating and bottom temperature conditions, as well as internal heating only shows that different heating modes may result in similar geotherms and thermal evolutions yielding comparable melting histories, geotherms, and lithosphere viscosity (Supplementary Figure S1 in Supporting Information S1), which is the focus of this work. A resolution test shows negligible variations in temperatures, velocities magnitudes, and melt volumes across the resolutions tested, while only the heat flow peak varies <∼20% when the resolution is doubled, although the overall evolutions show no dependence on resolution (Supplementary Figure S2 and S3 in Supporting Information S1). This shows that the FEM-PIC method yields accurate solutions for the stresses and melting in the lithosphere at the resolution of the models presented.
We present a selection of 29 models, which extend to the upper mantle only (e.g., Fischer & Gerya, 2016b;Liu et al., 2021;Perchuk et al., 2020;Piccolo et al., 2019;Wang et al., 2014). In this suite of numerical experiments, we model convection in the upper mantle in a space of 660 × 7,920 km, that is, aspect ratio 1 × 12, discretized in 64 × 768 elements, embedding a total of 983,040 Lagrangian particles. We present the models with wrap-around (periodic) boundary conditions on the side walls and no slip at the bottom, while the top is free slip. The list of common parameters used in the models is presented in Tables 1 and 2.
In the model suite presented, we use H = 0, 1, 2, and 3 × 10 −11 W kg −1 (Table 2) and bottom temperatures T b = 1,628°C, while we fix the top boundary temperature T 0 = 0°C. The range of internal heat production covers the range on Earth, which varied exponentially between present-day values of H = 3-7 × 10 −12 W kg −1 (Jaupart et al., 2015) to values in excess of 3 × 10 −11 W kg −1 , in the Hadean, for a bulk silicate Earth (BSE) (Turcotte & Schubert, 1982). The values of the cohesion are varied between 0 = 70 and 1 MPa (Table 2), covering the range from laboratory-constrained values for different rock compositions from dry to fluid filled (Gerya, 2009).
The initial configuration, common to all models, is achieved through 500 Myr of basally heated convection and no melting, at Ra i = ∼10 7 ( Figure 2). From this configuration, the models run with melting and basal and internal boundary conditions for a total time of 1 billion years (Gyr). The internal heat production and the basal temperature are kept constant in the models (Table 1), that is, no decaying of radiogenic heat is implemented (e.g., Weller & Kiefer, 2020). The convection regime, and whether this has high or low surface mobility, is not strongly affected by small, ≤∼3 times, variations of internal heat (Langemeyer et al., 2018). Our choice allows excluding  the effect of secular cooling due to radiogenic heat production decay and provides a simplified model that can be compared to scaled analysis of thermal evolutions.
Additional testing includes a suite of 32 models to address the role of free-slip boundary conditions, that is reflective boundaries (FS models, Supplementary Table ST1 in Supporting Information S1). These boundary conditions are commonly used to support scaling analysis, although they might not have a physical equivalent for the Earth and are not presented. Also, we ran 21 models to test the role of a lower mantle (LM models, Supplementary Table ST1 in Supporting Information S1). We use wrap-around side walls, while the bottom boundary is located at 2,970 km depth, and has a temperature T b = 2321°C, and a viscosity increase by a factor 30 in the lower mantle. The transition between upper and lower mantle is located at y = 660 km depth, whereas the height of the phase change is implemented following published work (Christensen & Yuen, 1985). The initial temperature distribution in these models is obtained as described above, yielding a Ra i = ∼10 9 . Finally we have ran 18 models varying the bottom temperature, T b = 1628, 1952, and 1302°C, and 6 purely internally heated models (TB models, Supplementary Table ST1 in Supporting Information S1), to address the role of the isothermal bottom assumption on the internal temperature.
Despite different boundary conditions, models achieve the same final potential temperatures, and comparable geotherms, which vary between ±25°C, and reproduce consistent maximum depletion degree, varying within ±2%, and comparable melt volumes, within ±12%, and more only in the reflective boundary walls model. This proves (a) the uniqueness of the mantle temperature-Melting degrees relation, while showing (b) the nonuniqueness of the temperature-Thermal boundary conditions. This supports the choice to present the models grouped by temperature ranges. We additionally discuss in Section 4.3 the limitation of this approach.

Results
Results are presented in two groups: in the first group, models reach mantle temperatures comparable to present day and melting is negligible, in the second, mantle temperatures are higher and melt degrees and volumes are larger. Within these groups, the thermal evolution differences are negligible and the regime is solely dependent on the lithospheric strength.
To compare the models we use the horizontally averaged mantle nonadiabatic internal temperature T i , the maximum potential temperature T P , the heat flow at the surface q = -k dT/dy, with the conductivity k = 0 c p , and the surface-averaged value , the root-mean-square (rms) velocity at surface , the surface mobility M = ∕ , and the total melt volume (area) in the top 250 km of the model V = ∫ 0 ∫ 250 0 ( ) dx dy, with L the width of the model space (Table 1). Additionally, we use the horizontally averaged lithospheric thickness h , defined by the depth of the isotherm T = 0.85 T i .

Convection Regimes With Low Depletion Degrees
Upper mantle models with internal heating rate comparable to present day, H ≤ 1 × 10 −11 W kg −1 and T b = 1,628°C show consistent mantle cooling throughout their evolution, evolving toward temperatures T P ≤ ∼1,430°C. Varying the yield strength in these models enables reproducing regimes that span from poor surface mobility (M ≤ 0.5), referred to as stagnant lid regime, to episodic mobility and mobile surface (M ≥ 1.0), with characteristics of present-day plate tectonics (e.g., Lenardic, 2018).
The models with high cohesion, 0 ≥ 20 MPa, develop a rigid lid above the convective mantle, with no major deformation and limited mobility. In models with cohesion 50 MPa (Figure 3a), surface velocities remain <0.5 cm yr −1 , as the lid develops downwellings, at locations where crust and lithosphere are shortened. Between these locations, large areas of thinned lithosphere and crust (Figure 3a, bottom panel, in grey) are found. The slow stretching of the thermal boundary allows thinning to be counterbalanced by cooling. This is shown by low surface heat flow of <40 mW m −2 measured above downwellings, and larger values, yet <∼80 mW m −2 , above stretching domains, indicated by divergent velocities (Figure 3a, top panel, in blue). The temperatures in these models remain consistently below the solidus, and no melting occurs.

Journal of Geophysical Research: Solid Earth
CAPITANIO ET AL.

10.1029/2021JB023911
10 of 26 Cohesions of 0 ≤ 15 MPa allow for the yielding of the lid and its increased mobility. This regime includes divergent and convergent margins with formation and recycling of lithosphere, akin to mid-ocean ridges and one-sided downwelling ( Figure 3b). The highest velocities reach ∼3 cm yr −1 , around convergent zones, and decrease to <1 cm yr −1 at sites of thickening. Lower yield strength favors thinning, decompression melting, and depletion; however, degrees remain small, F < 0.1, due to the relatively low mantle temperatures ( Figure 3b). Melt-depleted mantle is partly preserved in the lithosphere, due to its buoyancy. Heat flow varies largely between ∼150 mW m −2 and <10 mW m −2 , from divergent centers to downwelling areas. Models with lower yield strength 0 ≤ 5 MPa develop a similar regime to the previous model ( Figure 3c). Along downwellings, the down-going plate beneath a stationary upper plate accommodates convergence, therefore forming a subduction-like margin, at rates of 4-6 cm yr −1 .
Final geotherms of these models do not differ largely (Figures 3a-3c, right panel, solid red lines), and horizontally averaged viscosity (Figures 3b and 3c, right panel, solid black line) remain in all models <10 22 Pa s, in spite of local increase of maximum viscosity, due to minor depletion (dashed black line).
The evolution of the averaged mantle temperatures is similar for all models (Figure 4a), despite the different types of lid involvement in the convection. Temperature decreases at a rate of ∼100 ˚C/Gyr, attaining values of T P = ∼1,300°C by the end of the modeled time. In models with 0 ≥ 20 MPa, the averaged heat flow rapidly reaches values < 40 mW m −2 within ∼300 Myr due to the thickening of the conductive lid, then remains roughly constant throughout the modeled time ( Figure 4b). For lower cohesion values, the heat flow initially decreases, then increases to values >∼70-80 mW m −2 . The model's heat flow has episodic variations, occurring in shorter intervals for decreasing cohesion, reflecting episodes of lithospheric recycling. Surface rms velocities remain negligible throughout the evolution of the models with 0 ≥ 20 MPa (Figure 4c) and attain a mobility M < 0.5. These are features of a stagnant lid regime (see Lenardic, 2018) and references therein). The model with  (Moresi & Solomatov, 1998;Stein & Hansen, 2008). Models with cohesion ≤10 MPa, have similar features, although rapidly attain a mobility M = ∼1, due to the coupled overturn of lithosphere and mantle, which remain rather constant. Although some variations are observed, the consistent high mobility shows that these models are in a ML regime (Moresi & Solomatov, 1998;Stein & Hansen, 2008). Minor melt volumes are produced in these models (Figure 4e).
In summary, with mantle temperatures comparable to present day, melting is minor and the regimes depend mostly on the yield strength of the lithosphere. For cohesions above ∼20 MPa, the models reproduce sluggish to stagnant lid regimes. Instead, for values of cohesion <10 MPa, the models show features of a ML regime. For values in between, an episodic mode of convection emerges, where the stagnant lid is interrupted by short-lived events involving lithosphere mobility and recycling. Although the scaling of Earth-like stresses goes beyond the These results agree with what has been found in many different setups, from Cartesian to spherical, and from two-to three-dimensional (Crameri & Tackley, 2015;Langemeyer et al., 2018;Solomatov & Moresi, 2000;Stein & Hansen, 2008;Tackley, 1998Tackley, , 2000van Heck and Tackley, 2008).

Convection Regimes With High Depletion Degrees
With increased internal heat production H ≥ 2 × 10 −11 W kg −1 and T b = 1,628°C, higher mantle temperatures are achieved in the modeled period, resulting in large volumes of melt and degrees of melt depletion. The progressive thermochemical differentiation of the upper mantle into a stiffer lithosphere reduces mobility in low cohesion models, forcing the evolution of initial high mobility settings toward a low mobility lid (M ≤ 0.5) after ∼400 Myr. Yet, the mobile-to-stagnant lid evolution models have larger melt degrees and volumes than those in a stagnant lid.  The models with mantle cohesion of 0 ≤ 15 MPa allow mobility and the production of larger volumes of depleted mantle. However, larger volumes of dehydrated, stiffer mantle eventually suppress mobility, forming a rigid lid. Consequently, by the end of the modeled time, these models converge toward the same stable, thick, and stiff lithosphere, with low surface velocities and low heat flow. While these final characteristics are not different from the stagnant lid models, at higher strength, the low strength models differ from the latter for the higher depletion degrees and the formation of greater volumes of depleted lithospheric mantle. At the end of the modeled period, models with 0 = 10 and 15 MPa display surface velocities <0.5 cm yr −1 (Figure 5b and c) and a uniform, thick, and rigid lithosphere. Heat flow remains low, between ∼30 and 60 mW m −2 . The final configuration is indistinguishable from that achieved at higher strength (Figure 5a), although the melting is more pervasive and maximum depletion degree F is ∼0.36 and ∼0.44, for the two models shown, respectively, almost everywhere in the domain (Figures 5b and 5c).
The geotherm-viscosity comparisons in these models further illustrates the role of dehydration melting. While geotherms are comparable and hotter than the models presented in the previous section (Figures 5a-5c, right panel, red solid lines), the viscosity progressively increases from <10 22 Pa s to >10 24 Pa s (black solid and dashed line), leading to higher average lithosphere viscosities atop a hotter mantle.
The evolution of these models converges toward the same final configuration through different paths. Potential mantle temperatures increase at a constant rate of ∼200 ˚C/Gyr, reaching values T P > 1,600°C in all models (Figure 6a). At high lithospheric strengths, the models rapidly reach = ∼50 mW m −2 , then remain constant ( Figure 6b). Surface rms velocities are negligible throughout model evolution and mobility is <0.25 (Figures 6c  and 6d). Slow, rather constant velocities are indicative of a lithosphere prone to minor rifting only. Melt is produced at constant but at small rates throughout the simulations (Figure 6e). Models with low yield strength do not develop stable features and have a two-step evolution. In the first half of the modeled time, heat flow fluctuates between 90 and 70 mW m −2 , although the onset earlier, between ∼550 and 100 Myrs, and reach higher initial heat flow, 90 to 80 mW m −2 , for decreasing cohesion (Figure 6b). Surface rms velocities vary episodically and are up to ∼8 cm yr −1 , resulting in mobility variations between M = ∼1 and 0.5 (Figures 6c and 6d). Decreased yield strength allows for the faster onset of high mobility; however, this favors depletion and stiffening, and acts against sustained mobility. The negative feedback between yielding, depletion, and stiffening rapidly reduces mobility, which decreases to a constant value of ∼0.5, in the second half of the modeled time, while the heat flow drops to a lowered value of ∼50 mW m −2 in all models, then remains constant. This is further illustrated by the alignment of major melt production and mobility oscillations (Figures 6e and 6d). During episodes with increased mobility, features such as rifts and zones of downwelling form, where most of the depleted mantle is produced (Figures 5b and 5c). The growth in the volume of depleted mantle eventually stiffens the lithosphere above a rheological threshold, then mobility vanishes and ML-like features remain "frozen" in the lithosphere (Figures 5b and 5c, right panel).
The negative feedback leads models with high mantle temperatures to converge to the same endpoint: models with low strength evolve from an initial state of induced, short period episodicity to a poorly mobile, rigid lithosphere. The progressive thermochemical differentiation and stiffening reduces the difference between the models' surface mobility, although the melting and tectonic history are completely dissimilar for models in stagnant lid and ML regimes.

Convection Regimes and Depletion Degrees
Depletion degrees are largest in models with low cohesion, at all mantle temperatures achieved. In Figure 7 we show the depletion degree in the models in a cohesion-internal heat space. While detailed regime diagrams have been presented elsewhere (e.g., O'Neill et al., 2007), we show that only for large internal heat, H ≥ 2 × 10 −11 W kg −1 , depletion reaches F ≥ 0.36. Then, models with 0 ≥ 20 MPa reach low depletion degrees and produce negligible melting volumes (Figures 4-7). These models are all in a sluggish to stagnant lid regime (SL), for increasing cohesion, respectively. Lower cohesion allows larger melt production, reaching depletion degrees comparable to what is currently found at mid-ocean ridges, F = 0.1 ± 0.05, at low internal heat. The depletion degree increases substantially for internal heat ≥2 × 10 −11 W kg −1 , up to a maximum value of F = 0.44, favored by yielding.
Although an analysis of these regimes goes beyond the scope of the paper, all models with low strength display features of an episodic (EL) to ML, when internal heat is low, as presented elsewhere (see Lenardic, 2018, and

Journal of Geophysical Research: Solid Earth
CAPITANIO ET AL.

10.1029/2021JB023911
14 of 26 references therein). At higher internal heat, all models evolve through a phase of periodic mobility, first, before progressing to the same poorly mobile, thick, and rigid lithosphere ( Figure 6).
We do not constrain further the yield strength of the lithosphere, as yielding is strongly affected by the planform of mantle convection. Tractions beneath the lid increase in whole mantle convection, which helps, or sustains, rifting in the lithosphere (Lowman & Jarvis, 1996). On the other hand, mantle tractions decrease with increasing internal heat, making tectonics less viable (O'Neill et al., 2007). Additionally, model aspect ratios and dimensions affect the tractions (Guerrero et al., 2018;Lowman & Jarvis, 1996); however, these only shift the cohesion value at which the regime transitions (Langemeyer et al., 2018).
Additional tests show that the transition between regimes is consistent in upper mantle confined models, occurring at the same yield stress; however, in models with a lower mantle, the transition to a stagnant lid occurs at a yield stress that is a factor <∼2 higher than the upper mantle confined models (Supplementary Figure S4-S6 in Supporting Information S1).

Heat Flow Scaling and Implication for Thermal Evolution
In this section, we quantify the variation of lithospheric thickness and heat flow in the models, as a function of mantle internal temperature, and compare them to proposed scaling laws. The internal temperature, T i , the average lithospheric thickness h , and heat flow are determined by a numerical experiment.
To support the presentation of the results, we use parameterized convection scaling (e.g., Solomatov, 1995). The balance between internal heat and heat released through the surface is captured by the power-law relation between the surface heat flux and the vigor of convection: where the exponent β expresses the sensitivity of surface heat flux to convection. The scaling in (18) yields exponents that are weakly dependent on the choice of the reference viscosity (Korenaga, 2003).
Following Davies (1980), from Equation 18 the scaling of heat flow is then, the average thickness of the boundary layer h is Here, the parameters a and b are determined with readjustment of the heat loss equation to a reference heat flow, thickness, and temperature (e.g (Christensen, 1985): where q 0 = k T 0 /h 0 , T 0 = 1,400°C, h 0 = 111 km, and ( 0) is the internal viscosity at T 0 . Additionally, we define the depth of melting h sol , found setting T(y) = T sol (y), which depends on the solidus chosen (Katz et al., 2003), and show for reference the corresponding heat flow q sol and Nu sol .
The extrapolation of (18) to the early Earth has remained problematic and values of β > 0 or β < 0 have been proposed. Thermal boundary layer theory for isoviscous convection finds that β ≈ 1/3 and 1/4 for basally and internally heated fluid with free-surface boundary conditions, respectively (Turcotte & Oxburgh, 1967). Numerical modeling and scaling analysis extends this finding to cases with temperature-and stress-dependent rheologies, forming stagnant to ML regimes (Moresi & Solomatov, 1998), where β ≈ 1/3 and 1/4 were found for a stagnant lid to a transitional regime, and increasing viscosity, in basally heated convection (Solomatov, 1995). Lower values of β ≈ 0 to 0.15 show the lithosphere's thermal evolution independence on mantle temperatures, either due to a "surface effective viscosity" or to bending dissipation (Christensen, 1985;Conrad & Hager, 1999), thereby illustrating the control of visco-plastic rheology on the regime. When extrapolated to the early Earth, β > 0 must be discarded (see Christensen, 1985;Lenardic, 2006) and Korenaga (2003) proposed β = −0.15, when dehydration stiffening is considered and a thermochemical boundary layer forms. For the focus on the early Earth of this work, we do not aim at the fit of models and scaling laws, but rather use the scaling to support the models' outcomes presentation.
We first compare models with large strength, 0 > 15 MPa, which are all in a stagnant lid regime (SL). We plot all measurements for every time step (Figures 8a and 8b), starting from the time the heat flow has stabilized (see  Table 2.

Journal of Geophysical Research: Solid Earth
CAPITANIO ET AL.

10.1029/2021JB023911
16 of 26 Figures 6 and 4b). The models in an SL regime show an increasing thickness with decreasing internal temperatures for T i < 1,400°C and internal heat production H ≤ 2 × 10 −11 W kg −1 (Figure 8a, black, blue, and dark green circles). All these models follow the anticorrelated thickness-internal temperature trend, although showing a systematic thickness increase. This is likely due to the effect of the depleted mantle, which becomes increasingly thick. This trend clearly inverts for T i > 1,400°C (Figure 8a, brown circles). With increasing internal temperature the average lithosphere thickness increases, showing the effect of the large volumes of stiffer depleted mantle on the lithospheric thickness. The heat flow follows similar trends, with increasing surface heat flow for T i < 1,400°C (Figure 8b, black, blue, and green circles), the inverting at higher temperature, as shown by models with highest internal heat tested (Figure 8b, brown circles).
These measurements from the models follow the proposed scaling between T i and h and with a β > 0 for T i < ∼T 0 changing to β < 0 at higher temperatures, T i > ∼T 0 . The models with basal heating and with mixed heating follow a complex evolution aligning with trends defined by β = 0.3 and 0.25, and we do not attempt any further fit. At higher temperatures, the models are initially affected by the melt depletion, then follow the trend defined by β = −0.15, as suggested by Korenaga (2003). This work also provides an explanation for the systematic departure from the trends for intermediate temperatures (Figures 8a and 8b, green circles).
A compilation of all models shows how all regimes converge to a final stagnant lid configuration, as presented earlier (Figures 8c and 8d). Models are plot using time-averaged values in the last ∼400 Myr of the modeled time. While the SL regime models (circles) follow the trends shown before, the episodic and ML regimes (squares) deviate and the thickness and heat flow tend to be less dependent of internal temperature, when depletion is negligible. The models with negligible mantle depletion (black squares) align with the trend defined by the scaling with exponent β = 0, proposed for regimes controlled by a visco-plastic lithosphere. However, with increasing mantle depletion the difference between SL and E-ML models reduces. Models with a small melting production and mantle depletion (cyan squares) shows an intermediate behavior between the SL and the E-ML trends, whereas for T i > ∼1,300°C, the difference between the two endmembers regime reduces (light green squares), and for T i > ∼1,400°C the difference in the final configuration has vanished. This analysis shows the role of dehydration stiffening on the thickening and thermal stabilization of the lid, illustrating how this suppresses any difference in convection models' regimes at high temperatures. While the models initial stage reproduces SL or E-ML regime, these eventually converge toward a unique final configuration of a poorly mobile surface.

Controls on the Craton Formation
In this section, we compare available observations from Archaean cratons with the models, to test the role of depletion-dependent rheology in the formation of cratonic cores. The diagnostic features we focus on are the degree of melting, the volumes of depleted lithosphere, and the depth distribution.
Geochemical and petrological constraints on mantle temperatures, melt/depletion degrees and distribution are recovered using remnant basalts, picrites and komatiites in cratons, and complementary lithospheric mantle peridotites found in xenoliths (Griffin et al., 1999(Griffin et al., , 2003Herzberg et al., 2007;Lee et al., 2011;Pearson, 1999;Pearson & Wittig, 2008). Melting temperature of basalts infer Archaean mantle potential temperatures increasing between ca. 1,470 and 1,640°C (Herzberg et al., 2010), at rates >∼100 ˚C/Gyr, while higher potential temperatures are possibly implied by komatikic melts. Mantle potential temperatures steadily decrease since the Proterozoic at cooling rates of −50 to −100 ˚C/Gyr, from 1,600 to ∼1,450°C, then reach present-day mantle potential temperatures of 1,350 ± 50°C. These temperatures are reproduced by the models, with temperatures between 1,560 and 1,620°C and a heating rate of ∼170 ˚C/Gyr, for the largest internal heat, H = 3 × 10 −11 W kg −1 , and T b = 1,620°C, which fall in the range of the Archaean values. Models with moderate internal heat, H = 2 × 10 −11 W kg −1 , show temperatures clustered around ∼1,450°C, in agreement with those inferred for the end of the Proterozoic, while for H ≤ 1 × 10 −11 W kg −1 the models reproduce mantle temperatures in agreement with the Phanerozoic values. Models with H = 1-2 × 10 −11 W kg −1 allow reproduction of the cooling rates since the Proterozoic, varying from 130 to 0 ˚C/Gyr. These trends vary little among the different regimes. Although these temperatures can be achieved with a diverse range of boundary conditions and model geometry (Supplementary Figure S2 in Supporting Information S1), the results only depend on the temperature achieved, as these are sensitive to the melting.
High depletion degrees are characteristics of the cratonic cores and suggest that these formed in an environment with high, although ephemeral, mobility. Variable depletion degrees are constrained from the rock record, reaching maximum values of 0.3-∼0.45 (Lee et al., 2011) in Archaean cratons, then decreasing to <∼0.3 in the Proterozoic, reaching present-day values of ∼0.08 (Figure 9, color coded boxes). When compared to the models with temperature compatible with the Archaean mantle, only at lower yield strength, 0 < ∼30 MPa, melting degree reaches F > 0.30, while larger values, reaching F = 0.44, are reproduced with the lowest cohesion tested (Figure 9, red rim symbols). Similarly, models with internal heat H < 2 × 10 −11 W kg −1 yield temperatures compatible with the Proterozoic (Figure 9, green rim symbols); however, only models with 0 < ∼10-20 MPa achieve compatible depletion degrees of F ≤ 0.38. For mantle temperature comparable to the present day (Figure 9, blue rim symbols), values of F decrease to <0.26, and to ∼0.1 ± 0.05, for no internal heating (Figure 9, black rim symbols). These models recover realistic present-day temperatures and melting degrees (McKenzie & Bickle, 1988).
Additional constraints come from the comparison of depth-distribution of depleted mantle in the models and in the Archaean cratons, which support the idea of severe thinning in a rift-like environment (Lee & Chin, 2014). The models' cumulative composition of the top 250 km over the 7,920 km-wide domain at the end of the calculations (Figure 10) has two unique characteristics controlled by the yield stress: (a) larger volume of high-degree melting and (b) shallower depth distribution with decreasing lithospheric yield stress. Surprisingly, the variations due to potential temperature are secondary. The abundance of depleted mantle varies with increasing mantle temperatures; Figure 9. Maximum potential temperature T P versus largest depletion degree F MAX , from cratons observations and numerical models. Data for Archean (pink box), Proterozoic (green box) and Phanerozoic (blue box) from Herzberg et al. (2010) and Lee et al. (2011). Dashed line for the values inferred from komatiites. Models' potential temperature span the last 400 Myr of the modeled period.

Journal of Geophysical Research: Solid Earth
CAPITANIO ET AL.

10.1029/2021JB023911
18 of 26 however, for the highest yield strength, σ 0 = 50 MPa, maximum depleted mantle volume reached is <∼50%, and only ∼15% reaches higher degree melting F (Figures 10a, 10c and 10e). The cumulative composition at lower strength increases dramatically and depleted rocks make up from 60%-90% to 20%-60% of the lithosphere, for decreasing mantle temperatures. The vertical distribution of the depleted material is strongly affected by the yield stress. In high cohesion models and high internal heat (Figures 10a and 10c), the distribution reaches a maximum as deep as ∼80-90 km and has a smoother distribution. Instead, in the lower strength models, σ 0 = 5 MPa (Figures 10b and 10d), the distribution is skewed toward shallow depth, with pronounced peaks between 10 and 40 km. At lower internal heat, these characteristics are lost, and the distribution of depleted rocks is weakly depth dependent. Differences in models with values for σ 0 < 10 MPa are negligible (Figures 3 and 5). Although the volumes here are to be considered an upper bound (see above), the inference to a yield strength comparable to the present today remains valid.
In most of the Archaean cratons, the occurrence of peridodite with the largest degrees of depletion makes up to the largest component of the lithospheric volume and is largest at shallow depth (Griffin et al., 2003), further supporting the idea of a rift-like environment in a high mobility lid. Most of these cratonic keels have similar distributions, with largest volumes of depleted mantle with F ≥ 0.3 occurring at shallow depth, tailing off at ∼150 km depth (Pearson & Wittig, 2014). The Kaapval craton's keel shows pronounced stratification (Griffin et al., 2003;Kopylova & Russell, 2000), with a peak reaching 90% depletion found in the uppermost part (Figure 10g). The Slave craton has a ∼100 km-thick lithospheric keel, mostly (∼80%) made up of depleted mantle rocks, decreasing half way through its thickness (Figure 10g). Other domains, such as the Siberian craton, the Limpopo Belt, and the Lesotho terrane (Figure 10h), show an abundance of depleted mantle similar to the Kaapvaal and Slave keels', although their peaks in depletion occur at deeper levels. Volumes and shallow occurrence of depleted mantle lithosphere are most compatible with models with lower strength (Figures 10b, 10d and 10e), between 30 and 5 MPa. Finally, the composition of the Botswana and Gawler cratons indicates small volumes of highly depleted mantle in the keel and a smooth distribution with depth. Such small volumes are compatible with those in models with lower temperature in the mantle. Although speculative, these latter two cratons are younger than the other cratons (Handy et al., 2007;Smith et al., 2009), resulting in the secular cooling potentially biasing the cumulative distribution toward the characteristics of the cooler mantle models.
In summary, the preserved depleted keel distributions in Archaean cratons suggest a lithospheric yield strength low enough to allow high, although short-lived, surface mobility. While the volumes in the models are to be considered an upper bound (see above), degree of melting is strongly affected by the height of the geotherm above the solidus, maximized by thinning of the lithosphere, as suggested by the shallow emplacement depth. This supports the inferences to a ML and is contrary to the idea that rigid lithospheric keels are formed in a global stagnant lid regime (see Bédard, 2018;Stern, 2018).
A lithospheric strength low enough to allow for mobility is a condition necessary to reproduce the variety of tectonic features in the Archaean record. Features of vertical and horizontal tectonics are found in most cratons (Van Kranendonk, 2010). Archaean granite greenstone terrains and their deformation have been reconciled with a large-scale tectonic environment with negligible horizontal compression (Griffin et al., 1999(Griffin et al., , 2003Simon et al., 2007;van Hunen and Moyen, 2012;Van Kranendonk et al., 2007). Melting beneath very thin lithosphere is recorded in the Kaapvaal craton (Simon et al., 2007), with large volumes of depleted continental lithospheric mantle, formed at ∼3.5-3.2 Ga. Subsequently, short-lived subduction-like environments, at ∼2.9 Ga, are recorded (Simon et al., 2007). In the Pilbara craton (Van Kranendonk et al., 2007), a similar occurrence of mantle melting and depletion is recorded after ∼3.2 Ga, with arc-like magmatism, interpreted to reflect short-lived subduction and episodic rifting, whereas environments akin to rifting were widespread prior to this time, atop repeated mantle plume activity. Similar evolution is proposed for the Inukjuak domain, Québec (Caro et al., 2017), where evidence of Hadean recycling and Eoarchaean stabilization points toward a lithosphere with initial mobility, that later vanished.
Other mechanisms have been proposed to form cratonic cores and the structures they host. Compressional stresses are invoked to force thickening of lithosphere remnants into cratons (Beall et al., 2018;Wang et al., 2018). Similar structures may be also achieved during oceanic closure-like events by underthrusting and imbrication (Cooper & Miller, 2014), while plume impingement (Liu et al., 2021), rifting (Lee & Chin, 2014), and accretion of subduction arcs are invoked to explain high depletion degrees (see Lee et al., 2011, and references therein). These features are compatible with those generated in our models.

Speculations on the Geodynamic Regime of the Early Earth
Cratons currently cover a small percentage of the surface of the planet. Whether the mechanisms forming the cratons were relevant to the whole Archaean Earth is critical to the inferences on the global geodynamic regime. Here, we discuss the role of the craton-forming mechanism for the thermal evolution of the Earth.
The thermal evolution of the Earth is controlled by the balance of internal heat and the flux from its surface. In the parameterized approach, this is expressed by the relation Nu ∼ Ra β . When melting and dehydration stiffening are not considered in mantle convection models, the evolution follows the scaling with a β ≥ 0, ranging from values of β = 0.25-0.3 found for the stagnant lid regime and a range of temperature-dependent viscosity values (Christensen, 1985;Davaille & Jaupart, 1993;Moresi & Solomatov, 1998), to β = ∼0, when models have a lowered yield strength, and the lithosphere high mobility reduces the Nu-Ra dependence (Christensen, 1985). However, the extension of the scaling to higher mantle temperatures, comparable to the Archaean's, remains problematic. This scaling predicts a convergence toward a poorly ML regime, for higher temperatures, that is, no plate tectonics. However, while this is commonly invoked as the context for the early Earth tectonic remnants (see Bédard, 2018;Stern, 2018), it must be ruled out as a viable global heat-regulating regime as it predicts excess internal heat release, the "thermal catastrophe", which is not confirmed by the observations (Christensen, 1985;Davies, 1980).
Our results show that when mantle depletion and stiffening are considered, increasing mantle temperatures result in thicker lithosphere and reduced heat flow, aligning with the scaling defined by a negative β. At temperature comparable to present day, melting is negligible and the model reproduces the range of β from ∼0.3 to ∼0, that is from stagnant lid to ML regimes. However, for increasing temperatures, the mantle depletion volumes and degree increase, progressively narrowing the range of possible regimes, in spite of the low yield strength and the initial high mobility, converging to the same β = −0.15 trend. This trend results from the deepening of the geotherm -solidus intersection, allowing for melting in a thicker mantle layer, favoring the differentiation of a thicker lithosphere (Korenaga, 2003). The consequent stiffening enhances the control of a temperature-and depletion-dependent viscosity (Korenaga, 2006) over temperature dependence only, thereby breaking down the dependence of lithosphere thickening and heat flow's from internal temperature (Doin et al., 1997).
These results agree with other models proposed to solve the thermal paradox invoking a lithosphere substantially thicker than the TBL, a thermo-chemical boundary, supporting the idea that the mechanism proposed here is of relevance for the evolution of the whole Earth. Korenaga (2006) proposes the controls of dehydrated, stiffer lithosphere in a hotter Earth, thereby resulting in increasing lid thickness and viscosity with increasing temperature, still allowing reduced surface mobility. This yields a negative heat flow-internal temperature scaling, that is, β = −0.15, when dehydration stiffening becomes dominant. This configuration is compatible with a pervasive differentiation of the lithosphere and low heat flow (Figures 11a and 11e). However, as internal heat is lowered, mixed behavior emerges, with domains of lowered heat flow intervening with mobile environments where heat flow is higher (Figures 11a and 11d). This configuration agrees with that proposed in several studies (Lenardic, 2006;Lenardic et al., 2003Lenardic et al., , 2005 where different heat flux domains may modulate the integrated heat flow. This configuration helps mitigating the paradox of the Archaean thermal evolution (Lenardic, 2006), as well as relaxing the constraints proposed by the model of Korenaga (2006), which may lead to excessive stiffening and a shutdown of surface motions.
A limit of the modeling approach chosen here is the lack of latent heat in the conservation of energy, whereas this is only accounted for in the constitutive laws. Melt extraction and heat advection may have critically mitigated the excessive heating of the mantle (Moore & Webb, 2013) and facilitated mobility (Jain et al., 2019;Lourenço et al., 2018;Rozel et al., 2017) in the early Earth. However, while these models account for the inherent cooling of the mantle, the stiffening of the residual, cooler mantle is not considered. Here, we have shown that considering the impact of dehydration and cooling of the residual mantle, due to melting, leads to viscosity increase of large volumes of mantle. It remains to be tested whether the stiffening of the thick portion of lithosphere, where melt is extracted, would suppress further volcanism, stabilizing the lid, and suppress mobility providing rigid keels to squishy crusts, as shown in our models.
Another solution to the thermal paradox invokes intermittent subduction, extrapolating the Phanerozoic record to the Archaean (Silver & Behn, 2008), to reduce the predicted heating of the early Earth mantle. While we do not rule out subduction in our models, we do not support an intermittent fully fledged plate tectonic scenario to have existed throughout Earth history. The transition from a thermo-chemical boundary layer to a thermal boundary layer is associated with the emergence of a ML, therefore it may have triggered a major tectonic transition and the emergence of plate tectonics. Figure 11 shows that models with present-day mantle temperatures reproduce ML-like stable plate margins and heat flow distributions (Figures 11a-11c) and no melting. Instead, at higher mantle temperatures, mobility is hampered by the formation of larger mantle depleted volumes progressively embedding in the thermochemical layer (Figures 11a, 11d and 11e), while heat flow reduces as a consequence of its thickening. Other studies have illustrated the viability of an episodic regime for the early Earth (Höink et al., 2013;Jellinek & Jackson, 2015) and how increasing internal heat production (O'Neill et al., 2007;Stein et al., 2013) and depleted lithosphere rheology (Weller & Kiefer, 2020) favor a transition to a stagnant lid regime. However, the progressive stiffening of the whole lithosphere in our models makes this process irreversible, when conditions for the formation of cratons are met.

Journal of Geophysical
Constraining further when the transition might have occurred on Earth remains speculative. This is controlled by the dehydration stiffening and its waning, therefore by the yielding of the lithosphere. Scaling for the models' plane-layer geometry and aspect ratio imply mean temperatures are achieved with a much higher internal heat values in an Earth-like spherical shell (O'Farrell et al., 2013;O'Farrell & Lowman, 2010), although comparable temperatures and melt degrees between models and observations support the validity of the our outcomes. Additionally, a rearrangement to mantle convection to whole mantle and to wider aspect ratios might have critically increased the shear tractions beneath the lithosphere (Lowman et al., 2001;Lowman & Jarvis, 1996), thereby favoring yielding and melting, supporting, or forcing, the regime transition.

Limitations of the Modeling Approach
The comparable melt degrees between models and observations support the validity of the models' outcomes, although a note of caution is in order. The choice of the model geometry and aspect ratio, and the uncertainties on the early Earth mantle, do not allow tight constraints on realistic temperature evolution.
The choice of Cartesian two-dimensional modeling space is common in studies focusing on Precambrian lithospheric processes (e.g., Gerya, 2014) and in scaling analyses (e.g., Sleep, 2015); however, when applying the outcomes to the Earth, the planar 2D model space yields hotter geotherms. This setup greatly underestimates the internal heat required to achieve specific mantle temperatures, when compared to spherical geometries (O'Farrell et al., 2013;O'Farrell & Lowman, 2010). Additionally, the models' aspect ratio may strongly affect the temperatures and heat flow (Lowman et al., 2001); however, for a ratio above 6 (Lowman et al., 2001) and for large Ra the mean temperatures converge (O'Farrell & Lowman, 2010). We have neglected the role of three-dimensional structures, such as strike-slip tectonics (Langemeyer et al., 2021). However, comparisons of spherical and Cartesian models show that the tectonic regimes characterized by divergent and convergent margins are consistently reproduced in both geometries (Crameri & Tackley, 2014).
The choice of thermal boundary conditions may also limit the inferences on the early Earth. Mantle temperatures with an insulating bottom boundary are lower, ∼30%, than the isothermal boundary (Lowman et al., 2001). In this context, the lower mantle has a role in the thermal evolution (Supplementary Figure S4-S6 in Supporting Information S1), which is not fully explored here. We do not model plutonic emplacement and heat advection; however, our models are to be considered complementary to the squishy lid and plutonic-squishy lid regimes (Lourenço et al., 2018;Rozel et al., 2017): these studies focused on the role of the melt weakening, assuming variable to complete plutonism, whereas here we focus on the residue stiffening that follows melt extraction. The melt extraction and heat advection would contribute to further cool the residual mantle and stiffen the lid (Phipps-Morgan, 1997), therefore the implications for the formation of craton might be negligible.
The implications for the thermal evolution may also be affected by the lack of latent heat in the energy equation formulation, although it is accounted for in the rheology. This likely affects the volume of depleted mantle and might have a feedback on the heat-pipe mode (Moore & Webb, 2013). This requires a self-consistent melt transport model, otherwise done with ad hoc assumptions (Moore & Webb, 2013;van Thienen, 2004), which remains computationally challenging for large-scale models, to date.

Conclusions
Modeling upper mantle convection in 2D with internal heat values ranging from the present day to those of the hotter early Earth provides a tool to investigate the tectonics and geodynamic regimes that may have operated through Earth's evolution. We capture the role of melt extraction and associated cooling and stiffening implementing a temperature-and depletion-dependent rheology. At mantle potential temperatures similar to Phanerozoic, melting, melt extraction, and rheological stiffening of the mantle are negligible. Then, the Earth's regime and the viability of plate tectonics depends on the ability of the thermal boundary layer to yield, forming plate margins, such as subduction zones and ridges. In contrast, when internal heat production is higher and mantle temperatures reach over 200°C higher than the present day, melting increases substantially, leaving large volumes of depleted, stiffer mantle residue at shallow depth. In models with low lithospheric strength, yielding favors lithospheric mobility and thinning, with widespread rifting and localized subduction-like downwellings, where melting and larger depletion occur. In turn, larger depletion causes lithospheric mantle stiffening, which resists further yielding and mobility. This negative feedback initially induces oscillations in the mobility and then suppresses it completely. Despite cohesion and temperature conditions favorable to a ML regime, a rigid lithosphere forms. This end result is not different from the rigid lithosphere in the stagnant lid regime; however, the melting histories and depleted mantle volumes differ substantially, providing a diagnostic tool to probe Earth's cratons. When compared to the models, the record of shallow large melting and depletion observed in cratons is best reconciled with a tectonic environment in which cratons are formed as a consequence of lithospheric yielding, whereas the record of reworking, recycling, episodic subduction, and rifting support a high mobility lithosphere, eventually stabilizing to become a thick depleted cratonic keel. The models show a thermal regime transition between the