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

Domain: earth_science
Published: 2025-04-25
Editorial outcome: major_revision

Abstract
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 (87Sr/86Sr = 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 million years-long periods of massive impulsive subduction induced by mantle plumes.

Introduction
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 surface1,2. 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 crust3.

Several approaches have been used to provide estimates of continental crust volumes during the Hadean and Archean3–7. Studies assessing the degree of mantle depletion4,6,7 are viewed as the most reliable since they account for both crustal extraction and recycling3. Most of these studies rely on the bulk rocks’ Sm-Nd, Lu-Hf, and U-Pb isotopic systems3,4,6. 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)8. 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 recognized4. This uncertainty is also valid for the short-lived 146Sm-142Nd decay system9. Also, most ancient rocks are altered, which has been shown to affect Sm-Nd and Lu-Hf isotopic systems10,11. Further, whatever the volume of continental crust on early Earth, the tectonic regime accounting for it remains debated, with suggestions being stagnant-lid, plutonic-squishy-lid, episodic-lid, mobile-lid, ridge-only, and heat-pipe regimes12,13.

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 OIB8,14,15, 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 parentals) that form the bulk of continents14. This observation makes Nb/U and Ce/Pb ideal proxies to assess both the production of continental crust and recycling of restites through time4.

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 melts16 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 interpretations10,11,14.

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 time16–21. 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 crystallization17. 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 isotopes8,14,15 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
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 Hadean24 (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 systems24. This observation has been interpreted as resulting from elemental fractionation occurring during early Hadean magma ocean solidification25.

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 H2O 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 Cr2O3) 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 TiO2 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 TiO2 contents in the melt, as well as Al2O3 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 K2O/TiO2 (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 K2O/TiO2 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 melts26,27 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 between Cl/Ti and Rb/Sr (Fig. 2e), suggesting a gain of Rb as well as Na, K, Pb, U, H2O, 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.

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(χ2) 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 Cl/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), which is close to the average 95% confidence interval (±0.0010) analytical uncertainty of 87Sr/86Sr measurements in individual melt inclusions.

A statistical analysis of our Sr-isotope dataset reveals three statistically homogenous groups based on their initial 87Sr/86Sr (see Fig. 3b and Methods for details). The means of initial 87Sr/86Sr 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 (87Sr/86Sr) initial signature yielding a (87Sr/86Sr) initial weighted mean of 0.69932 ± 0.00024 (95% confidence interval, MSWD = 0.88, p(χ2) = 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 µ142Nd 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(χ2) = 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 km28 in a single mantle plume are focused in local magmatic systems of Mauna Loa and Kīlauea volcanos, which coexist in time and differ in compositions of Sr, Nd, Pb isotopes, and trace elements29. 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 system17,30. 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 142Nd/144Nd 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 extraction4. 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.48,31–33. 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.28,31–33. 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 komatiites25 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-perovskite34. In addition, similar and even much more significant decoupling of εNd and εHf is common in highly depleted abyssal and ophiolitic peridotites because of “following radiogenic ingrowth as a result of preferential fractionation of Lu/Hf during partial melting”35. The main mantle source component of studied komatiite was highly depleted mantle harzburgite, for instance, similar to refractory harzburgites from Archean cratons36, as follows from their contents of incompatible elements (including Sm-Nd and Lu-Hf), which are lower in parental melt than in BSE mantle16. 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 furthers substantiates independent observations obtained via thermodynamic modeling37.

Extraction of continental crust from mafic crust and/or peridotite in the presence of H2O followed by recycling of restite back to the mantle lead 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 melts14. 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 87Sr 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 studies4 using our Nb/U values and Sr geochronological constraints together with the estimated Nb and U contents of the mantle8,31–33 and continental crust38–41. 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 × 1022 kg as the mass of the present-day continental crust and 4.01 × 1024 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 history3,6,42.

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 zircons43,44, the radiogenic isotope compositions of bulk rocks45 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.

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 87Sr/86Sr, 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 StagYY46, 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 domain47, our models generate both basaltic and felsic melts48, include cooling of the core, and use pressure- and temperature-dependent water solubility maps for different mantle minerals12. For this study, we improved previous models12,48 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 studies49,50, 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) composition31 at 4.5 Ga with a low (0.01 wt.%) initial water content12, 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.151, 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 strenth 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 domain52. A mobility value above 1 means that the lithosphere moves horizontally relative to the underlying mantle with significant velocities, and subduction is active13.

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 regime13. 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 87Sr/86Sr 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, 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 (Supplementary 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.

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 subduction53, 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 ago42 as well as several recent models3,6.

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 suite54,55.

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 described16. The new samples 2216, 2217, 2218, and 2222 from Saw Mill and 1528 from Pioneer Complexes are similar to those in composition. These cumulates consist of partially altered (serpentinine 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-H2 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 protocol18. 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, H2O, 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 software56. Details on this approach are reported in earlier studies16. 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 elsewhere55. The Mg-Fe olivine-melt geothermometer was used for the estimation of komatiitic melt inclusion crystallization/entrapment temperature, with correction for the effects of measured water contents57.

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 published55. Average and 2SE uncertainties are reported in Supplementary Data 2–4.

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

Sr-isotope compositions were measured on a Neptune XT MC–ICP-MS (Thermo-Fisher Scientific) following conventional cup configurations59. 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 1013 Ω 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 1011 Ω amplifiers. The different dynamic response time between 1011 Ω amplifiers and 1013 Ω amplifiers was accounted by using the Thermo-Fisher software built-in tau correction60. 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-461. 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 run62), no REE2+ 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 (βSr) was then calculated for each individual analysis based on the measured (88/86) value and then 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 87Rb is then calculated using the measured 85Rb and the natural 87Rb/85Rb value (0.38571). Then, the amount of 87Sr is obtained by removing the contribution of 87Rb on the measured 87 signals (87 net signal = 87Rb + 87Sr). 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 2SD) 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 87Sr/86Sr values. Standardization was, however, required for 87Rb/86Sr 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 87Sr/86Sr and 87Rb/86Sr 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 REE2+ 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 analyses63 and consists of calculating the degree of homogeneity of the reference glasses over all sequence: if the MSWD value—or corresponding p(χ2) 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(χ2) value, indicating statistical homogeneity. No inter-sequence excess variance was required for the 87Sr/86Sr values: MSWD and p(χ2) 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 87Rb/86Sr, 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 87Sr/86Sr, 87Rb/86Sr, 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 confirm that reported uncertainties of individual inclusions are not underestimated.

All melt inclusions discussed in this work have 87Rb/86Sr values measured lower than the 87Rb/86Sr of KL2-G (0.0699, or 0.0248 as Rb/Sr), the reference glass with the most elevated 87Rb/86Sr. Melt inclusions displaying 87Rb/86Sr higher than KL2-G were disregarded since the accuracy of the 87Rb correction could not be assessed. In addition, inclusions with 2 standard errors of measured 87Sr/86Sr 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 87Sr/86Sr 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:
σ inclusion = σ RM × (Sr inc / Sr RM) × (t RM / t inclusion) (1)

Where σ inclusion and σ RM are 95% c.i. uncertainties of individual inclusions and reference material (glasses), respectively; Sr inclusion and Sr RM are Sr contents in inclusions and reference materials; t inclusion and t RM are 87Sr/86Sr acquisition times for inclusions and reference materials (Supplementary Table 1).
Predicted average σ inclusion calculated for all Weltevreden melt inclusions and for the 14 unradiogenic ones are 0.0009 ± 0.0001, as well as their average measured σ inclusion, are 0.0010 ± 0.0001, thus virtually identical (Supplementary Table 1). We, therefore, conclude that reported uncertainties in 87Sr/86Sr 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−2. Carrier gas was He (ca. 0.3–0.9 L/min) with the addition of N2 (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 44Ca2+/44Ca+ 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-G64 and NIST-614-G65) 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 87Sr/86Sr obtained on Weltevreden melt inclusions show a range of values that correspond to an MSWD of 2.3. Considering the number of individual data points—n = 137—this MSWD value translates to a p(χ2) of 4.4 × 10−16, i.e., largely below the accepted threshold of 5% typically used in geoscience to indicate homogeneity of a sample population66. Elevated MSWD ratios are common in Earth science geochemistry, and information from ‘dispersed’ datasets has geological meaning67. This significant heterogeneity is also clearly observed for (87Sr/86Sr) initial values that show an MSWD of 2.5, translating to a p(χ2) of 0. Further, it is also clear from Supplementary Fig. 6 that the range of (87Sr/86Sr) initial values is not resulting from under/over correction of the (87Sr/86Sr) measured values: the range of 87Rb/86Sr is rather narrow and all analyses with elevated 87Rb/86Sr, i.e. > (87Rb/86Sr) KL2-G, were rejected from the dataset, preventing inaccurate over correction of (87Sr/86Sr) initial. Further, since all potential sources of uncertainty have been examined and propagated in isotopic ratios, this range of (87Sr/86Sr) 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 observations17,19,30 for different radiogenic isotope systems, regardless of the analytical protocol (laser ablation, ion probe, or thermal ionization mass spectrometry). 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 (87Sr/86Sr) initial significantly, i.e., considering individual uncertainties (±2se) ‒ 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 (87Sr/86Sr) initial = 0.69932 ± 0.00024 (c.i. 95%) comprised of 14 inclusions (MSWD = 0.88, p(χ2) = 0.57), (ii) one with intermediate (87Sr/86Sr) initial = 0.700624 ± 0.000086 made of 116 inclusions (MSWD = 1.1, p(χ2) = 0.27) and one with elevated (87Sr/86Sr) initial = 0.70177 ± 0.00024 made of 7 inclusions (MSWD = 0.48, p(χ2) = 0.82).

Applying the t-test for independent populations with unequal variances68 reveals that the means of (87Sr/86Sr) 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) (87Sr/86Sr) initial and (ii) BSE evolution. The latter is defined as follows: initial 87Sr/86Sr = 0.698990 at 4.567 Ma age69, initial BSE Rb/Sr = 0.0331, and 87Rb decay constant of 1.3972 × 10−11 yr−1. 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 results are also robust considering the measured Rb/Sr ratios of trapped melts as a 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 spreadsheet70.

Contamination of Weltevreden komatiites
Studied inclusions in olivine yield significant intercorrelated variations of Cl, K, Rb, H2O, 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, H2O, 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 brines26,27. 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 H2O contents. Also note that the extent of contamination tends to increase with decreasing host olivine Fo content (Fig. 2d). The initial 87Sr/86Sr 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 87Sr/86Sr values.

Methods for geodynamic modeling
Essentials and setup. We study the thermochemical evolution of compressible mantle using the code StagYY46, which has been developed and widely used over several decades for global-scale modeling of Earth’s evolution spanning its age12,13,48–50. 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 functions12, the code forms basaltic and felsic melts and considers both intrusive (plutonic) and eruptive (volcanic) magmatism48. 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 rocks48,71. We employ temperature- and pressure-dependent water solubility maps for different mantle minerals, which control the water in-gassing and out-gassing12. We use a two-dimensional quadrant spherical annulus geometry47 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 Ei and activation volume Vi72,73. Following the Arrhenius law, the diffusion and dislocation creep (proxy) components of the temperature- and pressure-dependent viscosity are given as:
η(T,P) = η0 Δηi exp((Ei + P Vi(P)) / RT) (2)

where η0 is the reference viscosity (2:35 1019 Pas) at zero pressure and reference temperature T0 (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) = Vi exp(-P / Pi) (3)

where Pi is the pressure scale, which is different for each layer i. Plastic deformation in the lithosphere is accommodated by assuming plastic yielding52,74. The maximum stress that a material can sustain before deforming plastically is given by the yield stress σY, which has both brittle and ductile components:
σY = min(σY,ductile, σY,brittle) (4)

The ductile yield stress σY,ductile increases linearly with pressure as:
σY,ductile = σ0Y + σ'Y P (5)

where σ0Y 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:
σY,brittle = c + μP (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ϵ̇, where ϵ̇ 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 data75,76, 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−15 1/s and using dry olivine rheology parameters77 (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 sediments51. A similar effective friction coefficient was suggested based on the model of thermal cracking of oceanic lithosphere on early Earth78.

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 αw parameter (0 < α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
ρω = ρd (1 - a CH2O) (7)

where ρω, ρd are densities of water-containing and dry material, respectively; CH2O is the concentration of water in weight fraction, and a is a constant, typically between 1 and 2 estimated using Gibbs energy minimization code79 (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 studies12,48.

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 previously49,50, 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 melts 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 appear 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 (87Sr/86Sr) 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 model47, 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 models6, 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_X80 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.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This manuscript aims to elucidate the formation process of the continental crust using the interdisciplinary approach of\ngeochemical analysis and geodynamic modelling. The geochemical analysis focused on the abundance ratios of Nb/U and\nCe/Pb, and pointed out that 80% of the present continental crust was already formed about 4 billion years ago, and that the\ncontinental crust will be formed in 3.5 billion years in the same amount as at present. In addition to this, they also found the\npresence of unradiogenic Sr in the source-mantle. Geodynamics modelling with TTG crust formation was performed to\nprovide a direct explanation for the continental growth rates suggested by Nb/U and Ce/Pb and Sr isotopes from\ngeochemical analysis, which also attempted to find which types of mode of mantle convection would be relevant to explain\nthose geochemical analysis. The results in geodynamics modelling suggest that the mode of oscillating mobile lid (short-\nterm episodic lid) as a tectonic style of the Earth's surface controls the Hadean geodynamics.\nThe geochemical analysis has been examined in detail from various aspects of the data, and the confirmation of the\npossibility that the continental crust was already generated up to 80% of the present continent 4 billion years ago and grew\nto as much as the present continent 3.5 billion years ago from the results shown in Figure 5 is a significant advance for\ngeology. While this is a significant advance for geology, there are still considerable problems in using geodynamic\nmodelling to directly explain the geochemical analysis data.\nThe reason is that authors did not seem to discuss the model limitations in geodynamic modelling applied in this manuscript.\nFor instance, it is a long-term issue that the geochemical data extracted from geodynamic modelling could not explain the\nactual residence time of isotopes provided in experimental analyses (e.g., Xie and Tackley, 2004 in JGR). In this manuscript,\nauthors attempted using a criterion written in Line 730 to 754, which seems to be a bit odd to me because the experimental\nanalyses should be done with the actual rocks popped out from the deep interior but, by using geodynamics modelling,\nauthors extracted the chemical composition data from a certain depth range (400-600 km depth). I would suspect whether\nauthors could explain the residence time of Rb-Sr isotope system or not (without using a criterion quoted in Line 730 to 754\nin this manuscript).\nThe other aspect which authors should check is how TTG crust could be contributed to formation of the continental crust. In\nFigure 6, as far as I looked at, the continental crust would not be appeared as a result of accumulation of TTG crust. Authors\nshould discuss the difference between TTG crust and actual continental crust, in order to make sure if the formation rate of\nTTG crust is applicable to use the continental growth curve obtained in Figure 5 in this manuscript.\nTherefore, again, it seems reasonable for authors to clearly state the limits of application of the geodynamic model and to\nuse it to validate the interpretations obtained by experimental analysis (in this case, the interpretation on the continental\ngrowth curves in Fig. 5), with a clear statement of the model limitations before applying for geochemical data directly. In this\nsituation, this manuscript is not acceptable for publication in ‘Nature communication’.\nHere below is described for major issue in geodynamics modelling.\nThe model itself includes very complicated processes in mantle dynamics, notably the formation of the TTG crust, and it is\nhoped that it can adequately explain the results of geochemical analyses, but even in this situation a complete\nreconstruction of the thermal and material conditions in the Earth's deep interior is not possible. In order to achieve a\ncomplete reproduction, various observational and analytical data must be explained without over- or under-explaining the\nmodel assumptions. In current framework in geodynamics modelling, several assumptions are made during model\nconstruction. Hence, we could not fully explain geophysical and geochemical observations in modelling approach.\nHowever, it is possible to interpret the results of observations and analyses by clarifying the model application limits.\nHowever, the following issues need to be checked in detail\n1. The manuscript states that Perple_X was used to generate the water solubility limit in the mantle, but it is questionable\nwhether Perple_X can handle the water-saturated phases of the major mantle minerals. The fact that the water content of the\noceanic crust component exceeds 10% in the higher pressurehan 5 GPa is questionable when compared with Iwamori\n(2007 in Chemical Geology), who analysed detailed water solubility limit for the mantle peridotite and the basalt. There is\nalso a point in the mantle peridotite component where the water solubility limit becomes very low around 5 GPa, and after\nthis point, the mantle rocks can only transport the NAMs' worth of water. This point calls into question the estimate of\nwater content allowance using Perple_X.\n2. Regarding the Rb-Sr isotope system, the author describes in the manuscript that the Rb-Sr isotope system was introduced\ninto the simulation model. With regard to the introduced Rb-Sr isotope system, it is questionable whether the simulation case\nthat the authors refer to as the best-fit model is capable of reproducing Figures S6 and S7 to some extent. As pointed out by\nXie and Tackley (2004 in JGR), it seemed to be difficult to explain the trends in U-Pb and Sm-Nd isotopic ratios directly in a\nsimulation model in terms of residence time of isotopes. This means that the validity of the isotopic system model needs to\nbe verified by generating Figures S6 and S7 from the simulation model and comparing them with the results of geochemical\nanalyses. However, the following sampling criteria issues should also be addressed at the same time.\n3. For the geochemical analysis by numerical simulation (Figs 7 and 8), the sampling criteria described in lines 730-754\nneed to be verified whether they are geochemically correct or not. The samples used for geochemical analysis are rocks\nerupted at the surface. However, the sampling in the model reads as picking up data from a certain temperature range\n(condition (i)), depth range (condition (ii)) and Sr content (condition (iii)). It should be verified whether the results obtained\nfrom this sampling approach can be compared with the results of the geochemical analysis. There are no problems with\nconditions (i) and (iii), but I suspect the relevance on the depth range (taken from the mantle transition zone). A justification\non such a criterion would be needed.\n4. As described above, if the aim is to try to explain the experimental data directly, it is necessary to clearly discuss the\napplication limits of the model, taking points 2 and 3 into account. It is necessary to be a little more careful and to clearly\nstate what is explained and what is not as well as the model limitations. Next, it is clear that the TTG crust is formed, but it is\nnot formed as continental crust. If continental crust is formed, the dynamics should change a little more. Repeatedly, from this\npoint of view, the model limitations in authors’ geodynamics modelling should be clarified.\n5. Figure 6: Representative examples of Group 1 and Group 2 should be shown in the case using Figure 7 and Figure 8. In\naddition, the composition field in Fig. 6 could not be seen where the TTG crust is located, which is quite difficult to check due\nto the very small figure. This means that it is unclear whether the time series data on TTG production is correct or not.\nMoreover, although continental crust might be formed from TTG crust, a justification regarding the representation of\ncontinental crust production by TTG formation volumes may be necessary, which might ensure an applicability for\ncomparing between Fig. 5 and Fig 6B and F.\nDetails:\nLine 331: Oscillating mobile lid seems that this regime might have a certain frequency in the surface mobility but, looking at\nFig. 6A, the mobility data never has a certain frequency. Hence, it would be good to call ‘short-term episodic lid’.\nlines 365 and 370 – Information on the water content in the lithosphere would not be needed because authors never show\nthe water content profile as a function of time even in supplemental material. I recommend dropping such statements from\nthe revised manuscript.", "reviewer_2": "The manuscript has greatly improved since the first submission. Overall, I agree with the broad conclusions. The Hadean\nderived Sr isotope group is now very well justified and I agree with the interpretation that it represent a vestige of Hadean\ndepleted mantle. I agree that by extension that there should also be a complementary enriched continent-like crust in the\nHadean. While the method for deciphering this is novel, it should be pointed out that a similar conclusion can be made from\n142Nd variation in early Archean rocks.\nAs the author indicate, the variation in Sr isotopes from melt inclusions in olivine from komatiite flows implies the presence of\nanticrysts and xenocrysts that were transported from depth to erupt at the Earth’s surface. One thing that is majorly missing is\na discussion on how a komatiite magma transports anticrysts and/or xenocrysts to the surface. The general consensus at\npresent is that komatiites are derived from deep melting >2GPa and that the transect the lithosphere without stalling and\nfractionating. As they undergo decompression they cross out of the stability field for olivine and thus any olivine crystals\ncargo would react with the melt and be consumed. I took the Weltevreden komatiite chilled margin composition from\nConnolly et al 2011 and ran a pseudosection for it in MAGEMin-Thermocalc to see the change in the liquidus temperature\nwith pressure. At 2GPa the liquidus is above 1700° but at surface pressure it is at 1550°. For olivine anticrysts and/or\nxenocrysts to survive in a liquid so far about its liquidus temperature requires some explanation.\nLine 114: In the section on the Weltevreden komatiites there needs to be an introduction to the decoupled Nd and Hf\nisotopes in these samples. This is introduced in line 222 without context.\nLines 228 to 231: The recent work of Kamber and Tomlinson on deep melting beneath cratons should be referenced here.\nThe relevance of their observation of the residual mineralogy during the petrogenesis of komatiite melts would be a valuable\ncontribution here.\nFigure 6: Can you include colour bars in A,B,E & F for the Eo- to Mesoarchean and also in figure 7&8?\nFigure 7: I am confused as to what is being presented in panels a and d. In the figure caption in states the\n“Hadean fraction (with 87Sr/86Sr < 0.6997”. However, in the methodology it states that the Hadean fraction was selected for cells with an\ninitial 87Sr/86Sr < 0.6997. Please clarify what is being presented in the figure caption.\nIn addition, it is not clear to me why you are choosing this as a criterion. Because of the initial 87Sr/86Sr in the model of\n0.698990 initially most pixels will pass the above criteria and there will be a gradually decreasing number of pixels that pass\nthe criterion with time. Therefore, of course there would appear to be no effect of continental crust extraction (Line 356) until\nlate in the model. The pixels with 87Sr/86Sr < 0.6997 will include both depleted (most pixels) and undepleted pixels early in\nthe model and thus the ratios will be affected and later in the model when the proportion of pixels 87Sr/86Sr < 0.6997 are\nmuch smaller they will be dominated by high Ce/Pb & Nb/U. The change in the proportion of unradiogenic Sr pixels with\ntime can be inferred from panels c & f.\nI think it would be better to have mode boxes of proportions of high Ce/Pb & Nb/U pixels with a patterning for pixels with\n87Sr/86Sr < 0.6997. then you could visually relate the proportion of depleted pixels to the proportions of TTG in panels b and\ne.\nLines 360: The interpretation presented here is based on the confusing panels a & d in figure 7 and therefore, I am not\nconvinced it is correct.\nLine348: The impact of water on the lithosphere is mentioned here and in the final paragraph of the body\ntext. Yet, the data on the distribution of water in the lithosphere and how it changes with time is not presented in the body text or the\nsupplementary material."}, "author_response": "Response to the reviewers' comments to manuscript NCOMMS-24-70546-T.\nThe reviewers' comments are presented in blue italics, while our responses are in normal black font.\nReviewer #1 (Remarks to the Author):\nThe geochemical analysis has been examined in detail from various aspects of the data, and the\nconfirmation of the possibility that the continental crust was already generated up to 80% of the present\ncontinent 4 billion years ago and grew to as much as the present continent 3.5 billion years ago from the\nresults shown in Figure 5 is a significant advance for geology. While this is a significant advance for\ngeology, there are still considerable problems in using geodynamic modelling to directly explain the\ngeochemical analysis data. The reason is that authors did not seem to discuss the model limitations in\ngeodynamic modelling applied in this manuscript. For instance, it is a long-term issue that the\ngeochemical data extracted from geodynamic modelling could not explain the actual residence time of\nisotopes provided in experimental analyses (e.g., Xie and Tackley, 2004 in JGR).\nWe introduced the Rb-Sr isotopic system in our models, which was not done before in other geodynamic\nmodels, including Xie and Tackley (2004). All our models reproduce very well the Sr isotope composition\nof melt inclusions at the age of emplacement of Weltevreden komatiites (see detailed explanation in the\nanswer to point 2 below). We included new Figure 9, where we show the evolution of the simulated Sr\nisotope ratio in plume material with time and compare modeling results with geochemical data.\nIn current framework in geodynamics modelling, several assumptions are made during model\nconstruction. Hence, we could not fully explain geophysical and geochemical observations in modelling\napproach. However, it is possible to interpret the results of observations and analyses by clarifying the\nmodel application limits. However, the following issues need to be checked in detail\n1. The manuscript states that Perple_X was used to generate the water solubility limit in the mantle, but it\nis questionable whether Perple_X can handle the water-saturated phases of the major mantle minerals.\nThe fact that the water content of the oceanic crust component exceeds 10% in the higher pressure than\nwater solubility limit for the mantle peridotite and the basalt. There is also a point in the mantle peridotite\ncomponent where the water solubility limit becomes very low around 5 GPa, and after this point, the\nmantle rocks can only transport the NAMs' worth of water. This point calls into question the estimate of\nwater content allowance using Perple_X.\nPerple X with the updated database for deep mantle phases is the best available model to simulate\nproperties of materials in the P-T range of the entire mantle. We agree that this model is not perfect and\ndiscuss it in the revised Supplementary Material in the new sub-section “Model limitations” (lines 85-92).\nNote, however, that while maximum saturation water content in oceanic crust phases in our Perple_X-\nbased phase diagram is 8%, in the models presented in the paper, water content at the pressure above 5\nGPa never exceeds 6% and rarely exceeds 4% because of the limited input of the water from the surface.\n2. Regarding the Rb-Sr isotope system, the author describes in the manuscript that the Rb-Sr isotope\nsystem was introduced into the simulation model. With regard to the introduced Rb-Sr isotope system, it\nis questionable whether the simulation case that the authors refer to as the best-fit model is capable of\nreproducing Figures S6 and S7 to some extent. As pointed out by Xie and Tackley (2004 in JGR), it\nseemed to be difficult to explain the trends in U-Pb and Sm-Nd isotopic ratios directly in a simulation\nmodel in terms of residence time of isotopes. This means that the validity of the isotopic system model\nneeds to be verified by generating Figures S6 and S7 from the simulation model and comparing them\nwith the results of geochemical analyses. However, the following sampling criteria issues should also be\naddressed at the same time.\nTackley (2004) or other geodynamic models. All our models accurately replicate the mean strontium\nisotope composition of melt inclusions at the time of emplacement of Weltevreden komatiites, along with\nits standard deviation (refer to the new Figure 9). We assume that melts possess identical Sr isotope\nratios to their respective melt sources. The criteria for selecting numerical cells are conditions (i) and (ii).\nNote that a good fit of the simulated Sr isotope modal ages to the geochemical data provides an\nadditional argument for the validity of condition (ii) for selection of the numerical cells to simulate sources\nof Weltevreden komatiites. Figure S6 illustrates the computation of 87Sr/86Sr ratios of melt inclusions at\nthe time of emplacement (3.27 Ga) derived from measured values utilizing the decay of 87Rb to 87Sr.\nFigure S7 is a qualitative illustration depicting the evolution of the 87Sr/86Sr ratio of potential mantle\nsource components, informed by their Rb/Sr ratio and the decay of 87Rb to 87Sr. Our models inherently\nreplicate all the consequences of the nuclear decay law and thus automatically replicate Figures S6 and\nS7. In response to this comment, we have included in the manuscript a new Figure 9 and explanation on\nlines 422-428 of the revised MS.\n3. For the geochemical analysis by numerical simulation (Figs 7 and 8), the sampling criteria described in\nlines 730-754 need to be verified whether they are geochemically correct or not. The samples used for\ngeochemical analysis are rocks erupted at the surface. However, the sampling in the model reads as\npicking up data from a certain temperature range (condition (i)), depth range (condition (ii)) and Sr content\n(condition (iii)). It should be verified whether the results obtained from this sampling approach can be\ncompared with the results of the geochemical analysis. There are no problems with conditions (i) and (iii),\nbut I suspect the relevance on the depth range (taken from the mantle transition zone). A justification on\nsuch a criterion would be needed.\nIn order to compare models with the geochemical data, we model isotopes and trace elements in the\npotential komatiite source of Hadean age directly in the mantle, rather than the compositions of the melts\nintruded into the crust. The reason is that our (and any other) geodynamic model currently cannot\nreproduce the existence of primary melt from different mantle sources at the crustal level without mixing\nthem, which is achieved using the geochemical approach presented in this paper. However, we can\ndetect a potential Hadean melt source considering numerical cells in the mantle with high enough\npotential temperature (condition “i”) and Hadean Sr modal age (condition “iii”). The extraction of melt\nalters isotope and trace element compositions in heterogeneous sources due to varying fusion of\nheterogeneities; therefore, we must account for numerical cells at the depth prior to melt extraction.\nTherefore, we consider cells at the depth below 400 km (condition “ii”) before this hot plume material is\nascending to the upper mantle and komatiitic melt is extracted. We extended the explanation of modeling\nsampling criteria in the revised Supplementary Material on lines 30-41.\n4. As described above, if the aim is to try to explain the experimental data directly, it is necessary to\nclearly discuss the application limits of the model, taking points 2 and 3 into account.\nDone; see above.\nIt is necessary to be a little more careful and to clearly state what is explained and what is not as well as\nthe model limitations. Next, it is clear that the TTG crust is formed, but it is not formed as continental\ncrust. If continental crust is formed, the dynamics should change a little more. Repeatedly, from this\npoint of view, the model limitations in authors’ geodynamics modelling should be clarified.\nDone in the special sub-section “Model limitations” in the revised Supplementary Material on lines 65-98.\n5. Figure 6: Representative examples of Group 1 and Group 2 should be shown in the case using Figure\n7 and Figure 8.\nin Figures 7 and 8.\nIn addition, the composition field in Fig. 6 could not be seen where the TTG crust is located, which is quite\ndifficult to check due to the very small figure. This means that it is unclear whether the time series data on\nTTG production is correct or not. Moreover, although continental crust might be formed from TTG crust,\na justification regarding the representation of continental crust production by TTG formation volumes may\nbe necessary, which might ensure an applicability for comparing between Fig. 5 and Fig 6B and F.\nThe TTG content in the crust is not seen in the composition field of Fig. 6 because most of the generated\nTTG in our models is recycled back to the mantle. Although a similar process is suggested in\ngeochemical models (e.g., Rosas and Korenaga, 2018), it may be well that our models overestimate TTG\nrecycling rate (discussed in the new sub-section “Model limitations” in the revised Supplementary Material\non lines 66-83. Note that despite the possible overestimation of the TTG recycling rate, the amount of the\nrecycled TTG in our models is typically less than 0.1% in the Hadean fraction of the source (Sr ratio <\n0.6997), and therefore it does not significantly affect the Nb/U and Ce/Pb ratios.\nLine 331: Oscillating mobile lid seems that this regime might have a certain frequency in the surface\nmobility but, looking at Fig. 6A, the mobility data never has a certain frequency.\nWe agree; thanks for pointing that out.\nHence, it would be good to call ‘short-term episodic lid’.\nWe call it now “fluctuating mobile lid regime” because key feature of this regime is its rather high mobility\nlines 365 and 370 – Information on the water content in the lithosphere would not be needed because\nauthors never show the water content profile as a function of time even in supplemental material.\nI recommend dropping such statements from the revised manuscript.\nGood point. We added Figure S11 in the revised Supplementary Material, showing the distribution of\nwater in the same models as in Figure 6.\nReviewer #2 (Remarks to the Author):\nThe manuscript has greatly improved since the first submission. Overall, I agree with the broad\nconclusions. The Hadean derived Sr isotope group is now very well justified and I agree with the\ninterpretation that it represent a vestige of Hadean depleted mantle. I agree that by extension that there\nshould also be a complementary enriched continent-like crust in the Hadean.\nWe express our gratitude for this remark.\nWhile the method for deciphering this is novel, it should be pointed out that a similar conclusion can be\nmade from 142Nd variation in early Archean rocks.\nWe concur that the variation of 142Nd signifies mantle differentiation during the Hadean eon. Notably, the\nkomatiite flow, which yielded the majority of inclusions from our Hadean group (8 out of 14), exhibits the\ngreatest deviation in µ142Nd from terrestrial values (Figure S2). However, µ142Nd alone does not reliably\ndifferentiate between the formation of oceanic and continental crust Our ability to distinguish oceanic vs.\ncontinental contributions in the Hadean component constitutes our principal finding. We include this\ndiscussion on lines 60 and 222-227.\nimplies the presence of anticrysts and xenocrysts that were transported from depth to erupt at the Earth’s\nsurface. One thing that is majorly missing is a discussion on how a komatiite magma transports anticrysts\nand/or xenocrysts to the surface. The general consensus at present is that komatiites are derived from\ndeep melting >2GPa and that the transect the lithosphere without stalling and fractionating. As they\nundergo decompression they cross out of the stability field for olivine and thus any olivine crystals\ncargo would react with the melt and be consumed. I took the Weltevreden komatiite chilled margin composition\nfrom Connolly et al 2011 and ran a pseudosection for it in MAGEMin-Thermocalc to see the change in the\nliquidus temperature with pressure. At 2GPa the liquidus is above 1700° but at surface pressure it is at\n1550°. For olivine anticrysts and/or xenocrysts to survive in a liquid so far about its liquidus temperature\nrequires some explanation.\nThis point of Reviewer 2 is based on misunderstanding. In the paper, we do not suggest that antecrysts\nwere transported from depth. We just state:\n“This variability suggests that the studied olivine crystals are, in fact, antecrysts that grew and\ntrapped their melt inclusions in different places of a complex magmatic plumbing system.”\nFurthermore, we are confident that olivines crystallize and entrap inclusions in the crust under low\npressure, thereby negating the necessity for decompression melting. The low crystallization pressure of\nolivine is ascertained due to the contamination of all melt inclusions in olivine with seawater-altered high-\nsalinity brine and serpentinite prior to their entrapment. This indicates that antecrysts formed within\nseawater-altered shallow oceanic crust. We add this clarification to the text: lines 156-160 and 174-176.\nThe effect of such contamination on the 87Sr/86Sr ratio of melt inclusions is within analytical errors, as\nseen from no correlation between the 87Sr/86Sr ratio and contamination proxies. On lines 233-236 of the\nrevised manuscript, we note, “This feature is expected because the difference between the 87Sr/86Sr\nratios of 3.3 Ga old seawater and Weltevreden komatiite melt is only 0.0020 (Figure 2h), which is merely\ntwice the average 95% confidence interval analytical uncertainty of 87Sr/86Sr measurements in melt\ninclusions.”.\nLine 114: In the section on the Weltevreden komatiites there needs to be an introduction to the decoupled\nNd and Hf isotopes in these samples. This is introduced in line 222 without context.\nDone in revised manuscript. Lines 122-125.\nLines 228 to 231: The recent work of Kamber and Tomlinson on deep melting beneath cratons should be\nreferenced here. The relevance of their observation of the residual mineralogy during the petrogenesis of\nkomatiite melts would be a valuable contribution here.\nDone in revised manuscript. Lines 252-253.\nFigure 6: Can you include colour bars in A, B, E, & F for the Eo- to Mesoarchean and also in Figures 7 &\n8?\nGood idea, thanks. Will be added in the final version of the manuscript.\nFigure 7: I am confused as to what is being presented in panels a and d. In the figure caption in states the\n“Hadean fraction (with 87Sr/86Sr < 0.6997”. However, in the methodology it states that the Hadean fraction was selected for cells with an\ninitial 87Sr/86Sr < 0.6997. Please clarify what is being presented in the figure caption.\nThat is a misprint; of course, Hadean fraction is defined as consisting of cells with 87Sr/86Sr < 0.6997 at\nany model time, not initial.\nWe choose this criterion primarily because we want to test whether we can reproduce the composition of\nthe Hadean group of melt inclusions and find out the physical conditions of how these unradiogenic Sr\nvalues could survive over 700 Ma in convecting mantle.\nBecause of the initial 87Sr/86Sr in the model of 0.698990 initially most pixels will pass the above criteria\nand there will be a gradually decreasing number of pixels that pass the criterion with time.\nCorrect. Shown in panels c and f in Figures 7 and 8.\nTherefore, of course there would appear to be no effect of continental crust extraction (Line 356) until\nlate in the model.\nThere will be an effect if a large amount of depleted U and Pb Hadean material mixes into plumes. But\nyou are right; that is not happening in our models.\nThe pixels with 87Sr/86Sr < 0.6997 will include both depleted (most pixels) and undepleted pixels early in\nthe model and thus the ratios will be affected\nMost pixels with 87Sr/86Sr < 0.6997 early in the model have Nb/U and Ce/Pb ratios close to the BSE\nvalues as expected (see above).\nand later in the model when the proportion of pixels 87Sr/86Sr < 0.6997 are much smaller they will be\ndominated by high Ce/Pb & Nb/U.\nNot obligatory. That will not happen if (i) no continental crust is produced in the Hadean or (ii) if restites\nfrom continental crust production in the Hadean are not recycled into the mantle and mixed into plume\nsources. Increasing of Ce/Pb and Nb/U ratios in the Hadean group with time happens in models of Group\nI, but this effect is much less pronounced in models of Group II, i.e., no subduction, no this effect.\nConsider the following end-member cases.\nCase 1: No continental crust is produced in the Hadean. Obviously, in this case, Nb/U and Ce/Pb ratios in\nthe Hadean fraction will have BSE values.\nCase 2. Continental crust is produced in the Hadean, but restites are not recycled and remain in the\nlithosphere. In this case again, Nb/U and Ce/Pb ratios in the Hadean fraction sampled by plume will have\nBSE values.\nCase 3. Continental crust is produced in the Hadean, and restites of Hadean age are immediately\nrecycled into the lower mantle and mixed into the plume sources. In this case, Nb/U and Ce/Pb ratios in\nthe Hadean fraction sampled by the plume will increase during Hadean time proportionally to the\nproduction of the continental crust and will stay constant after Hadean time.\nFinally intermediate Case 4. Continental crust is produced in the Hadean, and restites of Hadean age\nare recycled into the lower mantle and mixed into the plume sources not immediately, but continuously\nduring some time, say 1 Gyr. In this case, Nb/U and Ce/Pb ratios in the Hadean fraction will increase\ncontinuously over time even after Hadean. Our models of Group I follow case 4, and models of Group II\nfollow case 2 with some addition of case 4.\nThis discussion was added in the revised Supplementary Material on lines 43-63.\nThe change in the proportion of unradiogenic Sr pixels with time can be inferred from panels c & f.\nIndeed, we concur, and this is precisely illustrated in panels c and f. Note that this proportion matches\ngeochemical data in the case of models Group I, but it is too small in the case of Group II models.\nfor pixels with 87Sr/86Sr < 0.6997. then you could visually relate the proportion of depleted pixels to the\nproportions of TTG in panels b and e.\nPanels “a” and “d” in Fig. 7 and much less in Fig. 8 do show increasing Nb/U and Ce/Pb ratios with time in\nthe cells of the Hadean age (87Sr/86Sr < 0.6997) because of increasing proportions of TTG production\nrestites in the source.\nLines 360: The interpretation presented here is based on the confusing panels a & d in figure 7 and\ntherefore, I am not convinced it is correct.\nThat is a misunderstanding. See explanation above\nLine 348: The impact of water on the lithosphere is mentioned here and in the final paragraph of the body\ntext. Yet, the data on the distribution of water in the lithosphere and how it changes with time is not\npresented in the body text or the supplementary material.\nValid observation. We add Figure S11 in the revised Supplementary Material. illustrating the distribution\nof water and its temporal evolution for representative models of both groups in the supplementary\nmaterial of the revised manuscript.\n6"}, {"round": 2, "reviewer_comments": {"reviewer_1": "Thank you for addressing my points in the earlier version of submitted manuscript.\nMost points addressed by authors would be sufficient to me as well as convince me.\nHowever, after looking through a revised version, here a couple of points are still raised up, which authors should address a\nbit more carefully.\nRegarding the residence time issue, authors did not use the similar way to Xie and Tackley (2004). In order to avoid such an\nissue, authors assume the sampling criteria quoted in ‘Method section’. As far as I read the supplemental material, the first\nand last paragraph seemed to quote this way but quite unclear. Please clarify a bit more on the difference between previous\ngeodynamics modelling and authors’ way (mainly the last paragraph of model limitations section), including more clear\nstatement on geochemical justification of authors’ criteria as well as numerical reasoning quoted in the manuscript.\nOn the solubility limit of water, please cite Iwamori (2007 in Chemical geology) and discuss a potential difference between\nauthors’ solubility and Iwamori’s one. Please expand the discussion done in l.84-l,92 in the supplemental material as well as\nalign l.416 to l.418 in the main text. This is because the solubility maps between two ways are quite different. The solubility\nmap using this manuscript would lead to overestimate the mantle water content and water flux into the mantle transition zone\nquoted in the revised manuscript, but this would be found for the more unrealistic situation regarding the effect of water in the\nupper mantle. In terms of such a difference from the experimental-based water solubility map, the authors should carefully\ndiscuss about the effect of water in the geodynamics modelling.\nI guess that this manuscript is almost there for a result which authors desire... but more careful discussions would be\nneeded.", "reviewer_2": "I am happy with most the revisions of the manuscript. The response on the origin of olivine anticyrst in a komatiite magma\nremains weak. The unradiogenic Sr isotope data will undoubtable come under scrutiny and so too will the hypothesis of\nolivine anticrysts in a komatiite magma.\nIf the olivine crystals are anticrysts, the way I envisage the origin of the unradiogenic Sr isotope melt inclusions is that they\nformed in the region of partial melting of the mantle where the komatiite magma originated where chemically discrete\nlithologies partially melt and the melt migrates into a melt network, accumulates and homogenises until it reaches sufficient\nvolume for the buoyancy to cause it to intrude through the lithosphere. It is in the early melt network at pressures of >2GPa ,\nprior to accumulation and homogenisation, that I envisage a scenario where an olivine could grow and trap a chemically\ndistinct magma.\nAll other scenarios I can think off are for xenocrysts of olivine, either lithospheric olivine from disaggregated lithospheric\nmantle or olivine in crustal lithologies that were assimilated by the komatiite magma.\nI appreciate that the authors may not be able to address the origin of these anticrysts/xenocrysts in this current paper. But the\nauthors should at least acknowledge that the mechanism for the komatiite flow to capture olivine remains unclear and\nrequires further study.\notherwise, I am happy for this manuscript to be published."}, "author_response": "RESPONSE TO REVIEWER COMMENTS\nReviewer comments are in blue italics; our response is in normal black.\nReviewer #1 (Remarks to the Author):\nThank you for addressing my points in the earlier version of submitted manuscript.\nMost points addressed by authors would be sufficient to me as well as convince me.\nWe are grateful to the reviewer for the prompt and positive response.\nHowever, after looking through a revised version, here a couple of points are still raised up,\nwhich authors should address a bit more carefully.\nRegarding the residence time issue, authors did not use the similar way to Xie and Tackley\n(2004). In order to avoid such an issue, authors assume the sampling criteria quoted in ‘Method\nsection’. As far as I read the supplemental material, the first and last paragraph seemed to quote\nthis way but quite unclear. Please clarify a bit more on the difference between previous\ngeodynamics modelling and authors’ way (mainly the last paragraph of model limitations\nsection), including more clear statement on geochemical justification of authors’ criteria as well\nas numerical reasoning quoted in the manuscript.\nWhile we use the same technique to compute isotope and trace element evolution in the StagYY\ncode as Xie and Tackley (2004) did, our models have several essential differences. First, the\nStagYY code has significantly evolved since 2004, and a key new feature in our models is the\nproduction of felsic melts. Continental crust production drastically affects the evolution of many\nisotopic systems and trace elements. Second, we introduce the Rb-Sr isotope system, which was\nnot included in previous studies. Our modeling aim is also different. We aim to compare the\nmodeling results with new geochemical data reporting the trace element composition of the\nHadean melt source in the Archean komatiites. To achieve this, we analyze the compositions of\npotential sources of komatiitic melts at the depth before these melts are separated from the\nsource rocks and before the melts from different sources are mixed.\nWe added this discussion in the Methods section of the revised manuscript.\nOn the solubility limit of water, please cite Iwamori (2007 in Chemical geology) and discuss a\npotential difference between authors’ solubility and Iwamori’s one. Please expand the discussion\ndone in l.84-l,92 in the supplemental material as well as align l.416 to l.418 in the main text.\nThis is because the solubility maps between two ways are quite different. The solubility map\nusing this manuscript would lead to overestimate the mantle water content and water flux into\nthe mantle transition zone quoted in the revised manuscript, but this would be found for the more\nunrealistic situation regarding the effect of water in the upper mantle. In terms of such a\ndifference from the experimental-based water solubility map, the authors should carefully\ndiscuss about the effect of water in the geodynamics modelling.\nThe phase diagram of water solubility in Iwamori's paper is based on the experimental study by\nSchmidt and Poli in 1998. The PerpleX tables are also derived from experimental data and\nthermodynamic models developed by Holland and Powell, with updates made in 2004. We\nutilize the most recent model.\nSince we are not the experts in mineral physics, we cannot be sure whether water solubility in\nbasalt composition at high pressures in the PerpleX model is overestimated, and we prefer not to\ndiscuss this in the paper. However, even if it is overestimated, the water input in our models is\nvery moderate. The water content in the upper mantle is also quite realistic, at a few hundred\nppm (Fig. S11), which does not contradict Iwamori (2007).\nWe added reference to Iwamori (2007) and a short discussion in the Methods section of the\nrevised manuscript.\nI guess that this manuscript is almost there for a result which authors desire... but more careful\ndiscussions would be needed.\nThanks again. The required discussions are added.\nReviewer #2 (Remarks to the Author):\nI am happy with most the revisions of the manuscript. The response on the origin of olivine\nanticyrst in a komatiite magma remains weak. The unradiogenic Sr isotope data will\nundoubtable come under scrutiny and so too will the hypothesis of olivine anticrysts in a\nkomatiite magma.\nIf the olivine crystals are anticrysts, the way I envisage the origin of the unradiogenic Sr isotope\nmelt inclusions is that they formed in the region of partial melting of the mantle where the\nkomatiite magma originated where chemically discrete lithologies partially melt and the melt\nmigrates into a melt network, accumulates and homogenises until it reaches sufficient volume for\nthe buoyancy to cause it to intrude through the lithosphere. It is in the early melt network at\npressures of >2GPa , prior to accumulation and homogenisation, that I envisage a scenario\nwhere an olivine could grow and trap a chemically distinct magma.\nAll other scenarios I can think off are for xenocrysts of olivine, either lithospheric olivine from\ndisaggregated lithospheric mantle or olivine in crustal lithologies that were assimilated by the\nkomatiite magma.\nI appreciate that the authors may not be able to address the origin of these anticrysts/xenocrysts\nin this current paper. But the authors should at least acknowledge that the mechanism for the\nkomatiite flow to capture olivine remains unclear and requires further study.\notherwise, I am happy for this manuscript to be published.\nWe are grateful to the reviewer for the prompt and positive response. However, we continue to\nbelieve that antecrysts in Weltevreden komatiites crystallized in shallow crust. The strongest\nevidence for this, as we already noted in the paper, is the assimilation of seawater-altered\nmaterials in ALL trapped melts, including the Hadean group (see Figures 2 d, e, f).\nWe include the following text in the revised version of the manuscript:\n“The contamination of all olivine-hosted melt inclusions from Weltevreden komatiites by\nseawater-altered crust (Figures 2 d-f) and low-density fluid bubbles in quenched inclusions\n(Figure 1c) indicates the crystallization of olivine antecrysts within shallow oceanic crust. The\nability of crystallizing olivine to sample isotopically diverse mantle-derived melts in shallow\ncrust implies that melts from deep heterogeneous mantle sources could be delivered and focused\nat crustal levels without complete mixing. Recent examples include the shield volcanoes of the\nHawaiian plume, where melts formed at depths of around 100 km [Herzberg et al., 2023] in a\nsingle mantle plume are focused in local magmatic systems of Mauna Loa and Kilauea volcanos,\nwhich coexist in time and differ in compositions of Sr, Nd, Pb isotopes, and trace elements\n[Pietruszka et al., 2013]. Furthermore, melts with various trace elements and Sr, Pb, and Nd\nisotope compositions are trapped in olivine crystallized in Mauna Loa volcano's shallow\nplumbing system [Sobolev et al., 2011; Reinhard et al., 2018].”\nFinally, the adjacent Weltevreden komatiite flows show significant differences, not only in long-\nlived isotopic systems Lu-Hf and Sm-Nd but also in the 142Nd/144Nd isotopic ratios (see Fig. S2).\nThese differences cannot be attributed to alteration; instead, they suggest that the primary\nmelts delivered to the surface were incompletely mixed. (These marked sentences we will\ninclude in the new revised version of MS."}, {"round": 3, "reviewer_comments": {"reviewer_1": "As I read the revised manuscript and authors' replies, I would not be still convinced for the water solubility issues but this\nwould not be an issue on the main stream on scientific content in the revised manuscript.\nTherefore, it is time that I should endorse the revised version of manuscript for publication.\nThank you for taking your time to address many of my comments.\nI hope that this manuscript would give a high impact to the community."}, "author_response": ""}]

Decision letter