Calcium isotope evidence for early Archaean carbonates and subduction of oceanic crust

Continents are unique to Earth and played a role in coevolution of the atmosphere, hydrosphere, and biosphere. Debate exists, however, regarding continent formation and the onset of subduction-driven plate tectonics. We present Ca isotope and trace-element data from modern and ancient (4.0 to 2.8 Ga) granitoids and phase equilibrium models indicating that Ca isotope fractionations are dominantly controlled by geothermal gradients. The results require gradients of 500–750 °C/GPa, as found in modern (hot) subduction-zones and consistent with the operation of subduction throughout the Archaean. Two granitoids from the Nuvvuagittuq Supracrustal Belt, Canada, however, cannot be explained through magmatic processes. Their isotopic signatures were likely inherited from carbonate sediments. These samples (> 3.8 Ga) predate the oldest known carbonates preserved in the rock record and confirm that carbonate precipitation in Eoarchaean oceans provided an important sink for atmospheric CO2. Our results suggest that subduction-driven plate tectonic processes started prior to ~3.8 Ga.

T he mechanisms and timing of continental crust formation are important to understand because they are intimately tied to (i) the history of plate tectonics, (ii) the chemical evolution of the atmosphere and oceans, and (iii) the proliferation of life on Earth 1-3 . There is still much controversy, however, surrounding the formation of continents and the timing for onset of subduction-driven plate tectonics, which is partly due to the difficulty of interpreting Earth's scarce early rock record. Archean continental crust (>2.5 Ga) accounts for only~5% of Earth's modern surface and is dominated by "granite-greenstone belts" that consist of subordinate metabasalts (~20%) and silicic plutonic rocks known as tonalite-trondhjemite-granodiorite (TTG) suites 4 . TTG suites no longer form in the present day as a consequence of lower mantle temperatures, yet adakites, which are formed during rare melting of subducting hot/young oceanic crust 5 , are the closest modern analogs for Archean TTGs 4,6 .
TTG suites unequivocally represent Earth's earliest preserved continental crust, yet whether or not subduction of oceanic crust is required for their formation is still debated. Other hypotheses suggest that TTGs derive from processes such as melting of hydrated mafic rocks at the base of thick oceanic plateaus 4,7 or extensive fractional crystallization of basaltic magmas in the midto-lower crust 8 . Using Si isotopes, recent studies suggest that surface-derived materials were recycled into the sources of TTGs 9,10 , thus favoring a horizontal tectonic scenario 11 . Yet, trace-element ratios, which are typically used to infer PT conditions for TTG petrogenesis, can lead to ambiguous results due to (i) the competing effects of temperature and pressure, (ii) the potential for open-system behavior 12 , and (iii) the similar traceelement signatures for restitic garnet (high-pressure) and hornblende (low-pressure), leading to significant debates regarding the depths of TTG generation 13 . Thus, it can still be questioned whether surficial materials (e.g., chert) were incorporated into TTG source rocks through subduction or through other processes such as burial during top-down construction of oceanic plateaus 7 .
Here, we develop a stable Ca isotope proxy that can constrain the apparent geothermal gradients (dT/dP) along which TTG magmas were generated and show that these can be used to discriminate between different geodynamic settings for the formation of ancient continental crust. Calcium isotopes are ideal for investigating TTG petrogenesis because they (i) are sensitive to magmatic processes and source variations, (ii) commonly equilibrate in plutonic settings, (iii) have well-defined equilibrium fractionation factors, and (iv) are insensitive to redox effects 14 .
Initial neodymium and hafnium isotopic compositions have limited ranges consistent with juvenile (mantle) sources and do not correlate with δ 44 Ca (Fig. 1b, c), suggesting that assimilation of older crust did not influence isotopic signatures. We find that a majority of TTGs and adakites have δ 44 Ca values that are anticorrelated with trace-element ratios associated with increasing residual garnet in source rocks at higher pressures (Supplementary Fig. 1). Since garnet is predicted to have the highest δ 44 Ca among Ca-rich minerals 14 , this observation suggests that Ca isotopes are sensitive to both temperature and to pressuredependent changes in residual mineralogy.
Phase-equilibrium modeling. In order to best address the opposite isotopic fractionation effects of increasing pressures and temperatures ( Supplementary Fig. 2a), and to account for the complex changes in mineral proportions and compositions in TTG source rocks during progressive melting ( Supplementary  Fig. 2b, c), we have incorporated equilibrium Ca isotope fractionation into phase-equilibrium models for TTG petrogenesis 12 (see "Methods," Supplementary Notes 2 and 3). In Fig. 2a show the evolving δ 44 Ca of silicate melts during progressive heating along three different geothermal gradients for two endmember scenarios (closed-system equilibrium vs. garnet fractionation, Supplementary Data 4). The overall trends are dominated by the competition between increasing temperature (decreasing isotopic fractionation) and increasing pressure (increasing isotopic fractionation due to increasing residual garnet), with the largest effects predicted along the lowest dT/dP path (500°C/GPa). We focus our discussion on models that use a depleted Archean tholeiite (DAT) protolith because they provide the best match to our measured major and trace-element data ( Supplementary Fig. 3) 12 . Broadly similar results, however, are obtained for enriched Archean tholeiite ( Supplementary Fig. 4), suggesting that our models are relatively insensitive to variations in protolith composition.
Geothermal gradient constraints in Archean granitoids. In Fig. 2b, we compare data and model results for δ 44 Ca vs. Dy/Yb (a trace-element proxy for restitic garnet). In concert with δ 44 Ca data and model results for other common trace-element proxies, we find that most TTGs, along with modern and Mesozoic adakites, are best explained by dT/dP paths between 500 and 750°C/GPa. The required geothermal gradients are similar to those recorded in hot "modern" (<750 Ma) subductionzone assemblages 16 but with overall higher temperatures that enable direct melting of mafic (oceanic) crust. Given that our data best agree with closed-system model results, this suggests that TTG melts may have initially encountered thermal or rheological impediments to melt loss 17 . Using the wet-basalt (DAT) solidus 12 as a lower temperature limit (Fig. 2c), our results agree with previous predictions for TTG petrogenesis (800-1000°C at 1-2 GPa) 18 . These results are also similar to PT estimates for modern adakites 6 , which are typically formed through subduction of hot/ young oceanic crust 5 20 ] therefore suggest that Archean TTGs most likely formed through (hot) subduction.
Although the results point toward melting of hydrated oceanic crust as the source of TTG magmas, slightly negative europium anomalies in a number of samples point toward a role for plagioclase fractionation within TTG plutons after emplacement 4,8 . Given that plagioclase is the only mineral that is isotopically lighter than silicate melt at equilibrium 14 , fractional crystallization of plagioclase during plutonic cooling at shallow depths (which would not affect the inherited Dy/Yb signatures) provides the most obvious means of increasing melt δ 44 Ca to values heavier than BSE (see Supplementary Note 5 and Supplementary Fig. 5), as observed in three of our TTG samples. Plagioclase fractionation must have also occurred in many of the other samples, suggesting that measured δ 44 Ca likely represent maximum values for parental TTG melts, prior to crystallization, and that our geothermal gradient predictions generally represent upper estimates (Supplementary Note 5).
Assimilation of carbonate sediments. Of the five samples that lie below our model predictions for 500°C/GPa, two have positive Eu/Eu* anomalies that can be explained through plagioclase accumulation (Supplementary Note 5) and one deviates only slightly from the predictions and could suggest formation along moderately lower geothermal gradients. The two NSB samples (~3.8 Ga, Fig. 1a), however, have unusually low δ 44 Ca that cannot be explained by equilibrium magmatic processes. This implies that these samples record either (i) kinetic isotope fractionations or (ii) incorporation of isotopically distinct Ca-rich materials, not recorded by other samples. In volcanic systems, large kinetic isotope effects can result from Ca diffusion during rapid crystal growth 21 . This mechanism cannot explain the composition of NSB samples, however, because a negative shift in bulk magma δ 44 Ca would require rapid growth and removal of minerals where Ca is strongly incompatible, which (by definition) would have little effect on the Ca budget. The TTG samples also (i) lack compositional banding, (ii) were collected far from lithological contacts, (iii) have similar initial Hf isotopic compositions in zircons and their host rocks 11 , and (iv) have no apparent correlations between δ 44 Ca and fluid mobile/immobile trace-element ratios ( Supplementary Fig. 6), suggesting that significant changes in bulk-rock chemistry associated with post-emplacement metamorphism/metasomatism did not take place. Although they are generally similar to other TTGs in terms of major/trace-element chemistry and Si isotopic compositions ( Supplementary Fig. 7), the NSB samples have much higher peraluminosity [denoted "A/ CNK," defined as molar Al 2 O 3 /(CaO+Na 2 O+K 2 O)] and oxygen isotope ratios (δ 18 O), which are typically attributed to assimilation of metasediments 22 (Fig. 3a). Thus, a diffusion-based fractionation mechanism, though possible, would imply that the high peraluminosity and elevated δ 18 O in NSB samples are merely coincidental. The most parsimonious explanation for low δ 44 Ca in these samples therefore is the incorporation of metasedimentary materials [e.g., carbonated ("marly") sediments] with high A/ CNK, high δ 18 O, and low δ 44 Ca values.
Calcium isotopes have been used as a proxy for recycled carbonates in a number of magmatic systems; however, many previous studies assume low δ 44 (Fig. 3b). Using a three endmember mixing model (average TTG, average shale, and average carbonate with δ 44 Ca = −1.25‰) shown in Fig. 3a, c, we find that NSB samples are best explained through incorporation of 30-50 wt% carbonated metasediments (shale and/or hyaloclastite containing 5-20 wt% carbonates, see "Methods"). Altered oceanic crust could also potentially serve as a source of carbonates in TTG protoliths, but given that these carbonates typically have higher δ 44 Ca values 25 , it is unlikely that they provide enough isotopic leverage to explain the NSB samples (Supplementary Note 6a).
As depicted in Fig. 4, our data provide further evidence for the incorporation of surficial materials (carbonates and shales) into TTG sources, while Si isotope data for the same samples suggest incorporation of chert 9,10 . Thus, several types of sediment [often found together in accretionary wedges (Supplementary Note 7)] could undergo subduction and melting in the early Eoarchean. The ubiquity of Si isotope signatures indicative of chert incorporation in TTGs, however, contrasts with the limited number of samples that require assimilation of carbonates, potentially indicating that carbonates were geographically limited, difficult to incorporate into melts during subduction, and/or generally had δ 44 Ca signatures similar to average TTGs (e.g., −0.3‰, which would represent slow/equilibrium precipitation from BSE-like seawater at~25°C, Supplementary Note 6). Thus, we cannot rule out that TTGs other than the NSB samples also incorporated carbonates. Given that their δ 44 Ca values can be reproduced through magmatic processes, however, carbonates are not necessary to explain the bulk of our data.
Our results constrain the apparent geothermal gradients along which oceanic crust (and metasediments) melted to produce TTG magmas, providing further evidence for the early operation of subduction-driven plate tectonics (Fig. 4) and contrasting with estimates relying on the appearance of blueschists or of paired metamorphic belts in the rock record 26 . Although we are analyzing the products of melting, as opposed to the metamorphic residues, our results suggest that relatively low dT/dP conditions indicative of subduction may be underrepresented in the ancient metamorphic record ( Supplementary Fig. 8), in agreement with previous predictions 27 . The data also provide independent evidence for the presence of carbonated sediments on the ocean floor in the early Eoarchean (>3.8 Ga), predating the oldest carbonate units preserved in the rock record 28 and suggesting that carbonates from earlier periods may be underrepresented. This implies that the silicate-carbonate cycle was in place by~3.8 Ga and provided a necessary sink for high levels of volcanic CO 2 outgassing [which would have otherwise led to a runaway greenhouse effect on Earth 29 ]. The results thus not only provide support for models that require elevated atmospheric pCO 2 to compensate for the faint young sun 30 but also have significant implications for continental emergence/weathering through time 31,32 . Although our data do not definitively prove the existence of a global network of interconnected plate boundaries (as required by some definitions of "plate tectonics"), we show that subduction events occurred repeatedly throughout the Archean, in agreement with a growing body of evidence suggesting that plate-tectonic processes started prior to 3.5 Ga 1,6,9,10,20,31-34 .

Methods
Chemical and isotopic analyses. In order to perform chemical and isotopic analyses relevant to this study, aliquots of powdered rock sample (~20-50 mg) are first dissolved in 2:1 mixtures of concentrated HF + 6 M HCl and refluxed at 130°C for 1 week. The samples are then evaporated to dryness and fully redissolved in 4 M HNO 3 (also refluxed at 130°C) for subsequent Ca isotope and major/traceelement analyses. Sample descriptions are available in Supplementary Note 1 and Supplementary Data 2.
Ca isotope analyses. Calcium (~50 μg) is separated from aliquots of the dissolved samples and standards using established column chemistry methods adapted from UC-Berkeley 21,35-38 and modified for multi-collector inductively coupled plasma mass spectrometry (MC-ICP-MS) analyses 39 . Eichrom ® DGA resin is thoroughly washed by alternating between ultrapure water and 4 M HNO 3 several times overnight before loading into home-made Teflon columns with~2 mL reservoirs. The resin (~250 μL) is loaded into the columns and washed with 2 mL H 2 O, 2 mL 4 M HNO 3 , 2 mL H 2 O, 1 mL 6 M HCl, 2 mL H 2 O, and 2 mL 4 M HNO 3 . Samples are then loaded onto columns in~1 mL of 4 M HNO 3 , and a total of 6.8 mL of 4 M HNO 3 (in progressively larger batches) is used to elute matrix elements. The Ca is then collected using 3.5 mL of ultrapure water (500, 1000, 2000 μL), evaporated to dryness, and then redissolved with one drop of~50% H 2 O 2 and concentrated HNO 3 (allowed to sit overnight to destroy any leached organics or resin that may have passed through the column frits). The samples are evaporated to dryness and redissolved in 1 mL of 4 M HNO 3 to repeat the entire process with fresh resin, ensuring the effective removal of matrix elements with isobaric interferences (e.g., Sr). After the second pass through the columns, the purified Ca from samples and standards is dissolved in 5 mL 0.1 M HNO 3 (to make~10 parts per million (ppm) Ca solutions) for isotopic analysis.
Calcium isotope compositions ( 44 Ca/ 42 Ca and 44 Ca/ 43 Ca) are measured via MC-ICP-MS using sample-standard bracketing on a Thermo Fisher Neptune instrument at IPGP. A majority of our measurements use SRM915b as a bracketing standard due to greater availability. To allow for easier comparison with published data and to best illustrate complementarity between terrestrial reservoirs, δ 44/42 Ca values are converted to δ 44 Ca by multiplying by 2.05 (the kinetic mass law) and are reported relative to BSE [with a recommended value of +0.95‰ relative to SRM915a 14 , which is +0.25‰ relative to SRM915b in this study]. As 40 Ca is not directly measured, the use of MC-ICP-MS for stable Ca isotope analyses precludes the need for radiogenic 40 Ca corrections, which could otherwise lead to significant variations in Archean-aged samples analyzed by thermal-ionization mass spectrometry 35,36,40 .
Samples are introduced into the mass spectrometer using an Apex desolvating nebulizer attached to an autosampler system with a probe aspiration rate of~100 μL/min.  reference materials W2-a (dolerite) and AGV-2 (andesite) are also measured alongside our samples and yield values of +0.79 ± 0.03‰ (2SE, n = 19) and +0.72 ± 0.02‰ (2SE, n = 9), relative to SRM915a, respectively, in agreement with previous studies 36,[41][42][43][44] . This suggests that our analytical methods are robust and suitable for the silicate rocks targeted in this study.
Major and trace-element analyses. Element concentration measurements are conducted using an Agilent 7900 quadrupole inductively coupled plasma mass spectrometer (Q-ICP-MS) at IPGP. Aliquots of~20-100 μL (representing 0.5 mg of rock) are taken from each dissolved sample and diluted prior to analysis in order to produce concentrations within the calibrated range of the instrument. The results are calibrated against multiple in-house standards, with calibration standards intermittently checked every ten samples to ensure accurate standardization of the analyses. Calibration standards and the internal element standards Sc and In are used to monitor and correct for drift and matrix effects. Major element concentrations are reported in wt% and trace-elements are reported in ppm (Supplementary Data 1). Our elemental data for USGS standards AGV-2 and W2-a agree with accepted values from the GeoREM database 45 , and TTG samples agree well with previously reported measurements [46][47][48][49][50] . Reproducibility for all elements is better than 5% (relative standard deviation).
Ca isotope phase-equilibrium modeling. In order to predict the stable Ca isotope composition of melts produced during progressive heating along different geothermal gradients, we incorporate ab initio estimates for equilibrium Ca isotope fractionation into the phase equilibrium models described above. This represents a significant improvement over previous modeling attempts where temperatures, pressures, and/or mineral proportions/compositions are fixed 15,51 and gives a more realistic picture of Ca isotopic fractionation during melt production in evolving magmatic systems. Phase equilibrium modeling methods are described in Supplementary Note 2 and are based on those of ref. 12 . Reduced partition function ratios (RPFRs) used in our Ca isotope models are discussed in Supplementary Note 3 and compiled in Supplementary Table 2.
Mineral-melt Ca isotope fractionation. To account for decreasing isotopic fractionation with increasing temperatures, the RPFR for each phase is recalculated at each temperature (where T is in Kelvin) according to: The Ca isotope fractionation factor between two phases a and b, (α a−b ), at a given temperature can then be calculated as the ratio of their temperature-dependent RPFRs: In the case of determining fractionation between melt and a set of minerals with changing proportions/compositions, the effective RPFR for the bulk solid phase is recalculated at each step by combining the RPFRs for each mineral and weighing them by their relative contributions to the total Ca budget of the solids, where X i Ca (P,T,x) is the fraction of total Ca in phase i (which varies with pressure, temperature, and bulk composition) and comes from the phase equilibrium modeling results.
Closed-system equilibrium calculations. In order to model the isotopic evolution of Ca during closed-system melting of a starting material with BSE δ 44 Ca, we start from the mass-conservation equation for equilibrium trace-element distributions between solids and melt in a closed system: In the closed-system scenario, we assume that δ 44 Ca bulk = 0‰. Comparing the results for Eq. 5 with a commonly used approximation of the form: We find that the results agree within <0.1 ppm (equivalent to <0.0001‰), which is negligible compared to typical analytical uncertainties (~0.1‰), and so we use the simpler equation (Eq. 7) in our models.
Garnet fractionation calculations. In order to simulate the progressive isolation of garnet cores from the rest of the bulk system, we evolve the bulk chemical and isotopic compositions of the system such that 5/6 of the garnet is removed every time 5 mol% garnet is reached [equivalent to~6 wt%, depending on the solid-solution compositions 12 ]. Afterwards, 1/6 of the garnet (e.g., the reactive rims) remains in the system and equilibrates chemically and isotopically with the other mineral phases and melt. The isotopic evolution of the system is thus calculated in an iterative fashion according to: where n is the number of steps (1 mol% of garnet growth, starting at n = 0) in the iterative calculation, δ 44 Ca nongrt (n) and δ 44 Ca grt (n) are the Ca isotopic compositions of all non-garnet phases (including melt) and of garnet, respectively (in the previous step), X i Ca (n) comes from the phase-equilibrium modeling results, and N 1 indicates natural numbers (beginning with 1). Note that removing larger batches of garnet (e.g., 10%) would decrease the δ 44 Ca effects predicted in our models, while smaller batches (e.g., 1%) would increase the predicted effects and approach those that would be predicted by a Rayleigh fractionation model. For simplicity, we only show the dependence of each above variable on n; however, we note that each of these variables also depends on pressure, temperature, and bulk system composition (P,T,x), as in the previous closed-system calculations (Eqs. 1-7). The Ca isotopic compositions of non-garnet (including melt) and garnet phases are calculated: where we assume that δ 44 Ca bulk (n = 0) is equal to BSE (=0‰) and where the Ca fraction of each phase i [X i Ca (n)] comes from the phase-equilibrium modeling results. The fractionation factor between garnet and non-garnet phases is calculated: After calculating the bulk composition of the system at each step of progressive garnet growth and sequestration (Eqs. [8][9][10][11], the δ 44 Ca of the melt can then be calculated in the same fashion as for the closed-system equilibrium case (Eq. 7) but with each variable depending on the evolving δ 44 Ca bulk (n) and bulk chemical composition of the system as garnet is removed (Eq. 12).

Data availability
All data are available in the main text, Supplementary Information,