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

Adrien Vezinet, Aleksandr Chugunov, Alexander Sobolev, Charitra Jain, Stephan Sobolev, Valentina Batanova, Evgeny V. Asafov, Alina Koshlyakova, Nicholas Arndt, Leonid Danyushevsky, et al.

► To cite this version:

Adrien Vezinet, Aleksandr Chugunov, Alexander Sobolev, Charitra Jain, Stephan Sobolev, et al. Growth of continental crust and lithosphere subduction in the Hadean revealed by geochemistry and geodynamics. Nature Communications, 2025, 16 (1), pp.3850. 10.1038/s41467-025-59024-6. hal-05050084

HAL Id: hal-05050084

https://hal.science/hal-05050084v1

Submitted on 22 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Received: 31 October 2024

Accepted: 9 April 2025

Published online: 25 April 2025

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

Check for updates

Adrien Vezinet Ⓓ<fnref n="1" />,⁷ Ⓒ, Aleksandr V. Chugunov Ⓓ<fnref n="1" />,⁷ Ⓒ,
Alexander V. Sobolev Ⓓ<fnref n="1" />,⁷ Ⓒ, Charitra Jain Ⓓ², Stephan V. Sobolev²,³ Ⓒ,
Valentina G. Batanova Ⓓ<fnref n="1" />, Evgeny V. Asafov⁴, Alina N. Koshlyakova Ⓓ⁴,
Nicholas T. Arndt<fnref n="1" />, Leonid V. Danyushevsky⁵ & John W. Valley Ⓓ⁶

The rates of continental crust growth and recycling on early Earth remain unclear due to the lack of information resulting from the extensive alteration of ancient rocks. Melt inclusions trapped and shielded from alteration in Archean high-Mg olivine crystals offer a solution to this problem. We report an unprecedented unradiogenic Sr mantle source component (⁸⁷Sr/⁸⁶Sr = 0.69932 ± 0.00024, 95% confidence interval) of melts included in olivine from 3.27 Ga komatiitic lava flows in the Barberton Greenstone Belt, South Africa. This component indicates a model age of 4.31 ± 0.19 Ga and significant chemical fractionation (Nb/U = 36.9 ± 1.5, Ce/Pb=16.7 ± 1.1), suggesting up to 80% ± 16% of the present-day continental crust's mass was extracted by the late Hadean from the whole mantle. Geodynamic models support this finding, explaining geochemical data by producing 40% to 70% of the present-day continental crust mass during the Hadean in a variable tectonic regime with tens of millions of years-long periods of massive impulsive subduction induced by mantle plumes.

The evolution of the early Earth and the 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 more than 100% of the volume of today's continental crust³.

Several approaches have been used to provide estimates of con-
tinental 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<sup>3,4,6</sup>. This approach, however, suffers from significant
problems, such as known mantle heterogeneity in isotope ratios and
systematic differences in these ratios for the mantle sources of mid-
ocean ridge basalts (MORB) and oceanic island basalts (OIB)<sup>8</sup>. 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<sup>4</sup>. This uncertainty is also
valid for the short-lived <sup>146</sup>Sm-<sup>142</sup>Nd decay system<sup>9</sup>. Also, most ancient
rocks are altered, which has been shown to affect Sm-Nd and Lu-Hf
isotopic systems<sup>10,11</sup>. Further, whatever the volume of continental crust

<footnote n="1">¹Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTerre, Grenoble, France. ²GFZ Helmholtz Centre for Geosciences, Geodynamic Modelling Section, Potsdam, Germany. ³University of Potsdam, Institute of Geosciences, Potsdam, Germany. ⁴Vernadsky Institute of Geochemistry and Analytical Chemistry, Russian Academy of Sciences, Moscow, Russia. ⁵Friendly Solutions, Sydney, NSW, Australia. ⁶WiscSIMS Lab, Dept. of Geoscience, Univ. of Wisconsin, Madison, WI, USA. ⁷These authors contributed equally: Adrien Vezinet, Aleksandr V. Chugunov, Alexander V. Sobolev.
✉e-mail: adrien.vezinet@univ-grenoble-alpes.fr; aleksandr.chugunov@univ-grenoble-alpes.fr; alexander.sobolev@univ-grenoble-alpes.fr; stephan@gfz-potsdam.de</footnote>
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) they are similar for present-day uncontaminated MORB and OIB⁸,¹⁴,¹⁵, suggesting mantle homogeneity of these parameters. On the other hand, 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 parents) 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⁴.

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, the subduction of continental crust or/and unmelted oceanic crust to the source of mantle-derived magmas will decrease the Nb/U and Ce/Pb ratios of these melts, leading to an underestimation of the extent of continental crust production. Third, the effect on these ratios of post-emplacement alteration, crustal assimilation, and/or metamorphic overprinting in rocks is difficult to identify and quantify, leading to inaccurate interpretations¹⁰,¹¹,¹⁴.

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 (Fig. 1c) may provide a reliable view of the geochemical composition of the mantle source of melts.

In this study, we discovered an 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. We also conducted advanced geodynamic modeling of the early Earth dynamics and extraction of the continental crust, coupled with the evolution of the 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 the 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 a fluctuating mobile-lid tectonic regime with several tens of million years-long periods of massive subduction induced by mantle plumes.

## Results and discussion

### Weltevreden komatiites

We determined the chemical and isotopic compositions of olivine-hosted melt inclusions in samples from komatiitic flows from the *ca.* 3.27 Ga refs. 22,23 Weltevreden Formation (Barberton Greenstone

Belt, South Africa, Supplementary Fig. 1). Whole-rock isotopic investigations indicate that the mantle sources of these komatiites experienced mantle depletion before the eruption, possibly during the Hadean²⁴ (Supplementary Fig. 2). Further, the high-precision isotopic data analyzed at the whole-rock scale on Weltevreden komatiites indicate the decoupling of the Hf and Nd isotopic systems²⁴. This observation has been interpreted as resulting from elemental fractionation occurring during early Hadean magma ocean solidification²⁵.

Here, using Sr-isotope signatures and trace-element contents measured directly in the homogeneous glass of melt inclusions (Fig. 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.

### Melt inclusions in komatiite olivine phenocrysts

More than 350 olivine-hosted melt inclusions were selected and prepared for geochemical analyses. The analytical workflow, detailed in Methods, Supplementary Tables 1–5 and Supplementary Data 1, 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 analyzed in single stream mode for trace elements. The measured compositions of melt inclusions and host olivines are presented in Supplementary Data 2. Figure 2a, 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 all trace-element contents (shown NiO and Cr₂O₃) of olivine differ both between and within flows. The observed variability suggests that the examined olivine crystals are antecrysts that formed and entrapped their melt inclusions at different sites within a complex magmatic plumbing system. The crystallization of antecrysts occurred in the shallow seawater-altered oceanic crust, as discussed in the subsequent section on contamination.

The TiO₂ contents of Weltevreden melt inclusions show a considerable variation, from ca. 0.14 wt. % to 0.21 wt. %, which is inversely correlated with host olivine Fo contents (Fig. 2c). This correlation follows a trend of *ca.* 33 wt% olivine extraction from the common parental melt (estimated from the ingrowth of Ti from the melts trapped in the most Fo-rich 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 (Fig. 2d–h) or K₂O/TiO₂ (Supplementary Fig. 3b–e), should remain constant in the melt during olivine crystallization because of the negligible contents of these elements in olivine. Clearly, this 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 (<1.5 wt%) amounts of seawater-derived ultra-saline brine before emplacement at *ca.* 3.27 Ga (see Methods): a process common for modern submarine basaltic melts²⁶,²⁷ and larger amounts (e.g., 10–20%) of seawater-altered serpentinite depleted in all incompatible elements similar to olivine cumulates of the same komatiites. This process influences all measured inclusions to varying extents (Fig. 2e), indicating that the crystallization of host olivines occurred at relatively shallow depths within the oceanic crust. Values of Cl/Ti can be used as a proxy for such contamination. A positive correlation is observed
<figure><img src="image_4.png" /><figcaption>**Fig. 1** | Olivine cumulates and olivine-hosted melt inclusions of komatiites from the Weltevreden Formation. a Olivine cumulate from the Weltevreden Formation (sample 1528 C) showing that although these cumulates are significantly altered, they still contain preserved unaltered olivine cores (plane-polarized light). b Partially crystallized natural (not annealed) olivine-hosted melt inclusion from sample 2217 before homogenization at high temperature. Inclusion consists of</figcaption></figure>

glass and acicular calcic pyroxene. c Heated and quenched olivine-hosted melt inclusion from sample 2216 consisting of glass, low-density gas bubble, and spinel (dark crystal). d, e Naturally quenched glassy inclusion in olivine from sample 2218. Inclusion consists of glass, low-density gas bubble, and small spinel and pyroxene crystals.

between Cl/Ti and Rb/Sr (Fig. 2e), suggesting a gain of Rb as well as Na, K, Pb, U, H₂O, and Sr during contamination (Fig. 2f, g, Supplementary Fig. 3). To avoid samples significantly affected by this process, we considered only melt inclusions with low Rb/Sr for geological interpretations in this study (Fig. 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, see Methods). Unlike Rb/Sr, other proxies, such as Ce/Pb and Nb/U, show a 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 and thus can be considered as a minimum estimate of original values.
<figure><img src="image_5.png" /><figcaption>**Fig. 2** | Contents of selected major, minor, and trace elements and Sr isotopes in komatiite melt inclusions (MI) and host olivine crystals. a, b Composition of host olivines, outlines indicate samples from the same flows, Fo = 100Mg/(Mg+Fe) molar. c, d Composition of melt inclusions and host olivine. e–h Composition of melt inclusions. Filled circles are those with a corresponding Sr-isotope analysis: blue for the main group of melt inclusions, yellow for Sr-radiogenic melt inclusions,</figcaption></figure>

and orange for Sr-unradiogenic melt inclusions (see below). Light gray circles are inclusions with no Sr data or rejected Sr isotope data because of Rb/Sr ratio over the maximum value of used reference materials (0.0244) or uncertainty (2SE) of 87Sr/86Sr value > 0.0015 (see Methods). The expected seawater Sr isotope composition at 3.3 Ga after$^{86}$Zr. Uncertainty bars are average individual uncertainty at a 95% confidence interval (c.i.). Source data are provided as a Source Data file.

## 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) 87Sr/86Sr values supported by p(χ²) values <<0.05, indicating statistical heterogeneity within the whole dataset (Fig. 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 the 87Sr/86Sr value of melts and CI/Ti, the chemical proxy of such contamination (see Fig. 2h). This feature is expected because the difference between the 87Sr/86Sr ratios of 3.3 Ga old seawater and Weltevreden komatiite melt is only 0.0020 (Fig. 2h),
<figure><img src="image_6.png" /><figcaption>**Fig. 3** | Sr-isotope of Weltevreden olivine-hosted melt inclusions. a Measured ⁸⁷Sr/⁸⁶Sr ratio. b Age-corrected ⁸⁷Sr/⁸⁶Sr ratio. MSWD - Mean Square Weighted Deviation, n- number of analyses, p(χ²) - chi-squared p-value for the population homogeneity test. Both ⁸⁷Sr/⁸⁶Sr values show low p(χ²) (<<0.05), indicating the lack of statistical homogeneity at the precision obtained for these Sr-isotope analyses.</figcaption></figure>

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

Three statistically homogenous groups can be distinguished in the data shown in (b): unradiogenic in orange, radiogenic in yellow, and neutral (main) in blue (see Methods for Statistical Treatment). Bars represent individual uncertainty at a 95% confidence interval. Source data are provided as a Source Data file.

which is close to the average 95% confidence interval (±0.0010) analytical uncertainty of ⁸⁷Sr/⁸⁶Sr measurements in individual melt inclusions.

A statistical analysis of our Sr-isotope dataset reveals three statistically homogenous groups based on their initial ⁸⁷Sr/⁸⁶Sr (see Fig. 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 Fig. 3b). Olivine antecrysts containing unradiogenic inclusions come from all studied flows (14 of a total of 137, thus 10%), but mostly (8 inclusions out of 14) from Keena's flow 2 (samples 1523B and 2217, Fig. 2a, b, Supplementary Fig. 1), where they compose 14% of 59 inclusions measured for Sr isotopes. The random olivine spinifex (quenched) zone of the same flow exhibits a maximum deviation of μ¹⁴²Nd from terrestrial values, as measured in the Weltevreden flows (Supplementary Fig. 2). This suggests that the component underwent Sm/Nd differentiation during the Hadean eon, potentially contributing to the mantle source of the komatiite from this flow. However, this did not elucidate the process that occurred in Hadean. The mean initial isotopic ratio of unradiogenic inclusions indeed transposes to a Sr-model age, assuming no Rb in the source, of 4.31 ± 0.19 Ga (95% confidence interval, p(χ²) = 0.57, Fig. 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 size of inclusion 2217-48-ol65, belonging to the unradiogenic group, four Sr-isotope & trace-element replicates were conducted. As demonstrated in Fig. 4d-f and Supplementary Table 4, 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.

The Weltevreden komatiites melt is contaminated by seawater-altered crust, as indicated by the composition of olivine-hosted melt inclusions (Fig. 2d-f). This and low-density fluid bubbles in quenched inclusions (Fig. 1c, e) suggest that olivine antecrysts formed in the shallow oceanic crust. The ability of olivine to sample isotopically

diverse mantle-derived melts during crystallization in shallow crust implies that melts from deep heterogeneous mantle sources could be delivered and focused at crustal levels without complete mixing. Recent examples include the shield volcanoes of the Hawaiian plume, where melts formed at depths of around 100 km²⁸ in a single mantle plume are focused in local magmatic systems of Mauna Loa and Kilauea volcanos, which coexist in time and differ in compositions of Sr, Nd, Pb isotopes, and trace elements²⁹. Furthermore, melts with various trace elements and Sr, Pb, and Nd isotope compositions are trapped in olivine crystallized in Mauna Loa volcano's shallow plumbing system¹⁷,³⁰. Finally, the adjacent Weltevreden komatiite flows show significant differences, not only in long-lived isotopic systems Lu-Hf and Sm-Nd but also in the ¹⁴²Nd/¹⁴⁴Nd isotopic ratios (see Supplementary Fig. 2). The observed differences cannot be attributed to alteration and are unlikely to result from contamination, as indicated by the weak correlation between Nd contents in melts and contamination proxies (Supplementary Table 5). Instead, they suggest that the parental melts delivered to the surface were incompletely mixed.

## 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, Fig. 4b). This minimum estimate of the original value, due to contamination (Fig. 2g), is significantly higher than all PM/BSE Nb/U estimates, which range between 26.7 and 32.4⁸,³¹⁻³³. A similar observation can be made for the mean Ce/Pb value of 16.73 ± 1.12 (MSWD = 2.3, *n* = 13, Fig. 4c) of this unradiogenic group, which is also the minimum estimate (Fig. 2f) but is significantly higher than the PM/BSE value of 9–11.2⁹,³¹⁻³³. The four replicates obtained on inclusion 2217-48-ol65 show the same tendency (Fig. 4d-f). The geochemical effect of Ca and Mg-perovskite cumulation proposed for the source of Weltevreden komatiites²⁵ 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 a higher partition coefficient for U than for Nb in the suggested mixture of Ca and Mg-perovskite³⁴. In addition, similar and even much more significant decoupling of εNd
## Main geochemical features of the Sr-unradiogenic sub-population

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

## Main geochemical features of 4 replicate analyses of a single melt inclusion

<figure><img src="image_9.png" /><figcaption>**Fig. 4** | Sr-isotope and trace-element signatures of the Weltevreden olivine-hosted melt inclusions of the Sr-unradiogenic group. a The Sr-model age (see Methods) of the unradiogenic group indicates mid-Hadean extraction from the Bulk Silicate Earth (BSE). b, c The Nb/U and Ce/Pb 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</figcaption></figure>

during the Hadean. d–f The results obtained on 4 replicate analyses of a single melt inclusion (MI) 2217-48-ol65 (belonging to the 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). Source data are provided as a Source Data file.

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”³⁵. The main mantle source component of studied komatiite was highly depleted mantle harzburgite, for instance, similar to refractory harzburgites from Archean cratons³⁶, as follows from their contents of incompatible elements (including Sm-Nd and Lu-Hf), which are lower in parental melt than in BSE mantle¹⁶. Thus, the decoupling of εNd and εHf in their sources can be explained by a mantle restite of the Eo-archean-Hadean age. Our interpretation of the nature of the source of Weltevreden komatiite further substantiates independent observations obtained via thermodynamic modeling³⁷.

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 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 a few tens of million years, which are well within our reported ± 0.19 Ga uncertainty of the model age (see Supplementary Fig. 7 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 minimum fraction of extracted crust to be 0.43% ± 0.09% of the processed mantle. Assuming 2.17 × 10²² kg as the mass of the present-day continental crust and 4.01 × 10²⁴ kg as the mass of the whole mantle, the fraction of continental crust extracted in the Hadean would be 80% ± 16% (Fig. 5, large orange dot) of the mass of the present-day continental crust. A similar calculation for individual melt inclusion 2217-48-ol65 yields 72% ± 24% (Fig. 5, hatched smaller orange dot). Such high values may be overestimated because they assume the involvement of the entire mantle. However, they still provide direct evidence supporting earlier models that argue for significant continental crust production and recycling early in Earth's history³,⁶,⁴².

Our 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's history, are consistent with the composition of ancient zircons⁴³,⁴⁴, the radiogenic isotope compositions of bulk rocks⁴⁵ and numerical modeling (see Geodynamic modeling). Elevated Nb/U and Ce/Pb values obtained in deeply sourced primary Weltevreden komatiitic melts indicate efficient, vertical transfer (subduction of oceanic crust residue after partial melting producing felsic crust) between the upper and the lower mantle in the Hadean.
<figure><img src="image_10.png" /></figure>

**Fig. 5 | Compilation of previously published continental growth curves³ and new estimates for the Hadean from this study (orange circles).** These estimates have been obtained using the measured average Nb/U value in the unradiogenic-Sr group (large orange dot) and individual inclusion 2217-48-ol65 (smaller orange hatched dot) and the average of Nb and U values from four mantle models⁸,³¹⁻³³ and four crustal models³⁸⁻⁴¹, assuming processing of the 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). The light blue field shows the range of continental crust fractions produced during the Hadean, which, according to the geodynamic models, fit geochemical data best (see the Geodynamic modeling section below). Source data are provided as a Source Data file.

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 the actual mass of mantle processed cannot be extracted from such a geochemical dataset. To finalize the analysis, we conducted geodynamic modeling, during which we monitored the geochemical proxies $^{87}$Sr/$^{86}$Sr, Nb/U, and Ce/Pb to determine the tectonic regime that most accurately reflects the measured geochemical data.

## Geodynamic modeling

We applied the mantle convection code StagYY⁴⁶, which has been extensively used over the last decade for modeling the coupled core-mantle-crust evolution of rocky planets. With a 2D spherical annulus domain⁴⁷, our models generate both basaltic and felsic melts⁴⁸, include cooling of the core, and use pressure- and temperature-dependent water solubility maps for different mantle minerals¹². For this study, we improved previous models¹²,⁴⁸ by including the effect of water on the density of mantle materials, incorporating a composite rheology (diffusion creep and dislocation creep proxy) for the crust and upper mantle based on experimental data, and by initializing the frictional strength of the early Earth's oceanic lithosphere. Following approaches from previous studies⁴⁹,⁵⁰, 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 the lithosphere's hydration and frictional strength (see Methods for model setup).

Our numerical experiments start with a homogeneous solid Earth of pyrolytic (BSE) composition³¹ at 4.5 Ga with a low (0.01 wt. %) initial water content¹², a mantle potential temperature of 1900 K, and a core temperature of 5000 K. They end after 1.5 billion years, i.e. at 3.0 Ga. Any cataclysm relating to the formation of the mantle or the Moon is assumed to be earlier than 4.5 Ga.

All models fall into two distinct groups that illustrate different tectonic regimes, as we show below Group I consists of models with

our preferred effective friction coefficient of the lithosphere for the early Earth of 0.1⁵¹, while Group II consists of the models with higher lithospheric strength (effective friction coefficient of 0.2). To better demonstrate the effect of the tectonic regime on geochemical parameters, we have selected parameters controlling surface water input in the models of both groups, ensuring they produce a similar amounts of continental crust. To achieve this, the high lithospheric strength models of Group II needed significantly more water input compared to those of Group I. The representative models from both groups are shown in Fig. 6. The typical behavior of all models is intense mantle convection during the first hundred million years (Fig. 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. During this stage, there is extensive production of both oceanic and continental crust (Fig. 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⁵². A mobility value above 1 means that the lithosphere moves horizontally relative to the underlying mantle with significant velocities, and subduction is active¹³.

In models of Group I, high mobility periods last several tens of million years and interchange with the low mobility periods (Fig. 6a). We call this tectonic regime a fluctuating 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 (Fig. 6d2). The plumes 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 reduced production of new crust (Fig. 6b, model time range 480-680 Myr). 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 breakthrough, approach the lithosphere, and induce new subduction zones, which produce and recycle new crust (Fig. 6d3).

Group II models are characterized by low mobility throughout their evolution and experience several mobility bursts that last a few million years (Fig. 6e). This type of tectonic regime is usually called an episodic-lid regime¹³. During these episodes of high mobility, the over-thickened lithosphere is rapidly subducted, resulting in resurfacing events (Fig. 6g3). Most of the continental crust is produced during these events.

In all models, at every time step, we identify regions of hot material (mantle plumes with 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 the Hadean fraction of material with the $^{87}$Sr/$^{86}$Sr ratio <0.6997. The evolution of the trace-element ratios in the Hadean fraction for a number of models from both groups is presented in Figs. 7a, d, 8a, d, and Supplementary Fig. 10.

The common feature of all models is that trace elements do not show any effect of continental crust production on the composition of mantle plumes till at least 500 Myr of model time (i.e., 4.0 Ga) even though the continental crust production during the first 500 Myr (i.e., 4.5-4.0 Ga) is more than twice that of the last 1000 Myr (i.e., 4.0-3.0 Ga) in most models. This time lag indicates the time required by the restites of hydrated basaltic crust, which were involved in the production of continental crust, to recycle and mix in substantial amounts with the plume source material in the lower mantle.

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, Fig. 5,
<figure><img src="image_11.png" /><figcaption>**Fig. 6** | Representative models from Group I and Group II. a, e Mobility of models Group I and Group II respectively over time with a moving average (thick line) over 11 mobility values (thin lines). b, f Mass of generated continental crust over time for models Group I and Group II, respectively, scaled with the present-day continental crust mass (CCM). c, d1–d4, and g1–g4 Snapshots of Earth's cross-section showing models' evolution with time (c, d1–d4 for Group I and e, g1–g4 for Group II), where four quadrants represent different fields (clockwise from top left: composition, Nb/</figcaption></figure>

U, viscosity, $^{87}$Sr/$^{86}$Sr). The color scales for each field are shown in (c). Snapshot
c displays the initial model state, the same for the models of both groups, corre-
sponding to an age of 4500 Ma and a model time of 0 Myr. Snapshots d1, d3,
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).
Source data are provided as a Source Data file.

and Fig. 7a, d). These models also predict a realistic amount of Hadean
component in plumes (Fig. 7c, f) and a rather moderate supply of water
into the lithosphere (0.5–0.9 present-day ocean mass (OM)) over the
period of 1500 Myr (Supplementary Fig. 11). In contrast, the models of
Group II where we chose the higher water input parameters show a
poor fit with Nb/U and Ce/Pb data (Fig. 8a, d), a too low amount of

Hadean component in plumes (Fig. 8c, f) and predict a water supply of
1.1–2.4 OM into the lithosphere over the period of 1500 Myr (Supple-
mentary Fig. 11). 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.
Group I: Models with high subduction activity

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

**Fig. 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 the mean trace-element ratio for the selected cells with uncertainty bars with two standard deviations of the mean value. **b, e** Mass of generated continental crust over time from the same models, scaled with the

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

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 uncertainty bars of 2 standard deviations of the mean value. Model parameters: f–effective friction coefficient of the lithosphere, alpha–water saturation fraction of the top 5 km of crust (see Methods section). Source data are provided as a Source Data file.

The difference between predictions by models of Group I and II is also evident from the calculated Sr isotope composition of potential komatiite sources (Fig. 9). The models with high subduction activity (models of Group I) accurately reproduce the mean strontium isotope composition of melt inclusions at the time of emplacement of Weltevreden komatiites, along with its standard deviation (Fig. 9a, b), reflecting the isotopic heterogeneity of potential melt sources. We assume that melts possess identical Sr isotope ratios to their respective melt sources. In contrast, models with low subduction activity (Group II) are less successful (Fig. 9c, d).

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–70% of CCM production during the Hadean in agreement with maximal estimates of 50 to 100% of CCM extraction from the whole mantle using a mass balance approach (Fig. 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 regimes 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 rapid, variable and short-lived. Moreover, our models show that in

Hadean-Archean Earth, in contrast to present-day Earth, large mantle plumes are of key importance to trigger such subduction, i.e. so-called plume-induced subduction⁵³, and when plume activity is diminished, the subduction activity is also diminished.

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³,⁶.

## Methods

### Methods for geochemical investigation

The extended analytical workflow developed for the geochemical study of Weltevreden melt inclusions is broadly similar to that presented in recently published geochemical studies of melt inclusion from Song Da ultramafic volcanic suite⁵⁴,⁵⁵.

### Samples

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 (Supplementary Fig. 1). Three samples (1521, 1522, and 1523) from the Saw Mill Complex were previously described¹⁶. The new samples 2216, 2217, 2218, and 2222 from Saw Mill and 1528 from Pioneer Complexes are similar to those in composition. These cumulates
Group II: Models with low subduction activity

<figure><img src="image_14.png" /><figcaption>**Fig. 8 | A subset of results of models of Group II compared with observations.**
The legend and panel description are the same as in Fig. 7. Note the poor fit of the geochemical data (red and blue diamonds) compared with the models of Group I</figcaption></figure>

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

shown in Fig. 7. Uncertainty bars are 2 standard deviations of the mean value.
Source data are provided as a Source Data file.

consist of partially altered (serpentine and magnetite) euhedral olivine crystals (more than 60 volume %) and interstitial groundmass consisting of spinel, magnetite, and a variety of pyroxenes (low-Ca, high Ca-pyroxenes, and pigeonite) enclosed in a matrix of altered glass (chlorite and magnetite). In the groundmass of sample 2222, we also found clinoenstatite with typical polysynthetic twinning and tremolite in intergrowth with orthopyroxene. The new samples are: 2216 (a few meters South of Keena's flow 1, 25.842597°S, 30.885031°E, olivine composition 92.8–93.1 mol.% Fo); 2217 (Keena's flow 2, 25.841475°S, 30.885683°E, olivine composition 93.2–93.8 mol.% Fo); 2218 (Keena's flow 1, 25.841803°S, 30.886147°E, olivine composition 94.6–95.1 mol.% Fo); 2222 (Gary's flow 2, 25.843578°S, 30.887742°E, olivine composition 94.8–95.2 mol.% Fo); and 1528 (Fig. 1a) from the Pioneer Complex (25.834308°S, 30.949083°E, olivine composition 92.5–93.5 mol.% Fo). Olivine crystals contain abundant melt inclusions, which are partly crystallized if they exceed 20 micrometers in diameter (Fig. 1b) or glassy if they are smaller (Fig. 1e).

Melt inclusions

Prior to polishing and mounting, melt inclusions hosted in olivine crystals were homogenized at high temperatures and 1 bar in the CO₂-H gas mixture corresponding to ΔQFM-1 (quartz-fayalite-magnetite oxygen buffer). 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 the previously established protocol¹⁸. Target temperatures for experiments range from 1280 °C to 1350 °C in order to complete the 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 that 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 the detection limit and were neglected in this study. Correction of the measured composition of glass in melt inclusions (see below) for post-entrapment olivine crystallization on the walls of inclusion's cavity and Fe-Mg exchange with host olivine (PEC) was performed using Petrolog 3 software⁵⁶. Details on this approach are reported in earlier studies¹⁶. PEC was used to correct the inclusion compositions presented in Supplementary Data 2 and 3.

Electron microprobe analyses

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

Raman microscope

The water content of melt inclusions was determined using a Raman microscope LABRAM Soleil, Horiba (ISTerre, Grenoble, France), equipped with 473 nm and 532 nm lasers and an optical microscope (Nikon Eclipse LVDIA-N) using an analytical protocol previously published⁵⁵. Average and 2SE uncertainties are reported in Supplementary Data 2–4.
<figure><img src="image_16.png" /><figcaption>**Fig. 9** | Calculated Sr isotope ratio in potential komatiite mantle sources versus time in typical models of Group I and Group II compared with observations. **a, b** Typical models of Group I. **c, d** Typical models of Group II. The observed mean Sr isotope ratio of melt inclusions in olivines in Weltevreden komatiites is indicated by red crosses. Blue solid circles indicate modeling results, and each circle represents the calculated mean isotope ratio for the selected mantle cells of hot plume</figcaption></figure>

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

material (see Methods for details). The solid blue line shows the Bulk Solid Earth (BSE) evolution trend. Both modeling results and data are shown with uncertainty bars with 2 standard deviations of the mean value. Model parameters: f-effective friction coefficient of the lithosphere, alpha-water saturation fraction of the top 5 km of crust (see Methods section). Source data are provided as a Source Data file.

## 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 laser ablation (LA) system (RESOLution SE, Applied Spectra) coupled with a multi-collector (Neptune XT, Thermo-Fisher Scientific) for Sr-isotope and a single-collector (8900, Agilent) inductively-coupled plasma mass spectrometer (ICP-MS), for TE. LA single stream and LASS protocols were used depending on the size of the melt inclusions: smaller inclusions were analyzed with LASS Sr-TE analytical protocol, while bigger ones had two laser analyses, one for Sr-isotopes and one for TEs. The smallest melt inclusions with diameter less than 30 micrometers, were analyzed only for TE. Uncertainties are reported at a 95% confidence interval unless stated otherwise. It is worth mentioning that Sr content in Weltevreden melt inclusions averages 36 ± 6 ppm, corresponding to a mass of total Sr of ca. 6 pg 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 this approach being 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 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 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¹¹ Ω amplifiers. The different dynamic response time between 10¹¹ Ω amplifiers and 10¹³ Ω amplifiers was accounted by using the Thermo-Fisher software built-in tau correction⁶⁰. 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 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 signals of masses 83.5 and 86.5 were always below the 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, 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 for each individual analysis based on the measured (88/86) value and the natural value (8.37520938). The Rb fractionation factor (βRb) was manually adjusted to provide the best accuracy for all three reference materials measured throughout the analytical run. This approach allows for accounting for day-to-day variations in βRb and is analogous to the approach used for 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 by 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 labeled 15XX, adjusting βRb value on a daily basis was sufficient to obtain accurate results (within 2 SD) for all three reference glasses covering all compositional ranges of measured inclusions used in each analytical run (KL2-G, NIST-614-G, and GOR128-G, Supplementary Fig. 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 labeled 22XX, standardization of both ⁸⁷Sr/⁸⁶Sr and ⁸⁷Rb/⁸⁶Sr was conducted with KL2-G as a calibration glass.
Readers can refer to the Sr-isotope DRS provided on the Iolite Github (https://github.com/iolite-LA-ICP-MS) for a 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 backgrounds, interference corrections, and standardization conducted. Note that including the REE²⁺ correction changes the final uncertainty to the sixth decimal only, i.e., two orders of magnitude lower than the 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 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 the 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 Data 1), while an inter-sequence excess variance of 5.5%, 6%, and 6% was propagated into ⁸⁷Rb/⁸⁶Sr, Rb/Zr and Sr/Zr values, respectively, for all analyses (reference glasses and melt inclusions). Inter-sequence excess variance likely reflects the contribution of elemental fractionation, which is unaccounted for in our data reduction protocol. Hence, the reported individual uncertainties, which include excess variance in ⁸⁷Sr/⁸⁶Sr, ⁸⁷Rb/⁸⁶Sr, and Sr/Zr values, are considered external uncertainties. Precision and accuracy of Sr isotope analyses of melt inclusions were double-checked by duplicate, triplicate, or quadruplicate analyses of sufficiently large inclusions (Supplementary Fig. 5 and Supplementary Table 3). The homogeneity of the replicates of 2, 3, or 4 analyses of the same melt inclusion presented in this diagram and table confirms that reported uncertainties of individual inclusions are not underestimated.

All melt inclusions discussed in this work have ⁸⁷Rb/⁸⁶Sr values lower than the ⁸⁷Rb/⁸⁶Sr of KL2-G (0.0699, or 0.0248 as Rb/Sr), the reference glass with the most elevated ⁸⁷Rb/⁸⁶Sr. Melt inclusions displaying ⁸⁷Rb/⁸⁶Sr higher than KL2-G were disregarded since the accuracy of the ⁸⁷Rb correction could not be assessed. In addition, inclusions with 2 standard errors of measured ⁸⁷Sr/⁸⁶Sr ratio over ±0.0015 (thin inclusions where the signal was collected for less than 8 seconds) were also excluded. Lastly, we checked the robustness of ⁸⁷Sr/⁸⁶Sr uncertainty in accepted individual inclusions by comparison with the uncertainties in the reference materials with similar content of Sr (GOR128 and NIST 614) using the following equation:

$$ \sigma_{\text{inclusion}} = \sigma_{\text{RM}} \times \sqrt{\frac{Sr_{\text{RM}} \times t_{\text{RM}}}{Sr_{\text{inclusion}} \times t_{\text{inclusion}}}} \quad (1) $$

Where $\sigma_{\text{inclusion}}$ and $\sigma_{\text{RM}}$ are 95% c.i. uncertainties of individual inclusions and reference material (glasses), respectively; $Sr_{\text{inclusion}}$ and $Sr_{\text{RM}}$ are Sr contents in inclusions and reference materials; $t_{\text{inclusion}}$ and $t_{\text{RM}}$ are ⁸⁷Sr/⁸⁶Sr acquisition times for inclusions and reference materials (Supplementary Table 1).

Predicted average $\sigma_{\text{inclusion}}$ calculated for all Weltevreden melt inclusions and for the 14 unradiogenic ones are 0.0009 ± 0.0001, as well as their average measured $\sigma_{\text{inclusion}}$, are 0.0010 ± 0.0001, thus virtually identical (Supplementary Table 1). We, therefore, conclude that reported uncertainties in ⁸⁷Sr/⁸⁶Sr values 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⁻². Carrier gas was He (ca. 0.3–0.9 L/min) with the addition of N₂ (ca.1 to 2.5 mL/min), which was mixed with Ar (ca. 0.6-1.1 L/min) before introduction into the spectrometer. The oxide production rate, monitored with ThO²/Th⁺, was <0.1%, and the doubly charged ratio monitored with ⁴⁴Ca²⁺/⁴⁴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 ablation to eliminate surface contamination, 30-40 s background measurement followed by 30–40 seconds 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 dwell times and TE data reduction are reported in Supplementary Table 2. Ce/Pb, Nb/U, and Rb/Sr of reference glasses (GOR128-G⁶⁴ and NIST-614-G⁶⁵) are shown in Supplementary Fig. 4 and reported in Supplementary Data 1. Reference 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 an 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⁶⁶. Elevated MSWD ratios are common in Earth science geochemistry, and information from ‘dispersed’ datasets has geological meaning⁶⁷. This significant heterogeneity is also clearly observed for (⁸⁷Sr/⁸⁶Sr)initial values that show an MSWD of 2.5, translating to a p(χ²) of 0. Further, it is also clear from Supplementary Fig. 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%) comprised of 14 inclusions (MSWD = 0.88, p(χ²) = 0.57), (ii) one with intermediate (⁸⁷Sr/⁸⁶Sr)initial = 0.700624 ± 0.000086 made of 116 inclusions (MSWD = 1.1, p(χ²) = 0.27) and one with elevated (⁸⁷Sr/⁸⁶Sr)initial = 0.70177 ± 0.00024 made of 7 inclusions (MSWD = 0.48, p(χ²) = 0.82).

Applying the t-test for independent populations with unequal variances⁶⁸ reveals that the means of (⁸⁷Sr/⁸⁶Sr)initial of all groups are different with a confidence level over 99.9% (see Supplementary Table 6).
## Model age calculation

Sr-model age of olivine-hosted melt inclusions was based on (i) ($^{87}$Sr/$^{86}$Sr)initial and (ii) 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 the uncertainty of our estimates. The calculation was conducted assuming no Rb left in the BSE source after melt extraction. This assumption is realistic because the source of Weltevreden komatiite has likely undergone significant melting and several stages of efficient melt extraction, resulting in severe depletion of most incompatible elements¹⁶. The relatively high Rb/Sr ratios of trapped melts (0.014–0.024 for inclusions measured for $^{87}$Sr/$^{86}$Sr) are apparently the result of their contamination during and before the crystallization of host olivine (Fig. 2e). Samarium–Neodymium (Sm-Nd) and Lutetium–Hafnium (Lu-Hf) model ages presented in Supplementary Fig. 2 were calculated using Ryan Ickert’s spreadsheet⁷⁰.

## Contamination of Weltevreden komatiites

Studied inclusions in olivine yield significant intercorrelated variations of Cl, K, Rb, H₂O, Pb, Sr, U, and Na, which are unrelated to olivine crystallization and extraction from the melt. This is well demonstrated by the strong correlation between ratios of elements excluded from the olivine structure, e.g., Cl, K, Rb, H₂O, Pb, Sr, Na, and Ti, which must stay constant during olivine fractionation (Supplementary Fig. 3). Similar variations are common for submarine glasses of OIB and MORB (Supplementary Fig. 3a) and have been explained by minor (less than 1 wt%) contamination of seawater-derived ultra-saline brines²⁶,²⁷. 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.20 wt% of Cl to Weltevreden melt is about or less than 1.5 wt%. We also cannot exclude the assimilation of larger amounts of serpentinite. This process does not affect major element concentrations considering their analytical uncertainty but is enough to be observed in trace elements and H₂O contents. Also note that the extent of contamination tends to increase with decreasing host olivine Fo content (Fig. 2d). The initial $^{87}$Sr/$^{86}$Sr isotope ratio of melt inclusions does not correlate with contamination proxies (see Fig. 2h); thus, we conclude that the effect of observed contamination is within the analytical uncertainties in $^{87}$Sr/$^{86}$Sr values.

## Methods for geodynamic modeling

**Essentials and setup.** We study the thermochemical evolution of compressible mantle using the code StagYY⁴⁶, which has been developed and widely used over several decades for global-scale modeling of Earth’s evolution spanning its age¹²,¹³,⁴⁸–⁵⁰. The model includes pressure- and temperature-dependent viscosity, plastic yielding, time-dependent radiogenic heating, core cooling, and phase changes. Following a two-step mantle differentiation and utilizing water-dependent solidi functions¹², the code forms basaltic and felsic melts and considers both intrusive (plutonic) and eruptive (volcanic) magmatism⁴⁸. 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⁴⁸,⁷¹. We employ temperature- and pressure-dependent water solubility maps for different mantle minerals, which control the water in-gassing and out-gassing¹².

We use a two-dimensional quadrant spherical annulus geometry⁴⁷ with an 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) the lower mantle (2), and the post-perovskite layer (3), with each layer having different values for activation energy $E_i$ and activation volume $V_i$⁷²,⁷³. 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 (2) $$

where $η_0$ is the reference viscosity (2.35 10¹⁹ Pa s) at zero pressure and reference temperature $T_0$ (1710 K), $Δη_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 (3) $$

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 $σ_Y$, which has both brittle and ductile components:

$$ \sigma_Y = \min(\sigma_{Y, \text{ductile}}, \sigma_{Y, \text{brittle}}) \quad (4) $$

The ductile yield stress $σ_{Y, ductile}$ increases linearly with pressure as:

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

where $σ_Y^0$ is the surface ductile yield stress and $σ_Y$ is the pressure gradient of the ductile yield stress. The brittle yield stress $σ_{Y, brittle}$ is calculated following the Byerlee approach as:

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

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

See Supplementary Table 7 for model parameters used in this study. The model uses a parameterization based on mineral physics data⁷⁵,⁷⁶, 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 initialized 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 ref. 12. Here, we describe only the model modifications introduced in this study.

## Models modifications

**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⁻¹⁵ 1/s and using dry olivine rheology parameters⁷⁷ (Supplementary
Table 7). We also considered frictional plasticity with a coefficient of friction of 0.1 corresponding to the properties of subduction channels, which lack continental sediments⁵¹. A similar effective friction coefficient was suggested based on the model of thermal cracking of oceanic lithosphere on early Earth⁷⁸.

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

**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 the water saturation fraction $\alpha_w$ parameter ($0 < \alpha_w < 1$) controlling the input of the surface water into the lithosphere.

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

$$ \rho_{\omega} = \rho_{d}(1 - aC_{H2O}) \quad (7) $$

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⁷⁹ (Supplementary Fig. 9) that we consider as 2.0 in this study.

**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 the maximum degree of melting of basalt to produce TTG is 30% instead of 10% used in previous studies¹²,⁴⁸.

**Isotopic systems and trace elements.** We have introduced the evolution of the Rb-Sr isotope system and the Nb, U, Ce, and Pb trace elements. The constant parameters and solid-melt partition coefficients used in the study are presented in Supplementary Tables 9 and 10, respectively. While we use the same technique to compute isotope and trace-element evolution in the StagYY code as done previously⁴⁹,⁵⁰, our models have several essential differences. First, the StagYY code has significantly evolved since 2004, and a key feature in our models is the production of felsic melts. Continental crust production drastically affects the evolution of many isotopic systems and trace elements. Second, the-Sr isotope system was not included in previous studies. Our modeling purpose is also different. We aim to compare the modeling results with geochemical data reporting the trace-element composition of the Hadean melt source in the Archean komatiites. To achieve this, we analyze the compositions of potential sources of komatiitic melts at the depth before these melts are separated from the source rocks and before the melts from different sources are mixed.

### Definition of sampling for Figs. 7–9

For the models presented in Figs. 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 100 K higher than the average temperature of all cells located at the same depth. This condition allows us to identify mantle plume material. In practice, the potential temperature of all selected cells in most of the models appears to be higher than 1700 °C. (ii) Cells have to be located within the depth range of 400–600 km. This condition is required to avoid cells from the top 300 km, from which the melt is extracted, and trace elements are partitioned. (iii) The average initial Sr isotope ratio ($^{87}$Sr/$^{86}$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-element 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 a value of 2 standard deviations of the mean Nb/U ratio (≥4), and with too large a value of 2 standard deviations of the 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 model output times, which we show in Figs. 7 and 8 of the main text.

To avoid misinterpretation, the sampling conditions behind modeling results presented in Figs. 7 and 8 require further elucidation. The rationale is that our geodynamic model, like others, currently fails to replicate the presence of primary melts from various mantle sources at the crustal level without their mixing, a process accomplished through the geochemical methodology outlined in this paper. However, we can detect a potential Hadean melt source considering numerical cells in the mantle with a high enough potential temperature (condition “i”) and Hadean Sr modal age (condition “iii”). Because the extraction of melt changes isotope and trace-element composition in the heterogeneous source due to different fusions of heterogeneities, we have to consider numerical cells at the depth before the melt is extracted. Therefore, we consider cells at a depth below 400 km (condition “ii”).

In the models depicted in Fig. 9 (main text), we apply cell selection criteria (i) and (ii), thereby sampling potential komatiite sources across all Sr-model ages. For presentation in Fig. 9, we also filter out values of Sr isotope ratios at times when there were less than 200 “plumes” cells in the mantle (i.e., low plume activity). The good fit of the simulated Sr isotope model ages to the geochemical data for Group I models (Fig. 9) further substantiates the legitimacy of condition (ii) for the selection of numerical cells to simulate the sources of Weltevreden komatiites.

Due to the 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 Fig. 7 comprises two sets of models (3 in each set) with the same frictional strength ($f=0.1$) but with different water saturation fractions ($\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 Fig. 8, also comprises two sets of models (2 in each set) with the same frictional strength ($f=0.2$) but with different water saturation fractions ($\alpha=0.40$ or $0.30$). Values of $\alpha$ were chosen such that all models of both groups produced a similar amount of continental crust during the Hadean.

### Explanation of the Nb/U and Ce/Pb trends in Figs. 7 and 8

To better explain trace-element trends in Figs. 7 and 8a, d, we consider the following end-member cases of continental crust production and recycling of the restites.

**Case 1:** No continental crust is produced in the Hadean. Obviously, in this case, Nb/U and Ce/Pb ratios in the Hadean fraction will not change from the BSE values.

**Case 2.** Continental crust is produced in the Hadean, but restites are not recycled and remain in the lithosphere. In this case, again, Nb/U and Ce/Pb ratios in the Hadean fraction sampled by the plume will not change from the BSE values.
Case 3. Continental crust is produced in the Hadean, and restites of the Hadean age are immediately recycled into the lower mantle and mixed into the plume sources. In this case, Nb/U and Ce/Pb ratios in the Hadean fraction sampled by the plume will increase proportionally to the production of the continental crust during Hadean time and will stay constant after Hadean time.

Case 4. Continental crust is produced in the Hadean, and restites of the Hadean age are recycled into the lower mantle and mixed into the plume sources not immediately but continuously during some time, say 1 Gyr. In this case, Nb/U and Ce/Pb ratios in the Hadean fraction will increase continuously over time, even after Hadean.

Our models of Group I (Fig. 7a, d) follow case 4, and models of Group II (Fig. 8a, d) follow case 2 with a minor addition of case 4.

## Parameter sensitivity

In the main text, we show the effect of lithospheric hydration and strength on the evolution of the trace-element ratios of the Hadean fraction of the potential komatiite source (Fig. 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 Supplementary Fig. 10 show our models that fit well with the geochemical observations computed with our preferred value of "a" of 2. Panels i-p in Supplementary Fig. 10 show modeling results obtained with the value of "a" parameter of 1.5, which is more suitable for rocks at relatively high pressures (Supplementary Fig. 9b) with all the other parameters the same as in the models with a=2. As we see from Supplementary Fig. 10, the variation of "a" parameter does not significantly change the modeling results.

## Models' limitations

It is clear that even advanced geodynamic models of the early Earth's evolution, which include the generation and recycling of oceanic and continental crusts, cannot avoid simplifications, and they have certain limitations. The main simplification of our global models is their 2D geometry. Although a 2D spherical annulus model setup replicates a 3D spherical case much better than a 2D cylindrical or a 2D axisymmetric model⁴⁷, some significant differences still remain. For instance, 2D plumes and subduction zones, by definition, have a global character, which leads to an exaggerated melt production and recycling of both oceanic and continental crust. As these artifacts accumulate over time, they are likely more significant for the modeling of the full 4.5 Gyr evolution of the Earth and are less important for the models of the early Earth evolution. Nevertheless, one possible artifact of our models is that most of the generated TTG is recycled back to the mantle during the 1.5 Gyr modeling time. Although the high rate of continental crust recycling in Hadean and Eo-Archean is also suggested in some geochemical models⁶, our models likely overestimate the TTG recycling rate. Despite the possible overestimation of the TTG recycling rate, the amount of recycled TTG in our models is typically less than 0.1% in the Hadean fraction of the source (Sr ratio <0.6997). Such a low amount of recycled TTG does not significantly affect the modeling results for Nb/U and Ce/Pb ratios in the Hadean melt source, and hence, our interpretation of geochemical data is valid.

In our simulations, we utilize thermodynamic models of Earth materials. Models such as Perple_X⁸⁰ are well constrained by experimental data for the uppermost mantle pressure-temperature conditions but may not be as reliable for the deeper mantle, particularly for water-containing phases. The water solubility map based on the Perple_X model for MORB composition shows higher water content at low temperatures and pressures above 5 GPa compared to maps based on experiments by ref. 81, as presented in ref. 82. The Perple_X calculations used in our models allow for a maximum of 8 wt% saturated water content in hydrous phases of the oceanic crust at pressures exceeding

5 GPa. However, in the models presented in this paper, the water content in the oceanic crust never exceeds 6 wt% and only rarely surpasses 4 wt% due to the limited input of water from the surface. The water content in the upper mantle is also quite realistic, at a few hundred ppm in the Hadean and Eo-Archean mantle (Supplementary Fig. 11), similar to the estimates by ref. 82 and ref. 83 for the present-day Earth.

Our geodynamic model, like others, cannot reproduce the presence of primary melts from various mantle sources at the crustal level without mixing them. This is achieved through the geochemical approach presented in this paper. Therefore, we must employ the modeling strategy detailed in the previous sections to replicate the geochemical observations, addressing not the compositions of melts, but rather their sources at depth.

## Data availability

All source data are provided with this paper and are available in the EarthChem library: https://doi.org/10.60520/IEDA/113703 Source data are provided with this paper.

## Code availability

The convection code StagYY is the property of Paul Tackley and Eidgenössische Technische Hochschule Zürich (ETH) and is available for collaborative studies from Paul Tackley (paul.tackley@eaps.ethz.ch).

## References

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

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

3. Korenaga, J. Crustal evolution and mantle dynamics through Earth history. *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* **376**, 20170408 (2018).

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).

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).

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).

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).

8. Hofmann, A., Class, C. & Goldstein, S. Size and composition of the MORB+ OIB mantle reservoir. *Geochem. Geophys. Geosyst.* **23**, https://doi.org/10.1029/2022GC010339 (2022).

9. Boyet, M., Garçon, M. & Avice, G. in *Treatise on Geochemistry (3rd Edition)* Vol. 1 (ed Elsevier) 73-118 (2025).

10. Vezinet, A., Pearson, D. G. & Thomassot, E. Effects of contamination on whole-rock isochrons in ancient rocks: a numerical modelling approach. *Lithos* **386-387**, 1060-1040 (2021).

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

12. Jain, C., Rozel, A. B., van Hunen, J., Chin, E. J. & Manjón-Cabeza Córdoba, A. Building Archean cratonic roots. *Front. Earth Sci.* **10**, 966397 (2022).

13. 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).
14. 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).

15. Reimink, J., Davies, J., Moyen, J. & Pearson, D. A whole-lithosphere view of continental growth. *Geochem. Perspect. Lett.* **26**, 45–49 (2023).

16. 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).

17. 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).

18. Sobolev, A. V. et al. Komatiites reveal a hydrous Archaean deep-mantle reservoir. *Nature* **531**, 628–632 (2016).

19. 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).

20. 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).

21. McDonough, W. & Ireland, T. Intraplate origin of komatiites inferred from trace elements in glass inclusions. *Nature* **365**, 432–434 (1993).

22. 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).

23. 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).

24. Puchtel, I. et al. Insights into early Earth from Barberton komatiites: evidence from lithophile isotope and trace element systematics. *Geochem. Cosmochim. Acta* **108**, 63–90 (2013).

25. Puchtel, I. S., Blichert-Toft, J., Horan, M. F., Touboul, M. & Walker, R. J. The komatiite testimony to ancient mantle heterogeneity. *Chem. Geol.* **594**, 120776 (2022).

26. Coombs, M. L., Sisson, T. W. & Kimura, J.-I. Ultra-high chlorine in submarine Kilauea glasses: evidence for direct assimilation of brine by magma. *EPSL* **217**, 297–313 (2004).

27. 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. *Geochem. Cosmochim. Acta* **123**, 150–165 (2013).

28. Herzberg, C. T., Asimow, P. D. & Hernández-Montenegro, J. D. The meaning of pressure for primary magmas: New insights from PRI-MELT3-P. *Geochem. Geophys. Geosyst.* **24**, e2022GC010657 (2023).

29. Pietruszka, A. J., Norman, M. D., Garcia, M. O., Marske, J. P. & Burns, D. H. Chemical heterogeneity in the Hawaiian mantle plume from the alteration and dehydration of recycled oceanic crust. *EPSL* **361**, 298–309 (2013).

30. 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).

31. McDonough, W. F. & Sun, S.-s The composition of the Earth. *Chem. Geol.* **120**, 223–253 (1995).

32. Hofmann, A. W. Chemical differentiation of the Earth: the relationship between mantle, continental crust, and oceanic crust. *EPSL* **90**, 297–314 (1988).

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

34. 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. *Geochem. Cosmochim. Acta* **69**, 485–496 (2005).

35. 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).

36. Kamber, B. S. & Tomlinson, E. L. Petrological, mineralogical and geochemical peculiarities of Archaean cratons. *Chem. Geol.* **511**, 123–151 (1999).

37. Tomlinson, E. L. & Kamber, B. S. Depth-dependent peridotite-melt interaction and the origin of variable silica in the cratonic mantle. *Nat. Commun.* **12**, 1082 (2021).

38. Wedepohl, K. H. The composition of the continental crust. *Geo-chim. Cosmochim. Acta* **59**, 1217–1232 (1995).

39. Rudnick, R. L. & Fountain, D. M. Nature and composition of the continental crust: a lower crustal perspective. *Rev. Geophys.* **33**, 267–309 (1995).

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

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

42. Armstrong, R. L. Radiogenic isotopes: the case for crustal recycling on a near-steady-state no-continental-growth Earth. *Philos. Trans. R. Soc. Lond. Ser. A* **301**, 443–472 (1981).

43. Harrison, T. M. *Hadean earth*. (Springer Cham, 2020).

44. Bauer, A. et al. Hafnium isotopes in zircons document the gradual onset of mobile-lid tectonics. *Geochem. Perspect. Lett.* **14**, 1–6 (2020).

45. Tusch, J. et al. Long-term preservation of Hadean protocrust in Earth's mantle. *Proc. Natl. Acad. Sci.* **119** https://doi.org/10.1073/pnas.2120241119 (2022).

46. 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).

47. Hernlund, J. W. & Tackley, P. J. Modeling mantle convection in the spherical annulus. *Phys. Earth Planet. Inter.* **171**, 48–54 (2008).

48. 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).

49. Xie, S. & Tackley, P. J. Evolution of helium and argon isotopes in a convecting mantle. *Phys. Earth Planet. Inter.* **146**, 417–439 (2004).

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

51. Sobolev, S. V. & Brown, M. Surface erosion events controlled the evolution of plate tectonics on Earth. *Nature* **570**, 52–57 (2019).

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

53. Gerya, T. V., Stern, R. J., Baes, M., Sobolev, S. V. & Whattam, S. A. Plate tectonics on the Earth triggered by plume-induced subduction initiation. *Nature* **527**, 221–225 (2015).

54. Kazzy, C. et al. Strontium isotope and trace element compositions of olivine-hosted melt inclusions from the Song Da ultramafic volcanic suite, northern Vietnam: Implications for chemical heterogeneity in mantle plumes. *Chem. Geol.* **122564** https://doi.org/10.1016/j.chemgeo.2024.122564 (2025).

55. 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
https://doi.org/10.1016/j.chemgeo.2024.122219 (2024).

56. Danyushevsky, L. V. & Plechov, P. Petrolog3: Integrated software for modeling crystallization processes. *Geochem. Geophys. Geosyst.* **12** https://doi.org/10.1029/2011GC003516 (2011).

57. Falloon, T. J. & Danyushevsky, L. V. Melting of refractory mantle at 1-5, 2 and 2-5 GPa under anhydrous and H2O-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).

58. 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).

59. Zhu, L. et al. Improved in-situ determination of sr isotope ratio in silicate samples using LA-MC-ICP-MS and its wider application for fused rock powder. *J. Earth Sci.* **31**, 262–270 (2020).

60. Craig, G. et al. Dynamic time correction for high precision isotope ratio measurements. *Thermo Fish. Sci. Bremen*, 1–8 https://assets.thermofisher.com/TFS-Assets/CMD/Technical-Notes/TN-30396-ICP-MS-Dynamic-Time-Correction-Neptune-Plus-TN30396-EN.pdf (2017).

61. Paton, C., Hellstrom, J., Paul, B., Woodhead, J. & Hergt, J. Iolite: freeware for the visualisation and processing of mass spectrometric data. *J. Anal. Spectrom.* **26**, 2508–2518 (2011).

62. Currie, L. A. Limits for qualitative detection and quantitative determination. Application to radiochemistry. *Anal. Chem.* **40**, 586–593 (1968).

63. 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).

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

65. 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).

66. Wendt, I. & Carl, C. The statistical distribution of the mean squared weighted deviation. *Chem. Geol.: Isot. Geosci.* Sect. **86**, 275–285 (1991).

67. Vermeesch, P. High MSWDs are not the problem, low ones are. *Natl. Sci. Rev.* https://doi.org/10.1093/nsr/nwaf036 (2025).

68. Welch, B. L. The generalization of 'STUDENT'S'problem when several different population variances are involved. *Biometrika* **34**, 28–35 (1947).

69. 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).

70. Ickert, R. B. Algorithms for estimating uncertainties in initial radio-genic isotope ratios and model ages. *Chem. Geol.* **340**, 131–138 (2013).

71. 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).

72. Karato, S.-i & Wu, P. Rheology of the upper mantle: a synthesis. *Science* **260**, 771–778 (1993).

73. 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).

74. 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).

75. 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).

76. 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).

77. Hirth, G. & Kohlstedf, D. In: *Inside the Subduction Factory*, Vol. 138 *Geophysical monograph-American Geophysical Union* 83–106 (2003).

78. Korenaga, J. Thermal cracking and the deep hydration of oceanic lithosphere: A key to the generation of plate tectonics? *J. Geophys. Res. Atmos.* **112** https://doi.org/10.1029/2006JB004502 (2007).

79. 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).

80. Connolly, J. A. D. The geodynamic equation of state: what and how. *Geochem. Geophys. Geosyst.* **10** https://doi.org/10.1029/2009GC002540 (2009).

81. Schmidt, M. W. & Poli, S. Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation. *EPSL* **163**, 361–379 (1998).

82. Iwamori, H. Transportation of H₂O beneath the Japan arcs and its implications for global water circulation. *Chem. Geol.* **239**, 182–198 (2007).

83. Hirschmann, M. M. Water, melting, and the deep Earth H₂O cycle. *Annu. Rev. Earth Planet. Sci.* **34**, 629–653 (2006).

84. Crameri, F. Scientific colour maps. *Zenodo* **10**, 5281 (2018).

85. Crameri, F., Shephard, G. E. & Heron, P. J. The misuse of colour in science communication. *Nat. Commun.* **11**, 5444 (2020).

86. 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).

## 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 management of University Grenoble Alpes provided invaluable assistance in establishing and maintaining the ISTerre Micro-Analytical Platform (IMAP), the core facility that enabled the results of this study. We are grateful to Sappi Forest for granting access to the sampling sites. For visualization of geodynamic results, scientific color maps⁸⁴ are used in this study to prevent visual distortion of the data and exclusion of readers with color-vision deficiencies⁸⁵. 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 were supported by the state assignment of the Vernadsky Institute of Geochemistry and Analytical Chemistry of the Russian Academy of Sciences (GEOHI RAS). The computing time for geodynamic modeling (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 sample 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 sample 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 a 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 a leading role in interpreting and visualizing results of geodynamic models and wrote most of the 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 the 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**

The authors declare no competing interests.

**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-025-59024-6.

**Correspondence** and requests for materials should be addressed to Adrien Vezinet, Aleksandr V. Chugunov, Alexander V. Sobolev or Stephan V. Sobolev.

**Peer review information** *Nature Communications* thanks David Murphy, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

© The Author(s) 2025