# Growth of continental crust and subduction of lithosphere in the Hadean revealed by geochemistry and geodynamics.

Alexander V Sobolev, Adrien Vezinet, Aleksandr V Chugunov, Charitra Jain, Stephan V Sobolev, Evgeny V Asafov, Valentina G Batanova, Nicholas T Arndt, Leonid V Danyushevsky, John W Valley, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Alexander V Sobolev, Adrien Vezinet, Aleksandr V Chugunov, Charitra Jain, Stephan V Sobolev, et al.. Growth of continental crust and subduction of lithosphere in the Hadean revealed by geochemistry and geodynamics.. 2024. hal-04848178

HAL Id: hal-04848178

https://hal.univ-grenoble-alpes.fr/hal-04848178v1

Preprint submitted on 19 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Growth of continental crust and subduction of lithosphere in the Hadean revealed by geochemistry and geodynamics.

**Alexander Sobolev**

alexander.sobolev@univ-grenoble-alpes.fr

University Grenoble Alpes https://orcid.org/0000-0002-1997-2032

**Adrien Vezinet**

University Grenoble Alpes https://orcid.org/0000-0001-9975-2760

**Aleksandr Chugunov**

University Grenoble Alpes

**Charitra Jain**

GFZ German Research Centre for Geosciences https://orcid.org/0000-0002-0562-2568

**Stephan Sobolev**

GFZ-Potsdam

**Evgeny Asafov**

Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Sciences

**Valentina Batanova**

ISTerre

**Nicholas Arndt**

University Grenoble Alpes

**Leonid Danyushevsky**

Norris Scientific

**John Valley**

University of Wisconsin, Madison

**Alina Koshlykova**

Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Science

## Article

### Keywords:

**Posted Date:** November 1st, 2024

**DOI:** https://doi.org/10.21203/rs.3.rs-3601806/v1
**License:** © Ⓟ This work is licensed under a Creative Commons Attribution 4.0 International License.

Read Full License

**Additional Declarations:** There is NO Competing Interest.
Growth of continental crust and subduction of
lithosphere in the Hadean revealed by geochemistry
and geodynamics

2

4

6

Adrien Vezinet¹*†, Aleksandr V. Chugunov¹*†, Alexander V. Sobolev¹*†, Charitra Jain²,
Stephan V. Sobolev²,³*, Valentina G. Batanova¹, Evgeny V. Asafov⁴, Alina N. Koshlyakova⁴,
Nicholas T. Arndt¹, Leonid V. Danyushevsky⁵, and John W. Valley⁶

8

**Affiliations:**

10

¹ Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel,
ISTerre; 38000 Grenoble, France.

12

² GFZ German Research Centre for Geosciences, Geodynamic Modeling Section; Potsdam,
Germany.

14

³ University of Potsdam, Institute of Geosciences; Potsdam, Germany.

⁴ Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Sciences; Kosygina str. 19, Moscow 119991, Russia.

16

⁵ Friendly Solutions; Australia.

18

⁶ WiscSIMS Lab, Dept. of Geoscience, Univ. of Wisconsin; Madison, WI, 53706, USA.

20

* Corresponding authors: adrien.vezinet@univ-grenoble-alpes.fr, aleksandr.chugunov@univ-grenoble-alpes.fr, alexander.sobolev@univ-grenoble-alpes.fr & stephan.sobolev@gfz-potsdam.de

22

24

† These authors contributed equally to this work

26
28

30

32

34

36

38

40

42

44

The rates of continental crust growth and recycling on early Earth remain controversial because materials in the ancient crust and mantle have been altered, or even erased, by ongoing geodynamical processes. Melt inclusions in minerals are pockets of magma that have been trapped and shielded from the external environment. Where found within Archean high-Mg olivine – the first mineral to crystallize in mantle-derived melts – these inclusions can provide an unaltered glimpse of the geochemical state of the early Earth’s mantle. We discovered an unprecedented unradiogenic Sr mantle source component ($^{87}$Sr/$^{86}$Sr=0.69932±0.00024, 95% c.i. here and below) in mantle-derived melts trapped in olivine from ca. 3.27 Ga komatiitic lava flows (Barberton Greenstone Belt, South Africa). This component shows a 4.31±0.19 Ga model age combined with significant chemical fractionation (Nb/U=36.9±1.5, Ce/Pb=16.7±1.1) translating to extraction, by the late-Hadean, of 80%±16% of mass of present-day continental crust assuming whole mantle processing. That agrees with the results of our geodynamic models explaining the Nb/U and Ce/Pb data by the production of 40 to 70% of the present-day continental crust mass during the Hadean in an oscillating mobile-lid tectonic regime with several tens of million years long periods of massive subduction induced by mantle plumes.
46

48

50

52

54

56

58

60

62

64

66

68

Evolution of the early Earth and dynamics of its geochemical reservoirs are not well understood due to the near absence of Hadean (> 4.0 Ga) rocks and minerals preserved at the present-day Earth's surface¹,². Consequently, the rates of continental crust growth and recycling on the early Earth are still debated: estimates for the fraction of continental crust formed at the Hadean-Archean transition range between 0% and 100% of the volume of today's continental crust³. Several approaches have been used to provide estimates of continental crust volumes during the Hadean and Archean³⁻⁷. Studies assessing the degree of mantle depletion⁴,⁶,⁷ are viewed as the most reliable since they account for both crustal extraction and recycling³. Most of these studies rely on the bulk rocks' Sm-Nd, Lu-Hf and U-Pb isotopic systems³,⁴,⁶. This approach, however, suffers from significant problems, such as known mantle heterogeneity on isotope ratios and systematic differences in these ratios of mid-ocean ridge basalts (MORB) and oceanic islands basalts (OIB) mantle sources⁸. In addition, the difference in isotopic effects of extraction and recycling of continental and oceanic crust or even fractionation effects at the magma ocean stage is not easily recognized⁴. Also, most ancient rocks are altered, which has been shown to affect Sm-Nd and Lu-Hf isotopic systems⁹,¹⁰. Further, whatever is the volume of continental crust on early Earth, the tectonic regime accounting for it remains debated, with suggestions being stagnant-lid, plutonic-squishy-lid, episodic-lid, mobile-lid, ridge-only and heat-pipe regimes¹¹,¹².

Among the varied geochemical proxies that have been used over the past decades to decipher the production and recycling of continental crust, Nb/U and Ce/Pb in fresh lavas or glasses have received the most attention. These ratios are regarded as canonical because (i) they do not fractionate when the mantle rocks melt and (ii) are similar for present-day uncontaminated MORB and OIB⁸,¹³,¹⁴ suggesting mantle homogeneity on these parameters. On the other hand,
70

72

Nb–U and Ce–Pb fractionate to give higher Nb/U and Ce/Pb values in the restite when hydrous mafic or ultramafic rocks melt in the crust or mantle, i.e., during the production of the felsic magmas (or their parentals) that form the bulk of continents¹³. This observation makes Nb/U and Ce/Pb ideal proxies to assess both the production of continental crust and recycling of restites through time⁴.

74

76

78

80

82

Three aspects limit the accuracy of crustal growth curves built from Nb/U and Ce/Pb ratios measured in mantle-derived rocks. First, the timing of the generation of these ratios in the mantle source is unknown and could be much older than the eruption age of measured magmas. This means that mantle-evolution curves based on Nb/U or Ce/Pb ratios in mantle-derived melts¹⁵ provide only lower estimates of the rate of crustal growth. Second, subduction of continental crust or/and unmelted oceanic crust to the source of mantle-derived magmas will decreases the Nb/U and Ce/Pb ratios of these melts, leading to underestimation of the extent of continental crust production. And third, the effect on these ratios of post-emplacement alteration, crustal assimilation and/or metamorphic overprinting is difficult to identify and quantify, leading to inaccurate interpretations⁹,¹⁰,¹³.

84

86

88

90

Since olivine is the first mineral to crystallize in mantle-derived magmas, melt inclusions in high-Mg olivine crystals—those with high Mg# or Fo (>=Mg/Mg+Fe molar%)—provide the most reliable information about parental melt compositions and the evolution of their mantle sources through time¹⁵⁻²⁰. This stems from the capacity of melt inclusions to acquire the elemental and isotopic compositions of the host magma and preserve this information long after crystallization¹⁶. Consequently, melt inclusions in olivine (Figure 1) may provide a reliable view of the geochemical composition of the mantle source of melts.
The novelty of this study's results is twofold. First, we discovered a unradiogenic Sr mantle component previously unknown on Earth with unequivocal fractionated Nb/U and Ce/Pb values. These ratios are more uniform in the present-day mantle than Nd, Hf, and Pb isotopes⁸,¹³,¹⁴ and explicitly mark the extraction of continental crust. Our Sr isotopes and Nb/U and Ce/Pb data came from pristine komatiite melt inclusions in olivine protected by the host mineral from alteration and suggest a significant event of continental crust extraction in Hadean time. Similar data on Nb/U and Ce/Pb of rocks are known only starting from 3.5 Ga and are compromised by their severe alteration. Second, we conducted advanced geodynamic modeling of the early Earth dynamics and extraction of the continental crust, coupled for the first time with the evolution of Rb-Sr isotope system and Nb, U, Ce, and Pb trace elements. This allows us to link geochemical data to the physical processes that might have been responsible for the observed evolution of chemical composition. In particular, using our models constrained by geochemical data we can discriminate between tectonic regimes that could have been active in Hadean and Eo-Archean time and suggest the preferred scenario of an oscillating mobile-lid tectonic regime with several tens of million years long periods of massive subduction induced by mantle plumes.
<figure><img src="image_3.png" /><figcaption>Figure 1: a. Olivine cumulate from Weltevreden Formation (sample 1528C) showing that although these cumulates are significantly altered, they still contain preserved unaltered olivine cores (plane polarized light). b. Partially crystallised natural (not annealed) olivine-hosted melt inclusion before homogenisation at high temperature. Inclusion consists of glass and acicular calcic pyroxene. c. Heated and quenched olivine-hosted melt inclusion consisting of glass, shrinkage bubble and spinel (dark crystal).</figcaption></figure>

108

110

112
114

Weltevreden komatiites

116

118

120

122

124

We determined the chemical and isotopic compositions of olivine-hosted melt inclusions in samples from komatiitic flows from the *ca.* 3.27 Ga<sup>ref. 21,22</sup> Weltevreden Formation (Barberton Greenstone Belt, South Africa, Fig. S1). Whole-rock isotopic investigations indicate that the mantle sources of these komatiites experienced mantle depletion before eruption, possibly during the Hadean<sup>23</sup> (Figure S2). Here, using Sr-isotope signatures and trace-element contents measured directly in the homogeneous glass of melt inclusions (Figure 1c), we provide solid geochronological constraints for the timing of both the isolation of mantle domains from the Bulk Silicate Earth (BSE, also referred to as primitive mantle, PM) and the acquisition of canonical ratios indicative of continental crust extraction.
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Contents of selected major, minor and trace elements and Sr isotopes in komatiite melt inclusions (MI) and host olivine crystals. Fo = 100Mg/(Mg+Fe) molar. Filled circles are those with a corresponding Sr-isotope analysis: blue for main group of melt inclusions, yellow for Sr-radiogenic melt inclusions, orange for Sr-unradiogenic melt inclusions (see below). Light grey circles are inclusions with no Sr data or rejected Sr isotope data because of Rb/Sr ratio over maximum value of used reference materials (0.0244) or uncertainty (2SE) of ⁸⁷Sr/⁸⁶Sr value > 0.0015 (SOM). Expected seawater²⁴. Errors bars are average individual uncertainty at 95% confidence interval (c.i.). Outlines on figures 2a,b indicate samples from the same flows.</figcaption></figure>

126

128

130

132
134

Melt inclusions in komatiite olivine phenocrysts

136

138

140

142

144

More than 350 olivine-hosted melt inclusions were selected and prepared for geochemical analyses. The analytical workflow, detailed in Methods, was as follows: (i) annealing of natural partly crystallized melt inclusions at high temperature and quenching to produce homogenous glass, (ii) electron microprobe analysis to determine major and minor elements in melt inclusions and host olivine, (iii) Raman-spectroscopy to determine H₂O contents of melt inclusions, and (iv) laser ablation split stream (LASS) ICP-MS analysis for Sr-isotope compositions and trace element contents of melt inclusions. Not all inclusions could be analyzed through the entire analytical workflow since LASS analysis requires a diameter ≥ 38 µm (optimal laser beam size) and only 195 inclusions satisfied this criterion. In addition, 19 more inclusions were analysed in single stream mode for trace elements.

146

148

150

152

154

156

**Figure 2a** and **b** present the composition of virtually unzoned olivine hosts of different isotopic groups of melt inclusions (see below) from five separate flows. The Fo and trace element contents (shown NiO and Cr₂O₃) of olivine differ both between and within flows. This variability suggests that the studied olivine crystals are, in fact, antecrysts that grew and trapped their melt inclusions in different places of a complex magmatic plumbing system. The TiO₂ contents of Weltevreden melt inclusions show a large variation, from ca. 0.14 wt. % to 0.21 wt. %, which is inversely correlated with host-olivine Fo contents (Figure 2c). This correlation follows a trend of ca. 33 wt% olivine extraction from the common parental melt (estimated from ingrowth of Ti from the melts trapped in the most Fo olivine) indicating that TiO₂ contents in the melt, as well as Al₂O₃ and CaO contents (not shown), are governed by olivine crystallization and extraction from parental melts with similar contents of these elements. Other proxies, such as Cl/Ti (Figure 2d-f), or K₂O/TiO₂ (Figure S3b-e), should remain constant in the melt during olivine
crystallization because of negligible contents of these elements in olivine. Clearly it is not the case. We attribute the large range of Cl/Ti and K₂O/TiO₂ for the same composition of host olivine to assimilation into the komatiitic melts of small (ca. 1-2wt%) amounts of seawater-derived ultra-saline brine before emplacement at ca. 3.27 Ga (see Methods): a process common for modern submarine basaltic melts²⁵,²⁶ or larger amounts (e.g. 10%) of seawater-altered depleted serpentinite. Values of Cl/Ti can therefore be used as a proxy for such contamination. A positive correlation is observed between Cl/Ti and Rb/Sr (Figure 2e) suggesting gain of Rb as well as Na, K, Pb, H₂O and Sr during contamination (Figure S3). To avoid samples significantly affected by this process, we considered only melt inclusions with low Rb/Sr for geological interpretations in this study (Figure 2e). A maximum value of 0.0244 was selected for Rb/Sr because it corresponds to the highest Rb/Sr value of the reference glasses used for LASS analysis (KL2-G, Figure S4). Unlike Rb/Sr, other proxies, such as Ce/Pb and Nb/U, show negative correlation with Cl/Ti, suggesting that the higher values measured for these two ratios do not result from the contamination of the komatiitic melts before entrapment.

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Sr-isotope of Weltevreden olivine-hosted melt inclusions. Both measured (panel a) and initial (panel b) ⁸⁷Sr/⁸⁶Sr values show low p(χ²) (<< 0.05) indicating the lack of statistical homogeneity at the precision obtained for these Sr-isotope analyses. Three statistically</figcaption></figure>

172

174
homogenous groups can be distinguished in the data shown in panel b: unradiogenic in orange, radiogenic in yellow and neutral (main) in blue (see Methods for statistical treatment).

## Sr isotopes and geochemical tracers in melt inclusions

The Sr-isotope signatures of Weltevreden melt inclusions show a large range of both measured and initial (back-calculated at 3266 Ma) ⁸⁷Sr/⁸⁶Sr values supported by $p(\chi^2)$ values << 0.05 indicating statistical heterogeneity within the whole dataset (Figure 3a,b). The effect of minor contamination of Weltevreden melts by seawater-derived brine or depleted serpentinite on their Sr isotope composition is within the analytical uncertainty, as seen by the lack of a significant correlation of ⁸⁷Sr/⁸⁶Sr value of melts and Cl/Ti, the chemical proxy of such contamination (see Figure 2h).

<figure><img src="image_6.png" /><figcaption>Main geochemical features of 4 replicate analyses of a single melt inclusion</figcaption></figure>

<figure><img src="image_7.png" /></figure>

**Figure 4.** Sr-isotope and trace element signatures of the Weltevreden olivine-hosted melt inclusions of Sr-unradiogenic group. The Sr model age (see Methods) of the unradiogenic group (panel a) indicates mid-Hadean extraction from the Bulk Silicate Earth (BSE) while the Nb/U (panel b) and Ce/Pb (panel c) values of these-inclusions are significantly higher than those of primitive mantle (PM) or BSE, interpreted as reflecting restites after extraction of continental crust in the presence of H₂O during the Hadean. Panels d, e and f show results obtained on 4

180

182

184

186

188

190

192
replicate analyses of a single melt inclusion (MI) 2217-48-ol65 (belongs to of Sr-unradiogenic group) showing the robustness of our analytical protocol. Uncertainty bars are reported at a coverage factor of 2 (i.e., 95% confidence interval), (see Methods and Fig. S5 to S7).

A statistical analysis of our Sr-isotope dataset reveals three statistically homogenous groups based on their initial ⁸⁷Sr/⁸⁶Sr (see Figure 3b and Methods for details). The means of initial ⁸⁷Sr/⁸⁶Sr values of these groups differ with confidence of over 99.9% (Supplementary Table 6). Out of these three groups, one, comprising 14 analyses, shows an extremely depleted (⁸⁷Sr/⁸⁶Sr)initial signature yielding a (⁸⁷Sr/⁸⁶Sr)initial weighted mean of 0.69932 ± 0.00024 (95% confidence interval, MSWD=0.88, p(χ²)=0.57, orange analyses in Figure 3b). Olivine containing unradiogenic inclusions come from all studied flows (14 of total 137, thus 10%), but mostly (8 inclusions out of 14) from Keena's flow 2 (samples 1523B and 2217, Figure 2a,b, Fig. S1), where they compose 14% of 59 measured for Sr isotopes inclusions. The mean initial isotopic ratio of unradiogenic inclusions transposes to a Sr-model age, assuming no Rb in the source, of 4.31 ± 0.19 Ga (95% confidence interval, p(χ²)=0.57, Figure 4a) which we interpret to indicate that some components that melted to form Weltevreden komatiitic melts were geochemically isolated since the Hadean. Owing to the large radius of inclusion 2217-48-ol65, belonging to the unradiogenic group, four Sr-isotope & trace elements replicates were conducted. As demonstrated in Figure 4d,e,f and Supplementary Table 3d, all four replicate analyses form a statistically homogenous population for Sr isotopes and canonical ratios of interest, hence further validating the analytical approach developed in this study.

## Continental crust extraction in the Hadean

Previous studies have used the geochemical proxy Nb/U to monitor the extent of continental crust extraction⁴. The weighted mean of the Nb/U values of the unradiogenic Sr melt inclusions is 36.9±1.5 (MSWD=0.71, n=14, Figure 4b). This value is significantly higher than all estimates
of the PM/BSE Nb/U values which lie between 26.7 and 32.4<sup>ref. 8,27-29</sup>. A similar observation can be made for the mean Ce/Pb value 16.73±1.12 (MSWD=2.3, n=13, Figure 4c) of this unradiogenic group, which is significantly higher than the PM/BSE value of 9-11.2<sup>ref. 8,27-29</sup>. The four replicate obtained on inclusion 2217-48-ol65 show the same tendency (Figure 4d,e,f). The geochemical effect of Ca and Mg perovskite cumulation proposed for the source of Weltevreden komatiites<sup>30</sup> to explain moderate decoupling of Sm-Nd and Lu-Hf isotopic systems in their composition will produce Nb/U ratios in the source that are lower, not higher than in BSE/PM as observed in melt inclusions due to higher partition coefficient for U than for Nb in suggested mixture of Ca and Mg-perovskite<sup>31</sup>. In addition, similar and even much more significant decoupling of εNd and εHf is common in highly depleted abyssal and ophiolitic peridotites because of "following radiogenic ingrowth as a result of preferential fractionation of Lu/Hf during partial melting"<sup>32</sup>. The main mantle source component of studied komatiite was highly depleted mantle harzburgite as it follows from their contents of incompatible elements (including Sm-Nd and Lu-Hf), which are lower in parental melt than in BSE mantle<sup>15</sup>. Thus, the decoupling of εNd and εHf in their sources can be explained by the mantle restite of Eoarchean-Hadean age.

Extraction of continental crust from mafic crust and/or peridotite in the presence of H₂O followed by recycling of restite back to the mantle leads to depletion in Rb compared to Sr, Pb compared to Ce and U compared to Nb in the mixed deep mantle source. This explains the low Rb/Sr and high Ce/Pb and Nb/U values measured in mantle-derived melts¹³. However, the presence of subducted continental crust and unmelted mafic crust in the mantle source may significantly affect all these elemental ratios. Fortunately, the amount of such crust in the recovered Hadean component in Weltevreden melt inclusions is shown to be negligible by the exceptionally low contents of ⁸⁷Sr of the unradiogenic component. Thus, high Nb/U and Ce/Pb
242

244

246

248

250

252

254

256

258

ratios of unradiogenic-Sr melt inclusions are fully attributed to the production of continental crust. However, the estimated Sr isotope model age of the mixed source could reveal younger dates by few tens of million years, which are well within our reported ± 0.19 Ga uncertainty of the model age (see Fig. S7 and Methods).

The fraction of extracted continental crust produced from a specific domain of the mantle can be evaluated following the mass balance approach that has been adopted in previous studies⁴ using our Nb/U values and Sr geochronological constraints together with the estimated Nb and U contents of the mantle⁸,²⁷⁻²⁹ and continental crust³³⁻³⁶. First, the fraction of crust extracted from any volume of processed mantle can be estimated. Using the Nb/U weighted average of 36.9 ± 1.5 obtained on inclusions from source-rocks differentiated in the mid-Hadean, we estimate the fraction of extracted crust to be 0.43%±0.09% of the processed mantle. Assuming now 2.17×10²² kg as the mass of the present-day continental crust, 4.01×10²⁴ kg as the mass of the whole mantle, the fraction of continental crust extracted in the Hadean, would be 80%±16% (Figure 5, orange dot) of the mass of the present-day continental crust. Similar calculation for individual melt inclusion 2217-48-ol65 yields 72%±24% (Figure 5, hatched orange dot). Such a high fraction, being the maximum value because it assumes the involvement of the whole mantle, still provides unprecedented evidence in support of earlier models arguing for massive continent production and recycling very early during Earth's history³,⁶,³⁷.
<figure><img src="image_8.png" /><figcaption>**Figure 5:** Compilation of previously published continental growth curves³ and new estimate for the Hadean from this study (orange circles). These estimates have been obtained using measured average Nb/U value in the Sr-unradiogenic group large (orange dot) and individual inclusion 2217-48-ol65 (orange hatched dot) and the average of Nb and U values from four mantle models,⁸,²⁷⁻²⁹ and four crustal models³³⁻³⁶ assuming processing of whole mantle mass. Red curves are models in agreement with our results indicating continental crust production and subsequent recycling of restites in the Hadean. Uncertainty bars are reported at a coverage factor of 2 (i.e., 95% confidence interval). Light blue field shows the range of continental crust fractions produced during Hadean which according to the geodynamic models fit best geochemical data (see below Geodynamic modelling section).</figcaption></figure>

Our new results interpreted as indicating that a large fraction of continental crust was extracted, and its residue was recycled into the mantle very early in Earth history is consistent with the composition of ancient zircons³⁸,³⁹, the radiogenic isotope compositions of bulk rocks⁴⁰ and numerical modelling (see Geodynamic modelling). Elevated Nb/U and Ce/Pb values obtained in deeply-sourced Weltevreden komatiitic melts indicates efficient, vertical, transfer (subduction of oceanic crust residue after partial melting producing felsic crust) between the upper and the
278

280

282

lower mantle in the Hadean. Although the geochemical results presented here indicate significant continental crust production in early Earth time, the tectonic regime accounting for the formation of this crust and actual mass of mantle processed cannot be extracted from such a geochemical dataset. To complete the picture, geodynamic modelling, within which $^{87}$Sr/$^{86}$Sr, Nb/U and Ce/Pb geochemical proxies were monitored in order to identify the configuration/regime best-reproducing the measured geochemical information.

## Geodynamic modelling

284

286

288

290

292

294

296

298

We used the mantle convection code StagYY$^{41}$, which has been extensively used over the last decades for the modelling of coupled core-mantle-crust evolution of rocky planets. With a 2D spherical annulus domain$^{42}$, our models generate both basaltic and felsic melts$^{43}$, include cooling of the core, and use pressure- and temperature- dependent water solubility maps for different mantle minerals$^{11}$. For this study, we improved previous models$^{11,43}$ by including the effect of water on the density of mantle materials, by incorporating a composite rheology (diffusion creep and dislocation creep proxy) for the upper mantle based on experimental data, and by initializing the frictional strength of the early Earth's oceanic lithosphere based on experimental data.

Following approaches from previous studies$^{44,45}$, we included the evolution of Rb-Sr isotopes as well as Nb, U, Ce and Pb trace elements in these geodynamic models (see Methods for further details). We computed a series of models by varying lithosphere's hydration and frictional strength (see Methods for model setup).

Our numerical experiments start with a homogeneous solid Earth of pyrolytic (BSE) composition$^{27}$ at 4.5 Ga with a low (0.01 wt. %) initial water content$^{11}$, mantle potential temperature of 1900 K and a core temperature of 5000 K. They end after 1.5 billion years, i.e. at
300

3.0 Ga. Any cataclysm relating to the formation of the mantle or the Moon is assumed to be
earlier than 4.5 Ga.

302

<figure><img src="image_9.png" /><figcaption>Figure 6: Representative models from Group I (A,B,D1-D4) and Group II (E,F,G1-G4). (A,E) Mobility of models over time with a moving average (thick line) over 11 mobility values (thin lines). (B,F) Mass of generated continental crust over time, scaled with the present-day</figcaption></figure>

304

306
continental crust mass(CCM). (C,D1-D4,G1-G4) Snapshots of Earth's cross sections showing the evolution of models with time where four quadrants represent different fields (clockwise from top left: composition, Nb/U, viscosity, $^{87}$Sr/$^{86}$Sr). Snapshot (C) shows initial model state corresponding to an age of 4500 Ma and model time of 0 Myr. Snapshots (D1,G1) show periods of active subduction, while snapshots (D2,D4,G2,G4) show periods of lulls in subduction activity where large mantle plumes (dark blue in the viscosity field) are trapped beneath the cold recycled material (yellow in the viscosity field). (E) shows a mobility burst at ca. 680 Myr of model time characterised by a resurfacing event (G3) and a spike in the production of continental crust (F).

All models fall into two different groups. The Group I consists of models with our preferred effective friction coefficient of the lithosphere for the early Earth of $0.1^{\text{ref.46}}$, while Group II consists of the models with higher lithospheric strength. The representative models from both groups are shown in Figure 6. Typical behavior of all models is intense mantle convection during the first hundred million years (Figure 6D1, G1), when multiple large hot mantle plumes form at the core-mantle boundary (CMB) and repeatedly approach the lithosphere. Plumes break the lithosphere and induce its subduction (Figure 6D1, G1). During this stage, there is extensive production of both oceanic and continental crust (Figure 6B, F). The following evolution depends on the lithospheric strength and is best described by the so-called mobility function which represents a ratio of the root-mean-square (rms) of the surface velocity averaged over the rms velocity of the entire computational domain $^{47}$. Mobility value above 1 means that lithosphere moves horizontally relative to underlying mantle with significant velocities and subduction is active$^{12}$.

In models of Group I, high mobility periods last several tens of million years and interchange with the low mobility periods (Figure 6A). We call this tectonic regime an oscillating mobile-lid regime. High mobility periods correspond to periods of extensive crustal production and lithospheric subduction and result in an accumulation of a large amount of cold recycled material at the CMB, which hinders the propagation of large mantle plumes (Figure 6D2). The plumes
336

that manage to rise through this overlying cold recycled material are small and are not able to
break the lithosphere or induce subduction, resulting in a reduced production of new crust (and
more recycling of existing crust) of the crust (Figure 6B, model time range 450-700 Myr). Then
it takes a few hundred million years to heat this cold material that sits above the CMB. Thus,
there is a period of quiescence before large plumes break through, approach the lithosphere and
induce new subduction zones, which produce and recycle new crust (Figure 6B, D3).

338

340

342

Group II models are characterized by low mobility throughout the evolution and experience
several mobility bursts that last several million years (Figure 6E). This type of tectonic regime is
usually called episodic-lid regime¹². During these episodes of high mobility, the overthickened
lithosphere is rapidly subducted resulting in resurfacing events (Figure 6G3). Most of the
continental crust is produced during these events. Note that despite showcasing different tectonic
regimes, the models from both groups may produce similar amounts of continental crust and
most of it during Hadean (Figure 6B,F), while models of Group II require for a larger input of
the surface water.

344

346

348

350

352

In all models, at every time step we identify regions of hot material (potential temperature above
1700°C) in the transition zone as the potential sources of komatiites (see Methods for details)
and select the model's cells that represent Hadean fraction of material with the ⁸⁷Sr/⁸⁶Sr ratio <
0.6997. The evolution of the trace element ratios in the Hadean fraction for a number of models
from both groups are presented in Figure 7a,e and Figure 8a,e. The common feature of all
models is that trace elements do not show any effect of continental crust production till at least
500 Myr of model time (i.e., 4.0 Ga) even though the continental crust production during the first
500 Myr is more than twice that of the last 1000 Myr in most models. This time lag indicates the

354

356

358
360

time required by the hydrated basaltic crust, which was involved in the production of continental crust to recycle and mix with the plume source material in the lower mantle.

362

The models of Group I that produce 40 to 70% of present-day continental crust mass (CCM) during the Hadean, show a good fit with the Nb/U and Ce/Pb data at the time close to the emplacement of Weltevreden komatiites at 3267 Ma (i.e., model time of 1230 Myr, Figure 5 and Figure 7a,d). These models also predict a realistic amount of Hadean component in plumes (Figure 7c,f) and a rather moderate supply of water into the lithosphere (0.5-0.9 present-day

364

ocean mass (OM)) over the period of 1500 Myr. In contrast, the models of Group II where we chose the water input parameters such that the models produce the similar amount of continental crust during the Hadean as the models of Group I, show a poor fit with Nb/U and Ce/Pb data (Figure 8a,d), a too low amount of Hadean component in plumes (Figure 8c,d) and predict a

366

368

water supply of 1.1-2.5 OM into the lithosphere over the period of 1500 Myr. The reason behind these different model predictions is the much lower recycling rate of the lithosphere in the models of Group II where the mobility, and therefore the subduction activity, are lower than in Group I models.

370

372

374

376
Group I: Models with high subduction activity

<figure><img src="image_10.png" /></figure>

**Figure 7:** A subset of results of models of Group I compared with observations. (a,d)- Evolution of trace element ratios Nb/U and Ce/Pb in the Hadean fraction (with $^{87}$Sr/$^{86}$Sr < 0.6997) of hot plume material in the transition zone from multiple models with different lithospheric hydration (see Methods for details). Each colored solid circle represents mean trace element ratio for the selected cells with error bars with 2 standard deviations of mean value. (b,e)- Mass of generated continental crust over time from the same models, scaled with the present-day continental crust mass. (c,f)-Amount of Hadean fraction in the hot plume material from the same models. Red solid diamonds (main Hadean group data) and blue diamonds (single inclusion data) show geochemical observations with error bars of 2 standard deviations of mean value.

Our models infer that in order to fit the geochemical observations, the following two conditions must have been fulfilled. First, the water supply into the lithosphere must have been high enough to allow for ca 40 to 70% of CCM production of continental crust during the Hadean in agreement with maximal estimates of 50 to 100% of CCM using mass balance approach
assuming whole mantle processing (Figure 5). Second, the tectonic regime in the Hadean and the Eo-Archean time must have been sufficiently mobile, i.e. with periods of extensive subduction.

Classical episodic regimes with rare partial resurfacing and sagduction or stagnant lid regime are not consistent with geochemical observations. We note that “subduction” in our models, defined as sinking into the mantle of the pieces of the entire upper thermal boundary layer (lid) including its surface, does not look like the present-day continuous and stable subduction that is

characteristic of Phanerozoic plate tectonics. It is much more variable and short-lived. Moreover, our models show that in Hadean-Archean Earth in contrast to the present-day Earth, large mantle plumes are of key importance to trigger such subduction and when plume activity is diminished, the subduction activity is also diminished.

404
Group II: Models with low subduction activity

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

**Figure 8: A subset of results of models of Group II compared with observations. The legend is the same as in Figure 7. Note poor fit of the geochemical data (red and blue diamonds) compared with the models of Group I shown in Figure 7.**

Altogether, the geochemical data and geodynamic models presented in this study indicate that large volumes of both oceanic and continental crust had already formed and recycled by the late-Hadean, thereby favoring the—still contentious—model presented by Richard Armstrong more than four decades ago³⁷ as well as several recent models³,⁶
418

# Methods

## Methods for geochemical investigation

420

### • Samples

422

424

426

428

430

432

434

436

438

440

The eight olivine-cumulates were sampled in the best-preserved parts of individual komatiite flows in the Saw Mill and Pioneer Complexes of the Weltevreden Formation of the Barberton greenstone Belt, South Africa, Fig. S 1). These cumulates consist of partially altered olivine crystals, and ground mass spinel and clinopyroxene grains enclosed in a matrix of altered glass (chlorite and amphibole). Three samples (1521, 1522 and 1523) from the Saw Mill Complex were previously described¹⁵. Sample 2216 (Gary's flow 1, 25.842597°S, 30.885031°E) is an olivine cumulate that contains large (up to 2 mm in diameter) partially serpentinized euhedral olivine grains with high Fo contents (92.8–93.1 mol.% Fo). The interstitial groundmass of the rock is made of acicular clinopyroxene (0.5-1 mm long), equant or skeletal spinel (up to 50 µm in diameter), and altered volcanic glass. Sample 2217 (Keena's flow 2, 25.841475°S, 30.885683°E) is an olivine cumulate consisting of best preserved out of all samples, but still partially serpentinized, euhedral olivine grains of different size (mainly 0.6–1 mm, with individual grains up to 1.5 mm in diameter; 93.2–93.8 mol.% Fo). The interstitial groundmass is made of acicular clinopyroxene (up to 1 mm long), equant and skeletal spinel crystals (up to 100 µm in diameter) and altered volcanic glass. Sample 2218 (Keena's flow 1, 25.841803°S, 30.886147°E) is an olivine cumulate consisting of partially serpentinized euhedral olivine grains of different size (mainly 0.5–2 mm, with exceptional grains up to 4 mm in length; 94.6–95.1 mol.% Fo). The interstitial groundmass is made of acicular clinopyroxene (up to 0.5 mm long), equant and skeletal spinel crystals (up to 100 µm in diameter) and altered volcanic glass. Sample 2222 (Gary's flow 2, 25.843578°S, 30.887742°E) is an olivine cumulate that contains large (up
442

444

446

448

to 2 mm in diameter) partially serpentinized euhedral olivine grains with high Fo contents (94.8–95.2 mol.% Fo). The interstitial groundmass of the rock is made of clinopyroxene (up to 1 mm long), equant spinel (up to 100 µm in diameter) and altered volcanic glass. Sample 1528 from the Pioneer Complex (25.834308°S, 30.949083°E) is an olivine cumulate consisting of partially serpentinized euhedral olivine grains or olivine pseudomorphs of different size (mainly 0.5–1.5 mm; 92.5–93.5 mol.% Fo, Figure 1, main text). The interstitial groundmass is made of acicular clinopyroxene (up to 1 mm long), equant and skeletal spinel crystals (up to 100 µm in diameter) and altered volcanic glass.

## Melt inclusions

450

452

454

456

458

460

462

464

Prior to polishing and mounting, melt inclusions hosted in olivine crystals were homogenized at high temperature and 1 atm. To homogenize partially crystallized melt inclusions, olivine crystals were placed in open platinum ampoules and heated in a vertical high-temperature furnace at Vernadsky Institute of Geochemistry (Moscow) following previously established protocol¹⁷. Target temperatures for experiments range from 1280°C to 1350°C in order to complete melting of crystalline phases in inclusions, and to minimize diffusional loss of hydrogen from inclusions through the host olivine. Upon heating completion, instantaneous quenching of the melt inclusion is achieved via immersing ampoules in water at room temperature. Because inclusions were heated at 1 atm pressure under reduced gas flow, all inclusions which were opened even partly lost their volatiles: e.g. Cl, S, H₂O as well as part of Na and most of Pb. Such inclusions, which can be identified by electron microprobe analysis, have Cl and S typically below detection limit, and were neglected in this study. Correction of the measured composition of glass in melt inclusions (see below) for post entrapment olivine crystallisation on the walls of inclusion's cavity and Fe-Mg exchange with host olivine (PEC) was performed using Petrolog 3 software software⁴⁸. Details on this approach are reported in
earlier studies¹⁵. PEC corrected inclusions compositions presented in Supplementary Tables 4 and 5.

### Electron microprobe analyses

The chemical composition of the olivine and melt inclusions were analysed using a JEOL JXA-8230 electron microprobe and a JEOL FEG JXA-iHP200F electron microprobe at ISTerre Micro-Analytical Platform (IMAP), Université Grenoble Alpes (France) using an analytical protocol published elsewhere⁴⁹. The Mg-Fe olivine-melt geothermometer was used for the estimation of komatiitic melt inclusion crystallization/entrapment temperature, with correction for the effects of measured water contents⁵⁰.

### Raman microscope

Water content of melt inclusions was determined using Raman microscope LABRAM Soleil, Horiba (ISTerre, Grenoble, France), equipped with 473 nm and 532 nm lasers and optical microscope (Nikon Eclipse LVDIA-N) using an analytical protocol previously published⁴⁹. Average and 2SE uncertainties are reported in the Supplementary Tables 3, 4 and 5.

### Laser Ablation Split Stream (LASS) ICP-MS analyses

LA-ICP-MS analyses were conducted at IMAP (Université Grenoble Alpes, France). Strontium (Sr) isotope signatures and trace element (TE) contents were measured directly in melt inclusions with a laser ablation (LA) system (RESOLution SE, Applied Spectra) coupled with multi-collector (Neptune XT, Thermo Fisher Scientific) for Sr-isotope and single-collector (8900, Agilent) for TE inductively-coupled plasma mass spectrometer (ICP-MS). LA single stream and LA split stream (LASS) protocols were used depending on the size of the melt inclusions: smaller inclusions were analyzed with our LASS Sr-TE analytical protocol while bigger ones had two laser analyses, one for Sr-isotopes and one for TE's. Uncertainties are reported at 95%
488 confidence interval unless stated otherwise. It is worth mentioning that Sr content in
Weltevreden melt inclusions average 36±6 ppm, corresponding to a mass of total Sr of ca. 6 pg
490 in most melt inclusions. Currently, such small amounts of Sr cannot be accurately measured via
the traditionally more precise TIMS approach owing to the lowest Sr blanks obtained through
492 such approach a few tenths of pg⁵¹.

Sr-isotope compositions were measured on a Neptune XT MC-ICP-MS (Thermo-Fisher
Scientific) following conventional cup configurations. Additionally, we also measured mass 90
494 in order to evaluate the Sr/Zr of Weltevreden melt inclusions and compare these ratios with the
trace element data obtained concurrently. Measurements were conducted in static mode and all
496 signals were collected in Faraday cups. Four 10¹³ Ω amplifiers were used for the measurement of
masses 84, 85, 86 and 87. All other masses (83, 83.5, 86.5, 88 and 90) were measured with 10¹¹
498 Ω amplifiers. Gain calibration was conducted overnight before each analytical run of melt
inclusions. Integration time was set to 1.049 seconds. The instrument was run in low-resolution
500 mode (mass resolution of ca. 400). Ni Jet-Sampler and Ni X-Skimmer cones were used. Sr-
isotope analyses were reduced offline using an in-house modified python code compatible with
Iolite-4⁵². All channels were first background subtracted. Since 83.5 and 86.5 signals were
502 always below detection limit (calculated as 3.29 × standard deviation of the background over the
entire analytical run⁵³), no REE²⁺ correction was conducted. Calcium-argides (Ca–Ar) species,
504 which can cause interference on masses 84, 86 and 88) are corrected using the background
subtracted signal of mass 83. A power-law based Sr-fractionation factor (β⁵¹) was then calculated
506 for each individual analysis based on the measured (88/86) value and the natural value
(8.37520938). The Rb fractionation factor (β⁶²) was manually adjusted to provide the best
508 accuracy for all three reference materials measured throughout the analytical run. This approach
allows to account for day-to-day variations in β⁶² and is analogous to the approach used for
510
correcting Yb interference in zircon Hf-isotope analysis. The amount of ⁸⁷Rb is then calculated using the measured ⁸⁵Rb and the natural ⁸⁷Rb/⁸⁵Rb value (0.38571). Then, the amount of ⁸⁷Sr is obtained via removing the contribution of ⁸⁷Rb on the measured 87 signals (87 net signal = ⁸⁷Rb + ⁸⁷Sr). Isotopic ratios are then calculated including a power-law based mass bias correction. For all melt inclusions labelled 15XX, adjusting βRb value on daily basis was sufficient to obtain accurate results (within 2SD) for all three reference glasses covering all compositional range of measured inclusions used in each analytical run (KL2-G, NIST-614-G and GOR128-G, see Fig. S 4), no standardization of the data was conducted here for ⁸⁷Sr/⁸⁶Sr values. Standardization was however required for ⁸⁷Rb/⁸⁶Sr, and Sr/Zr values. This is likely the result of minor elemental fractionation causing inaccuracy for corrected-only ratios. For melt inclusions labelled 22XX, standardization of both ⁸⁷Sr/⁸⁶Sr and ⁸⁷Rb/⁸⁶Sr was conducted with KL2-G as calibration glass. Readers can refer to the Sr-isotope DRS provided on the Iolite Github (https://github.com/iolite-LA-ICP-MS) for full view of the python code. Individual analysis uncertainties are calculated by the in-house built Iolite code. Hence, exports from Iolite include uncertainties of all background, interference corrections and standardization if conducted. Note that including REE²⁺ correction changes the final uncertainty to the sixth decimal only, i.e., 2 orders of magnitude lower than typical precision obtained with our LASS protocol for such Sr-depleted melt inclusions. Since several analytical runs were conducted for analyzing these melt inclusions, an evaluation of the inter-sequence (i.e. long-term) excess variance (i.e. long-term variance) was conducted. This approach follows the widely accepted workflow used in U–Pb isotope analyses⁵⁴ and consists of calculating the degree of homogeneity of the reference glasses over all sequence: if the MSWD value—or corresponding p(χ²) value—is within acceptable range of 5–95%, then no inter-sequence excess variance is propagated into individual analyses. If not, then uncertainties are expanded up to the level where all validation materials (here glasses) return MSDW value or
p(χ²) value indicating statistical homogeneity. No inter-sequence excess variance was required for the ⁸⁷Sr/⁸⁶Sr values: MSWD and p(χ²) values for the three reference glasses GOR128-G, KL2-G and NIST-614-G are: (i) 1.16 and 8.9%, (ii) 0.41 and 100%, (iii) 1.02 and 43% respectively (see data in Supplementary Table 1), while an inter-sequence excess variance of 5.5%, 6% and 6% was propagated into ⁸⁷Rb/⁸⁶Sr and Sr/Zr values respectively for all analyses (reference glasses and melt inclusions). Similar to the standardization mentioned above, the need of inter-sequence excess variance for these ratios likely reflects the contribution of elemental fractionation, which is unaccounted for in our data reduction protocol; an effect absent from ⁸⁷Sr/⁸⁶Sr calculation, hence explaining the better accuracy over several analytical runs for this latter ratio. Consequently, we consider the uncertainty reported in individual ⁸⁷Sr/⁸⁶Sr, ⁸⁷Rb/⁸⁶Sr and Sr/Zr values as external uncertainty. Precision and accuracy of Sr isotope analyses of melt inclusions were double checked by duplicate, triplicate or quREPlicate analyses of sufficiently large inclusions (Supplementary Fig. S 5 and Table 3c). The homogeneity of the replicates of 2, 3 or 4 analyses of the same melt inclusion presented in this diagram confirms that reported uncertainties of individual inclusions are not underestimated.

Lastly, we checked robustness of ⁸⁷Sr/⁸⁶Sr errors of accepted individual inclusions by comparison with the errors of reference materials using the following equation:

$$ \sigma_{inclusion} = \sigma_{RM} \times \sqrt{\frac{Sr_{RM} \times t_{RM}}{Sr_{inclusion} \times t_{inclusion}}} $$
Where $\sigma_{inclusion}$ and $\sigma_{RM}$ are 95% c.i. uncertainties of individual inclusions and reference material (glasses) respectively; $Sr_{inclusion}$ and $Sr_{RM}$ are Sr contents in inclusions and reference materials; $t_{inclusion}$ and $t_{RM}$ are $^{87}$Sr/$^{86}$Sr acquisition times for inclusions and reference materials (Supplementary Table 2).

Both predicted average $\sigma_{inclusion}$ calculated for all Weltevreden melt inclusions and for the 14 unradiogenic ones as well as their measured $\sigma_{inclusion}$ are 0.0010, thus virtually identical (Supplementary Table 2). We therefore conclude that reported errors of $^{87}$Sr/$^{86}$Sr ratios for individual melt inclusions are neither over nor underestimated and thus are valid for statistical treatment.

TE analyses were conducted with 20-µm spots/10-Hz pulse frequency and 38-µm spots/15-Hz pulse frequency for inclusions in single and split stream mode respectively and a laser fluence of 4 J.cm$^{-2}$. Carrier gas was He (ca. 0.3-0.9 L/min) with addition of N$_{2}$ (ca. 1 to 2.5 mL/min) which were mixed with Ar (ca. 0.6-1.1 L/min) before introduction into spectrometer. Oxide production rate, monitored with ThO$^{+}$/Th$^{+}$ was < 0.1%, and the doubly charged ratio monitored with $^{44}$Ca$^{2+}$/ $^{44}$Ca$^{+}$ was < 0.1%. The U/Th ratio ranged between 98% and 102%. Analyses were conducted in time-resolved acquisition mode (TRA) and included 1 second of pre-ablation to eliminate surface contamination, 30-40 s background measurement followed by 30-40 s sample ablation and signal measurement. Dwell time was 10-100 ms for different elements. All spectra were inspected in LADR software to define intervals for integration and exclude remaining surface contamination, if any. Concentrations were quantified from the measured ion yields normalized to Ca, previously measured on EPMA. Details on the varied sweep times and TE data reduction are reported in Supplementary Table 3a. Ce/Pb, Nb/U and Rb/Sr of reference glasses (GOR128-G$^{55}$ and NIST614-G$^{56}$) are shown in Fig. S 4 and reported in Supplementary Table 1. Reference-glass values agree with accepted values. Where required, an inter-sequence
excess variance was propagated into our analyses to account for minimal heterogeneity of
reference glasses between the varied analytical runs (up to 10%).

• Statistical analysis of the Sr-isotope dataset

Measured ⁸⁷Sr/⁸⁶Sr obtained on Weltevreden melt inclusions show a range of values that correspond to a MSWD of 2.3. Considering the number of individual datapoints—n=137—this MSWD value translates to a p(χ²) of 4.4×10⁻¹⁶, i.e., largely below the accepted threshold of 5% typically used in geoscience to indicate homogeneity of a sample population⁵⁷. This significant heterogeneity is also clearly observed for (⁸⁷Sr/⁸⁶Sr)initial values that show a MSWD of 2.2 translating to a p(χ²) of 9.8×10⁻¹⁴. Further, it is also clear from Fig. S 6 that the range of (⁸⁷Sr/⁸⁶Sr)initial values is not resulting from under/over correction of the (⁸⁷Sr/⁸⁶Sr)measured values: the range of ⁸⁷Rb/⁸⁶Sr is rather narrow and all analyses with elevated ⁸⁷Rb/⁸⁶Sr, i.e. > (⁸⁷Rb/⁸⁶Sr)KL2-G, were rejected from the dataset, preventing inaccurate over correction of (⁸⁷Sr/⁸⁶Sr)measured. Further, since all potential sources of uncertainty have been examined and propagated in isotopic ratios, this range of (⁸⁷Sr/⁸⁶Sr)initial values is interpreted to reflect Sr-isotope heterogeneity of the melts that were trapped in Weltevreden olivine crystals. The observed heterogeneity in elemental and isotopic compositions in Weltevreden melt inclusions is not a unique feature since a multitude of melt inclusion studies have reported similar observations¹⁶,¹⁸,⁵⁸ for different radiogenic isotope systems, regardless of the analytical protocol (laser ablation, ion probe or thermal ionization mass spectrometry).

Hence, we can discriminate statistically homogeneous populations in this dataset. One way to achieve this is to isolate individual analyses that are in the tails of this distribution: those with (⁸⁷Sr/⁸⁶Sr)initial significantly-i.e., considering individual uncertainties (±2σ)-outside the 95% confidence interval around the weighted mean of the dataset. Doing so yields three homogeneous populations of melt inclusions: (i) one with low (⁸⁷Sr/⁸⁶Sr)initial=0.69932±0.00024 (c.i. 95%)

584

586

588

590

592

594

596

598

600

602

604

606
comprised of 14 inclusions (MSWD = 0.88, p(χ²) = 0.57), (ii) one with intermediate ($^{87}$Sr/$^{86}$Sr)initial=0.700624±0.000086 made of 116 inclusions (MSWD = 1.1, p(χ²) = 0.27) and one with elevated ($^{87}$Sr/$^{86}$Sr)initial=0.70177±0.00024 made of 7 inclusions (MSWD = 0.48, p(χ²) = 0.82).

Applying of the t-test for independent populations with unequal variances⁵⁹ reveals that the means of ($^{87}$Sr/$^{86}$Sr)initial of all groups are different with the confidence over 99.9% (see Supplementary Table 6a,b).

## Model age calculation

Sr model age of olivine-hosted melt inclusions was based on (i) ($^{87}$Sr/$^{86}$Sr)initial and (ii) Bulk Silicate Earth (BSE) evolution. The latter is defined as follows: initial $^{87}$Sr/$^{86}$Sr = 0.698990 at 4.567 Ma age⁶⁰, initial BSE Rb/Sr=0.03²⁷ and $^{87}$Rb decay constant of 1.3972×10⁻¹¹ yr⁻¹. Even considering a 15% uncertainty in the BSE Rb/Sr value Sr-model ages remain far within uncertainty of our estimates. Calculation was conducted assuming no Rb left in the BSE source after melt extraction. This assumption is quite realistic since (i) degree of mantle melting is high in the hot early Earth and (ii) Rb is highly incompatible during melting. Samarium-Neodymium (Sm-Nd) and Lutetium-Hafnium (Lu-Hf) model ages presented in Fig. S 2 were calculated using Ryan Ickert's spreadsheet⁶¹.

## Contamination of Weltevreden komatiites

Studied inclusions in olivine yield significant intercorrelated variation of Cl, K, Rb, H₂O, Pb, Sr and Na, which is unrelated to olivine crystallization and extraction from the melt. This is well demonstrated by the strong correlation between ratios of elements excluded from olivine structure, e.g. Cl, K, Rb, H₂O, Pb, Sr, Na, and Ti, which must stay constant during olivine fractionation (Fig. S 3). Similar variations are common for submarine glasses of OIB and MORB
630

(Fig. S 3) and have been explained by minor (less than 1 wt%) contamination of seawater-derived ultra-saline brines<sup>25,26</sup>. The same explanation is also well suited for Weltevreden melt inclusions. The maximum amount of brine with a concentration of Cl of 10-30 wt% required to explain the addition of ca 0.20wt% of Cl to Weltevreden melt is about or less than 2wt%. This amount will not affect major element concentrations considering their analytical errors, but is enough to be observed in trace elements and H<sub>2</sub>O contents. Also note that the extent of contamination tends to increase with decreasing host olivine Fo content (Figure 2c, Fig. S 3).

632

634

636

The initial <sup>87</sup>Sr/<sup>86</sup>Sr isotope ratio of melt inclusions does not correlate with contamination proxies (see Figure 2h); thus, we conclude that the effect of such small amounts of assimilated brine is within the analytical errors of <sup>87</sup>Sr/<sup>86</sup>Sr determination.

638

640

**Methods for numerical investigation**

• Geodynamic models: essentials and setup

642

We study the thermochemical evolution of compressible mantle using the code StagYY<sup>41</sup>, which has been developed and widely used over several decades for global-scale modelling of Earth's evolution spanning its age<sup>11,43-45</sup>. The model includes pressure- and temperature-dependent viscosity, plastic yielding, time-dependent radiogenic heating, core cooling and phase changes.

644

Following a two-step mantle differentiation and utilizing water-dependent solidi functions<sup>43</sup>, the code forms basaltic and felsic melts and considers both intrusive (plutonic) and eruptive (volcanic) magmatism<sup>43</sup>. The mass ratio of erupted to intruded material follows a specified constant value of 30:70 (corresponding to 30% eruption efficiency), which has been previously shown to be important for producing Archean TTG rocks<sup>43,58,62</sup>. We employ temperature- and pressure- dependent solubility maps for different mantle minerals, which control the water in-gassing and out-gassing<sup>11</sup>.

646

648

650

652
We use a two-dimensional quadrant spherical annulus geometry⁴² with a Eulerian mesh, whose resolution varies radially, and is higher at the core-mantle boundary, 660 km phase transition and at the surface. The computational domain consists of 512 (lateral) times 128 (radial) cells and ~ 1.3 million Lagrangian tracers carrying various quantities (temperature, composition, water content, isotopes, trace elements, etc.) are advected through it. Free slip boundary conditions are used at the core mantle boundary and at the surface. Side boundary conditions are periodic.

A visco-plastic rheology is considered where the viscous deformation is accommodated by grain-size independent diffusion creep. The mantle is divided into 3 different layers *i*: upper mantle (1), lower mantle (2) and post-perovskite layer (3), with each layer having different values for activation energy $E_i$ and activation volume $V_i^{63,64}$. Following the Arrhenius law, the diffusion and dislocation creep (proxy) components of the temperature- and pressure-dependent viscosity are given as:

$$ \eta(T, P) = \eta_0 \Delta\eta_i \exp\left(\frac{E_i+PV(P)}{RT} - \frac{E_i}{RT_0}\right) \quad (1) $$

where $\eta_0$ is the reference viscosity (2.35 10¹⁹ Pa s) at zero pressure and reference temperature $T_0$ (1710 K), $\Delta\eta_i$ is the viscosity offset between layer *i* and the reference viscosity, *P* is the pressure, *R* is the gas constant and *T* is the absolute temperature. The activation volume decreases exponentially with increasing pressure in each layer *i* according to the relation:

$$ V(P) = V_i \exp\left(-\frac{P}{P_i}\right) \quad (2) $$

where $P_i$ is the pressure scale, which is different for each layer *i*. Plastic deformation in the lithosphere is accommodated by assuming plastic yielding⁴⁷,⁶⁵. The maximum stress that a material can sustain before deforming plastically is given by the yield stress $\sigma_Y$, which has both brittle and ductile components:
$$ \sigma_Y = \min(\sigma_{Y,\text{ductile}}, \sigma_{Y,\text{brittle}}). \qquad (3) $$

The ductile yield stress $\sigma_{Y,\text{ductile}}$ increases linearly with pressure as:

$$ \sigma_{Y,\text{ductile}} = \sigma_Y^0 + \sigma_Y P, \qquad (4) $$

where $\sigma_Y^0$ is the surface ductile yield stress and $\sigma_Y$ is the pressure gradient of the ductile yield stress. The brittle yield stress $\sigma_{Y,\text{brittle}}$ is calculated following the Byerlee approach as:

$$ \sigma_{Y,\text{brittle}} = c + \mu P, \qquad (5) $$

where c is the surface cohesion and $\mu$ is the friction coefficient. If the convective stresses exceed the yield stress, the viscosity is reduced to the yielding viscosity $\eta_Y = \sigma_Y/2\dot{\epsilon}$, where $\dot{\epsilon}$ is the 2nd invariant of the strain-rate tensor.

See Supplementary Table 7 for model parameters used in this study. The model uses a parameterisation based on mineral physics data$^{66,67}$, in which the mantle minerals are divided into olivine (ol), pyroxene-garnet (px-gt), TTG (tonalite-trondhjemite-granodiorite), and melt phase systems. Solid-solid phase transitions are assumed (see Supplementary Table 8) and the mantle is initialised with a pyrolytic composition: 80% harzburgite and 20% basalt, being a mixture of 60% olivine and 40% pyroxene-garnet. The adiabatic temperature, density, thermal conductivity, thermal expansivity and heat capacity are pressure-dependent following a third-order Birch–Murnaghan equation of state. Further details of the model can be found in$^{11}$. Here, we describe only the model modifications introduced in this study.

## Geodynamic models: modifications

1. Rheology: We improved our rheological model by adding a proxy for the pressure-, temperature- and stress-dependent dislocation creep rheology assuming stress corresponding to a geological strain rate of $10^{-15}$ 1/s and using dry olivine rheology parameters$^{68}$ (Supplementary Table 7). We also considered frictional plasticity with the
coefficient of friction of 0.1 corresponding to the properties of subduction channels, which lack continental sediments⁴⁶. Similar effective friction coefficient was suggested based on the model of thermal cracking of oceanic lithosphere on early Earth⁶⁹.

1. **Water solubility maps:** We modified water solubility map for px-garnet phases to better fit the map calculated by Perple_X (see Fig S8).

2. **Water in-gassing:** In this study, we assume that at every time-step, the mantle minerals only in the top 5 km of the computational domain (as opposed to 10 km done previously) are partially saturated with water in accordance with their pressure- and temperature-dependent solubility maps with water saturation fraction $\alpha_w$ parameter ($0 < \alpha_w < 1$) controlling the input of the surface water into the lithosphere.

3. **Density:** We corrected density of felsic material (TTG) by considering phase transformations of feldspars and quartz at high pressures (see Fig. S 9). We have also slightly modified density changes of olivine and pyroxene-garnet compositions at phase boundaries to decrease density difference between pyrolite and pyroxene-garnet in lowermost mantle to 1.7 %, simultaneously fitting density-depth distribution in PREM model along 1600K adiabatic geotherm. We introduced the effect of water on density using a simplified relation

$$ \rho_{\omega} = \rho_{d}(1 - aC_{H2O}), \qquad (6) $$

where $\rho_{\omega}$, $\rho_{d}$ are densities of water-containing and dry material respectively, $C_{H2O}$ is the concentration of water in weight fraction and $a$ is a constant, typically between 1 and 2 estimated using Gibbs energy minimization code⁷⁰ (see Fig. S 9) that we consider as 2.0 in this study.

4. **Conditions for felsic crust (TTG) formation from hydrated basalt:** In the present study, we assume that TTG can only be produced from basalt if the latter's water content is above 0.5
wt. (%) (instead of 0.05 wt. % used in earlier studies) and that maximum degree of melting of basalt to produce TTG is 30% instead of 10% used in previous studies¹¹,⁴³.

5. **Isotopic systems and trace elements:** We have introduced the evolution of Rb-Sr isotope system as well as Nb, U, Ce and Pb trace elements. Constant parameters and solid-melt partition coefficients used in the study are presented in Supplementary Tables 9, 10 respectively.

## Geodynamic models: definition of sampling for Figures 7,8 (main text)

For the models presented in Figures 7 and 8 (main text), we show the evolution of the mean Nb/U and Ce/Pb ratios with time for a selection of the model's cells fulfilling the following conditions: (i) Average cell temperature has to be more than 100K higher than average temperature of all cells located at the same depth. This condition allows us to identify mantle plume material. In practice, potential temperature of all selected cells in most of the models appear to be higher than 1700°C. (ii) Cells have to be located within the depth range of 400-600km. This condition is required to avoid cells from the top 300 km, from which the melt is extracted and trace elements are partitioned. (iii) Average initial Sr isotope ratio (⁸⁷Sr/⁸⁶Sr) of each cell has to be less than 0.6997 to identify the Hadean fraction. Following this procedure, we obtain mean values of trace elements ratios in the Hadean fraction of the potential komatiite source and standard deviations of mean values at every output model time. To avoid unrepresentative and too uncertain values, we further filter out values at times with too few Hadean cells (<=20), with too large value of 2 standard deviation of mean Nb/U ratio (>=4) and with too large 2 standard deviation of mean Ce/Pb ratio (>=2). As a result, we obtain mean values of trace element ratios and their standard deviations (multiplied by 2) for each model at every model output time, which we show in Figure 7& 8 (main text).
746

748

750

752

754

756

Due to inherent randomness of convection processes arising from the initial thermal perturbations and initial tracer positions, the model results vary even with the exact same parameters. Group I presented in Figure 7 comprises of two sets of models (3 in each set) with the same frictional strength ($f=0.1$) but with different water saturation fraction ($alpha=0.25$ or $0.20$). The higher the value of alpha, the higher is the input of water into the lithosphere. Group II presented in Figure 8 also comprises of two sets of models (2 in each set) with the same frictional strength ($f=0.2$) but with different water saturation fraction ($alpha=0.40$ or $0.30$). Values of alpha were chosen such that all models of both groups produce similar amount of continental crust during the Hadean.

## Geodynamic models: parameter sensitivity

758

760

762

764

766

In the main text we show the effect of lithospheric hydration and strength on the evolution of the trace elements ratios of the Hadean fraction of potential komatiite source (Figure 7). Here, we additionally show the influence of the variation of the parameter “*a*” (eq. 6) responsible for decreasing the density of rocks due to water-containing minerals. Panels a-h in Figure S10 show our models that fit well with the geochemical observations computed with our preferred value of “*a*” of 2. Panels i-p in Figure S10 show modelling results obtained with the value of “*a*” parameter of 1.5, which is more suitable for rocks at relatively high pressures (see Figure S9b) with all the other parameters same as in the models with $a=2$. As we see from Figure S10, the variation of “*a*” parameter does not change significantly the modeling results.
**Acknowledgements:**

We thank A. H. Wilson and G.R. Byerly for their help in fieldwork in Barberton Greenstone

Belt, South Africa in 2015, and Axel Hofmann and M. Elburg in 2022. M.V. Portnyagin for
providing reference glasses for calibration of Raman H₂O analysis. F. Brunet, B. Lanson, P.

Roux and the direction of University Grenoble Alpes provided invaluable help in establishing
and maintaining the ISTerre Micro-Analytical Platform (IMAP).

**Funding:**

This project is supported by a grant of the European Research Council (ERC) under the

European Union's Horizon 2020 research and innovation program (Synergy Grant MEET,
grant agreement no.856555) to A.V.S, S.V.S and J.W.V. ISTerre is part of the Laboratoires

d'Excellence (LabEx) OSUG@2020 (ANR10 LABX56). Participation of E.V.A. and A.N.K and
annealing of melt inclusions was supported by grant of Russian Research Foundation no 22-77-

00081 to E.V.A. The computing time for geodynamic modelling (project bbk00014) was
provided by the administration of the high-performance computer cluster "Lise" at the NHR

Centre NHR@ZIB jointly supported by the Federal Ministry of Education and Research and the
state government.

**Author contributions:**

A.V. participated in samples collection, developed the protocol for Rb-Sr isotope analyses and

LA-ICP-MS Rb-Sr/TE split stream analytical mode, conducted analyses, processed Rb-Sr
isotopes data, interpreted geochemical data and wrote most of the manuscript text; A.V.C.

participated in samples collection, found and prepared all melt inclusions and host minerals for
analytical work, conducted EPMA and Raman analyses, processed trace elements data and
participated in geochemical data interpretation; A.V.S. conceived the idea and designed the
project, assisted on the collection of samples and played leading role in interpreting geochemical
data; C.J. developed the geodynamic models, wrote and modified code and performed
computing; S.V.S. conceived the idea and designed the geodynamic model, played leading role
in interpreting and visualising results of geodynamic models and wrote most part of geodynamic
modeling section; E.V.A. performed inclusions homogenization, and assisted in the collection of
samples; A.N.K performed inclusions homogenization; V.G.B. developed analytical protocols for
EPMA and LA-ICP-MS analyses of melt inclusions and host olivines; N.T.A. governed in the
collection of samples and shared his experience in the study of komatiites; L.V.D. directed the
tuning of RESolution SE laser and Agilent 8900 ICP-MS and participated in preparation of LA-
ICP-MS analytical protocols; J.W.V. participated in the design of the project. All co-authors
contributed to discussing and interpreting data and writing the manuscript.

Competing interests

802 Authors declare that they have no competing interests

804
References

806 1 O'Neil, J., Carlson, R. W., Francis, D. & Stevenson, R. K. Neodymium-142 Evidence for Hadean Mafic Crust. *Science* **321**, 1828-1831 (2008).

808 https://doi.org:https://doi.org/10.1126/science.1161925

810 2 Bowring, S. A. & Williams, I. S. Priscoan (4.00–4.03 Ga) orthogneisses from northwestern Canada. *Contrib. Mineral. Petrol.* **134**, 3-16 (1999).

812 3 Korenaga, J. Crustal evolution and mantle dynamics through Earth history. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **376**, 20170408 (2018). https://doi.org:https://doi.org/10.1098/rsta.2017.0408

814 4 Campbell, I. H. Constraints on continental growth models from Nb/U ratios in the 3.5 Ga Barberton and other Archaean basalt-komatiite suites. *Am. J. Sci.* **303**, 319-351 (2003).

816 5 Dhuime, B., Hawkesworth, C. J., Cawood, P. A. & Storey, C. D. A change in the geodynamics of continental growth 3 billion years ago. *Science* **335**, 1334-1336 (2012).

818 6 Rosas, J. C. & Korenaga, J. Rapid crustal growth and efficient crustal recycling in the early Earth: Implications for Hadean and Archean geodynamics. *EPSL* **494**, 42-49 (2018). https://doi.org:https://doi.org/10.1016/j.epsl.2018.04.051

820 7 McCoy-West, A. J. et al. Extensive crustal extraction in Earth's early history inferred from molybdenum isotopes. *Nat. Geosci.* **12**, 946-951 (2019).

824 8 Dhuime, B., Hawkesworth, C. J., Cawood, P. A. & Storey, C. D. A change in the geodynamics of continental growth 3 billion years ago. *Science* **335**, 1334-1336 (2012).

826 9 Hofmann, A., Class, C. & Goldstein, S. Size and composition of the MORB+ OIB mantle reservoir. *Geochem. Geophys. Geosyst.* **23** (2022).

828 10 Whitehouse, M. Granulite fades Nd-isotopic homogenization in the Lewisian complex of northwest Scotland. *Nature* **331**, 705-707 (1988).

830 11 Jain, C., Rozel, A. B., van Hunen, J., Chin, E. J. & Manjón-Cabeza Córdoba, A. Building archean cratonic roots. *Frontiers in Earth Science* **10**, 966397 (2022).

832 12 Vezinet, A., Pearson, D. G. & Thomassot, E. Effects of contamination on whole-rock isochrons in ancient rocks: a numerical modelling approach. *Lithos* **386-387**, 106040 (2021). https://doi.org:https://doi.org/10.1016/j.lithos.2021.106040

834 13 Hofmann, A., Class, C. & Goldstein, S. Size and composition of the MORB+ OIB mantle reservoir. *Geochem. Geophys. Geosyst.* **23** (2022).

838 14 Jain, C., Rozel, A. B., van Hunen, J., Chin, E. J. & Manjón-Cabeza Córdoba, A. Building archean cratonic roots. *Frontiers in Earth Science* **10**, 966397 (2022).

840 15 Lourenço, D. L., Rozel, A. B., Ballmer, M. D. & Tackley, P. J. Plutonic-squishy lid: A new global tectonic regime generated by intrusive magmatism on Earth-like planets. *Geochem. Geophys. Geosyst.* **21**, e2019GC008756 (2020).

842 16 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

844 17 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

846 18 Reimink, J., Davies, J., Moyen, J. & Pearson, D. A whole-lithosphere view of continental growth. *Geochemical Perspectives Letters* **26**, 45-49 (2023).

848 19 Sobolev, A. V. et al. Deep hydrous mantle reservoir provides evidence for crustal recycling before 3.3 billion years ago. *Nature* **571**, 555-559 (2019).

850 20 Sobolev, A. V. et al. Deep hydrous mantle reservoir provides evidence for crustal recycling before 3.3 billion years ago. *Nature* **571**, 555-559 (2019).

852 21 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

854 22 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

856 23 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

858 24 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

860 25 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

862 26 Reimink, J., Davies, J., Moyen, J. & Pearson, D. A whole-lithosphere view of continental growth. *Geochemical Perspectives Letters* **26**, 45-49 (2023).

864 27 Hofmann, A. W., Jochum, K. P., Seufert, M. & White, W. M. Nb and Pb in oceanic basalts: new constraints on mantle evolution. *EPSL* **79**, 33-45 (1986).

868 28 Sobolev, A. V. et al. Deep hydrous mantle reservoir provides evidence for crustal recycling before 3.3 billion years ago. *Nature* **571**, 555-559 (2019).

870 29 Sobolev, A. V. et al. Deep hydrous mantle reservoir provides evidence for crustal recycling before 3.3 billion years ago. *Nature* **571**, 555-559 (2019).
852 16 Sobolev, A. V., Hofmann, A. W., Jochum, K. P., Kuzmin, D. V. & Stoll, B. A young source for the Hawaiian plume. *Nature* **476**, 434-437 (2011).
https://doi.org:https://doi.org/10.1038/nature10321

854 17 Sobolev, A. V. et al. Komatiites reveal a hydrous Archaean deep-mantle reservoir. *Nature* **531**, 628-632 (2016). https://doi.org:https://doi.org/10.1038/nature17152

856 18 Saal, A., Hart, S., Shimizu, N., Hauri, E. & Layne, G. Pb isotopic variability in melt inclusions from oceanic island basalts, Polynesia. *Science* **282**, 1481-1484 (1998).
https://doi.org:https://doi.org/10.1126/science.282.5393.1481

858 19 Jackson, M. G. & Hart, S. R. Strontium isotopes in melt inclusions from Samoan basalts: Implications for heterogeneity in the Samoan plume. *EPSL* **245**, 260-277 (2006).
https://doi.org:https://doi.org/10.1016/j.epsl.2006.02.040

862 20 McDonough, W. & Ireland, T. Intraplate origin of komatiites inferred from trace elements in glass inclusions. *Nature* **365**, 432-434 (1993).
https://doi.org:https://doi.org/10.1038/365432a0

864 21 Connolly, B. D. et al. Highly siderophile element systematics of the 3.3 Ga Weltevreden komatiites, South Africa: implications for early Earth history. *EPSL* **311**, 253-263 (2011).
https://doi.org:https://doi.org/10.1016/j.epsl.2011.09.039

868 22 Hofmann, A., Anhaeusser, C. & Li, X. Layered ultramafic complexes of the Barberton Greenstone Belt-age constraints and tectonic implications. *S. Afr. J. Geol.* **124**, 7-16 (2021). https://doi.org:https://doi.org/10.25131/sajg.124.0002

872 23 Puchtel, I. et al. Insights into early Earth from Barberton komatiites: evidence from lithophile isotope and trace element systematics. *Geochim. Cosmochim. Acta* **108**, 63-90 (2013). https://doi.org:https://doi.org/10.1016/j.gca.2013.01.016

874 24 Roerdink, D. L., Ronen, Y., Strauss, H. & Mason, P. R. Emergence of felsic crust and subaerial weathering recorded in Palaeoarchaean barite. *Nat. Geosci.* **15**, 227-232 (2022).
https://doi.org:https://doi.org/10.1038/s41561-022-00902-9

878 25 Coombs, M. L., Sisson, T. W. & Kimura, J.-I. Ultra-high chlorine in submarine Kılaua glasses: evidence for direct assimilation of brine by magma. *EPSL* **217**, 297-313 (2004).
https://doi.org:https://doi.org/10.1016/S0012-821X(03)00631-9

880 26 Kendrick, M. A., Arculus, R., Burnard, P. & Honda, M. Quantifying brine assimilation by submarine magmas: Examples from the Galápagos Spreading Centre and Lau Basin. *Geochim. Cosmochim. Acta* **123**, 150-165 (2013).
https://doi.org:https://doi.org/10.1016/j.gca.2013.09.012

884 27 McDonough, W. F. & Sun, S.-s. The composition of the Earth. *Chem. Geol.* **120**, 223-253 (1995). https://doi.org:https://doi.org/10.1016/0009-2541(94)00140-4

886 28 Hofmann, A. W. Chemical differentiation of the Earth: the relationship between mantle, continental crust, and oceanic crust. *EPSL* **90**, 297-314 (1988).
https://doi.org:https://doi.org/10.1016/0012-821X(88)90132-X

888 29 Lyubetskaya, T. & Korenaga, J. Chemical composition of Earth's primitive mantle and its variance: 1. Method and results. *Journal of Geophysical Research: Solid Earth* **112** (2007). https://doi.org:https://doi.org/10.1029/2005JB004223

892 30 Puchtel, I. S., Blichert-Toft, J., Horan, M. F., Touboul, M. & Walker, R. J. The komatiite testimony to ancient mantle heterogeneity. *Chem. Geol.* **594**, 1207-1276 (2022).
https://doi.org:https://doi.org/10.1016/j.chemgeo.2022.120776

894 31 Corgne, A., Liebske, C., Wood, B. J., Rubie, D. C. & Frost, D. J. Silicate perovskite-melt partitioning of trace elements and geochemical signature of a deep perovskitic reservoir.
Geochim. Cosmochim. Acta 69, 485-496 (2005).

898 https://doi.org:https://doi.org/10.1016/j.gca.2004.06.041

32 Tilhac, R., Begg, G. C., O'Reilly, S. Y. & Griffin, W. L. A global review of Hf-Nd isotopes: New perspectives on the chicken-and-egg problem of ancient mantle signatures. *Chem. Geol.* 609, 121039 (2022).

900 https://doi.org:https://doi.org/10.1016/j.chemgeo.2022.121039

33 Wedepohl, K. H. The composition of the continental crust. *Geochim. Cosmochim. Acta* 59, 1217-1232 (1995).

34 Rudnick, R. L. & Fountain, D. M. Nature and composition of the continental crust: a lower crustal perspective. *Rev. Geophys.* 33, 267-309 (1995).
https://doi.org:https://doi.org/10.1029/95RG01302

35 Taylor, S. R. & McLennan, S. M. *The continental crust: its composition and evolution*. (Blackwell Scientific Publications, 1985).

910 36 McLennan, S. M. Relationships between the trace element composition of sedimentary rocks and upper continental crust. *Geochem. Geophys. Geosyst.* 2 (2001).
https://doi.org:https://doi.org/10.1029/2000GC000109

37 Armstrong, R. L. Radiogenic isotopes: the case for crustal recycling on a near-steady-state no-continental-growth Earth. *Phil. Trans. R. Soc. London, Ser. A* 301, 443-472 (1981).
https://doi.org:https://doi.org/10.1098/rsta.1981.0122

916 38 Harrison, T. M. *Hadean Earth*. (Springer Cham, 2020).

39 Bauer, A. *et al*. Hafnium isotopes in zircons document the gradual onset of mobile-lid tectonics. *Geochemical Perspectives Letters* 14, 1-6 (2020).
https://doi.org:https://doi.org/10.7185/geochemlet.2015

918 https://doi.org:https://doi.org/10.1073/pnas.2120241119

40 Tusch, J. *et al*. Long-term preservation of Hadean protocrust in Earth's mantle. *Proceedings of the National Academy of Sciences* 119 (2022).

922 https://doi.org:https://doi.org/10.1016/j.pepi.2008.07.037

41 Tackley, P. J. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid. *Phys. Earth Planet. Inter.* 171, 7-18 (2008).
https://doi.org:https://doi.org/10.1016/j.pepi.2008.08.005

924
926
42 Hernlund, J. W. & Tackley, P. J. Modeling mantle convection in the spherical annulus. *Phys. Earth Planet. Inter.* 171, 48-54 (2008).
https://doi.org:https://doi.org/10.1016/j.pepi.2008.07.037

43 Jain, C., Rozel, A. B., Tackley, P. J., Sanan, P. & Gerya, T. V. Growing primordial continental crust self-consistently in global mantle convection models. *Gondwana Res.* 73, 96-122 (2019).
https://doi.org:https://doi.org/10.1016/j.gr.2019.03.015

928
930
43
Xie, S. & Tackley, P. J. Evolution of helium and argon isotopes in a convecting mantle. *Phys. Earth Planet. Inter.* 146, 417-439 (2004).
https://doi.org:https://doi.org/10.1016/j.pepi.2004.04.003

45 Xie, S. & Tackley, P. J. Evolution of U-Pb and Sm-Nd systems in numerical models of mantle convection and plate tectonics. *Journal of Geophysical Research: Solid Earth* 109 (2004).
https://doi.org:https://doi.org/10.1029/2004JB003176

934
936
45
Sobolev, S. V. & Brown, M. Surface erosion events controlled the evolution of plate tectonics on Earth. *Nature* 570, 52-57 (2019).
https://doi.org:https://doi.org/10.1038/s41586-019-1258-4

46
938
47
Tackley, P. J. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations. *Geochem. Geophys. Geosyst.* 1 (2000).
https://doi.org:https://doi.org/10.1029/2000GC000036

940 https://doi.org:https://doi.org/10.1038/s41586-019-1258-4

942
47
Tackley, P. J. Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations. *Geochem. Geophys. Geosyst.* 1 (2000).
https://doi.org:https://doi.org/10.1029/2000GC000036
944 48 Danyushevsky, L. V. & Plechov, P. Petrolog3: Integrated software for modeling crystallization processes. *Geochem. Geophys. Geosyst.* **12** (2011). https://doi.org:https://doi.org/10.1029/2011GC003516

946 https://doi.org:https://doi.org/10.1029/2011GC003516

948 49 Kazzy, C. et al. Study of olivine-hosted melt and spinel inclusions from the Song Da ultramafic volcanic suite, northern Vietnam: Compositions, crystallization temperatures, and origin of the low-Ti komatiite-like and high-Ti primary melts. *Chem. Geol.*, 122219 (2024). https://doi.org:https://doi.org/10.1016/j.chemgeo.2024.122219

950 50 Falloon, T. J. & Danyushevsky, L. V. Melting of refractory mantle at 1·5, 2 and 2·5 GPa under anhydrous and H₂O-undersaturated conditions: implications for the petrogenesis of high-Ca boninites and the influence of subduction components on mantle melting. *J. Petrol.* **41**, 257-283 (2000). https://doi.org:https://doi.org/10.1093/petrology/41.2.257

954 https://doi.org:https://doi.org/10.1093/petrology/41.2.257

956 51 Koornneef, J. M. et al. TIMS analysis of Sr and Nd isotopes in melt inclusions from Italian potassium-rich lavas using prototype 1013 Ω amplifiers. *Chem. Geol.* **397**, 14-23 (2015). https://doi.org:https://doi.org/10.1016/j.chemgeo.2015.01.005

958 52 Paton, C., Hellstrom, J., Paul, B., Woodhead, J. & Hergt, J. Iolite: Freeware for the visualisation and processing of mass spectrometric data. *J. Anal. At. Spectrom.* **26**, 2508-2518 (2011). https://doi.org:https://doi.org/10.1039/C1JA10172B

960 53 Currie, L. A. Limits for qualitative detection and quantitative determination. Application to radiochemistry. *Analytical chemistry* **40**, 586-593 (1968).

964 54 Horstwood, M. S. et al. Community-Derived Standards for LA-ICP-MS U-(Th-)Pb Geochronology–Uncertainty Propagation, Age Interpretation and Data Reporting. *Geostand. Geoanal. Res.* **40**, 311-332 (2016). https://doi.org:https://doi.org/10.1111/j.1751-908X.2016.00379.x

966 https://doi.org:https://doi.org/10.1111/j.1751-908X.2016.00379.x

968 55 Jochum, K. P. et al. MPI-DING reference glasses for in situ microanalysis: New reference values for element concentrations and isotope ratios. *Geochem. Geophys. Geosyst.* **7** (2006). https://doi.org:https://doi.org/10.1029/2005GC001060

970 56 Jochum, K. P. et al. Determination of reference values for NIST SRM 610–617 glasses following ISO guidelines. *Geostand. Geoanal. Res.* **35**, 397-429 (2011). https://doi.org:https://doi.org/10.1111/j.1751-908X.2011.00120.x

972 57 Wendt, I. & Carl, C. The statistical distribution of the mean squared weighted deviation. *Chemical Geology: Isotope Geoscience Section* **86**, 275-285 (1991). https://doi.org:https://doi.org/10.1016/0168-9622(91)90010-T

976 https://doi.org:https://doi.org/10.1016/0168-9622(91)90010-T

978 58 Reinhard, A. et al. Sr and Nd isotopic compositions of individual olivine-hosted melt inclusions from Hawai'i and Samoa: Implications for the origin of isotopic heterogeneity in melt inclusions from OIB lavas. *Chem. Geol.* **495**, 36-49 (2018). https://doi.org:https://doi.org/10.1016/j.chemgeo.2018.07.034

980 https://doi.org:https://doi.org/10.1016/j.chemgeo.2018.07.034

982 59 Welch, B. L. The generalization of ‘STUDENT’S’problem when several different population variances are involved. *Biometrika* **34**, 28-35 (1947). https://doi.org:https://doi.org/10.1093/biomet/34.1.2.28

984 60 Wasserburg, G., Papanastassiou, D. & Sanz, H. Initial strontium for a chondrite and the determination of a metamorphism or formation interval. *EPSL* **7**, 33-43 (1969). https://doi.org:https://doi.org/10.1016/0012-821X(69)90008-9

986 https://doi.org:https://doi.org/10.1016/0012-821X(69)90008-9

988 61 Ickert, R. B. Algorithms for estimating uncertainties in initial radiogenic isotope ratios and model ages. *Chem. Geol.* **340**, 131-138 (2013). https://doi.org:https://doi.org/10.1016/j.chemgeo.2013.01.001

988 https://doi.org:https://doi.org/10.1016/j.chemgeo.2013.01.001
990 62 Rozel, A., Golabek, G. J., Jain, C., Tackley, P. J. & Gerya, T. Continental crust formation on early Earth controlled by intrusive magmatism. *Nature* **545**, 332-335 (2017). https://doi.org/https://doi.org/10.1038/nature22042

992 https://doi.org/https://doi.org/10.1038/nature22042

994 63 Karato, S.-i. & Wu, P. Rheology of the upper mantle: A synthesis. *Science* **260**, 771-778 (1993). https://doi.org/https://doi.org/10.1126/science.260.5109.771

996 64 Yamazaki, D. & Karato, S.-I. Some mineral physics constraints on the rheology and geothermal structure of Earth's lower mantle. *Am. Mineral.* **86**, 385-391 (2001). https://doi.org/https://doi.org/10.2138/am-2001-0401

998 65 Moresi, L. & Solomatov, V. Mantle convection with a brittle lithosphere: thoughts on the global tectonic styles of the Earth and Venus. *Geophys. J. Int.* **133**, 669-682 (1998). https://doi.org/https://doi.org/10.1046/j.1365-246X.1998.00521.x

1000 66 Irifune, T. & Ringwood, A. Phase transformations in subducted oceanic crust and buoyancy relationships at depths of 600–800 km in the mantle. *EPSL* **117**, 101-110 (1993). https://doi.org/https://doi.org/10.1016/0012-821X(93)90120-X

1002 67 Ono, S., Ito, E. & Katsura, T. Mineralogy of subducted basaltic crust (MORB) from 25 to 37 GPa, and chemical heterogeneity of the lower mantle. *EPSL* **190**, 57-63 (2001). https://doi.org/https://doi.org/10.1016/S0012-821X(01)00375-2

1006 68 Hirth, G. & Kohlstedf, D. in *Inside the Subduction Factory* Vol. 138 *Geophysical monograph-American Geophysical Union* 83-106 (2003).

1008 69 Korenaga, J. Thermal cracking and the deep hydration of oceanic lithosphere: A key to the generation of plate tectonics? *Journal of Geophysical Research: Solid Earth* **112** (2007). https://doi.org/https://doi.org/10.1029/2006JB004502

1010 69 Sobolev, S. V. & Babeyko, A. Y. Modeling of mineralogical composition, density and elastic wave velocities in anhydrous magmatic rocks. *Surv. Geophys.* **15**, 515-544 (1994). https://doi.org/https://doi.org/10.1007/BF00690173
Supplementary Files

This is a list of supplementary files associated with this preprint. Click to download.

• 20231120122ArticleandMethodsv0.pdf

• 20231120122ExplanatoryFileforSOMTablesv0.pdf

• ExplanatoryFileforSOMTablesv0.docx

• SOMTables.xlsx

• 20231120122TablesSOMv0.xlsx