Statistical geochemistry reveals disruption in secular lithospheric evolution about 2.5 Gyr ago

The Earth has cooled over the past 4.5 billion years (Gyr) as a result of surface heat loss and declining radiogenic heat production. Igneous geochemistry has been used to understand how changing heat flux influenced Archaean geodynamics, but records of systematic geochemical evolution are complicated by heterogeneity of the rock record and uncertainties regarding selection and preservation bias. Here we apply statistical sampling techniques to a geochemical database of about 70,000 samples from the continental igneous rock record to produce a comprehensive record of secular geochemical evolution throughout Earth history. Consistent with secular mantle cooling, compatible and incompatible elements in basalts record gradually decreasing mantle melt fraction through time. Superimposed on this gradual evolution is a pervasive geochemical discontinuity occurring about 2.5 Gyr ago, involving substantial decreases in mantle melt fraction in basalts, and in indicators of deep crustal melting and fractionation, such as Na/K, Eu/Eu* (europium anomaly) and La/Yb ratios in felsic rocks. Along with an increase in preserved crustal thickness across the Archaean/Proterozoic boundary, these data are consistent with a model in which high-degree Archaean mantle melting produced a thick, mafic lower crust and consequent deep crustal delamination and melting—leading to abundant tonalite–trondhjemite–granodiorite magmatism and a thin preserved Archaean crust. The coincidence of the observed changes in geochemistry and crustal thickness with stepwise atmospheric oxidation at the end of the Archaean eon provides a significant temporal link between deep Earth geochemical processes and the rise of atmospheric oxygen on the Earth.

The Earth has cooled over the past 4.5 billion years (Gyr) as a result of surface heat loss and declining radiogenic heat production. Igneous geochemistry has been used to understand how changing heat flux influenced Archaean geodynamics 1,2 , but records of systematic geochemical evolution are complicated by heterogeneity of the rock record and uncertainties regarding selection and preservation bias [3][4][5] . Here we apply statistical sampling techniques to a geochemical database of about 70,000 samples from the continental igneous rock record to produce a comprehensive record of secular geochemical evolution throughout Earth history. Consistent with secular mantle cooling, compatible and incompatible elements in basalts record gradually decreasing mantle melt fraction through time. Superimposed on this gradual evolution is a pervasive geochemical discontinuity occurring about 2.5 Gyr ago, involving substantial decreases in mantle melt fraction in basalts, and in indicators of deep crustal melting and fractionation, such as Na/ K, Eu/Eu* (europium anomaly 4 ) and La/Yb ratios in felsic rocks. Along with an increase in preserved crustal thickness across the Archaean/Proterozoic boundary 6,7 , these data are consistent with a model in which high-degree Archaean mantle melting produced a thick, mafic lower crust and consequent deep crustal delamination and melting-leading to abundant tonalite-trondhjemitegranodiorite magmatism and a thin preserved Archaean crust. The coincidence of the observed changes in geochemistry and crustal thickness with stepwise atmospheric oxidation 8 at the end of the Archaean eon provides a significant temporal link between deep Earth geochemical processes and the rise of atmospheric oxygen on the Earth.
Although it is generally assumed that mantle potential temperatures were higher in the Archaean, the ratio of terrestrial heat production to heat loss through time is very poorly constrained 8 . The resulting lack of confidence in models of secular thermal evolution has contributed to uncertainty regarding the style of Archaean geodynamics, given the temperature dependence of rock deformation mechanisms, mantle and crustal melting, and magma geochemistry 1 . The existence of highdegree-melt rocks such as komatiites, the predominance of the highsodium tonalite-trondhjemite-granodiorite (TTG) series, and the purported structural evolution of granite-greenstone terranes have all been used to suggest that plate tectonics may have been different in style or even nonexistent during the Archaean 9 . However, thermobarometric estimates of crustal geotherms similar to those for modern continents, geochronological evidence for tectonic-scale underthrusting, and the reported occurrence of subduction signatures such as high-field-strength element depletion have led other authors to argue that plate tectonics has existed in its present state since at least 3.0 or as early as 3.8 Gyr ago 9,10 .
Much like those generated in modern tectonic environments, rocks from preserved Archaean crust are characterized by extreme geochemical heterogeneity. To obtain systematic trends from this record, studies of igneous geochemistry over Earth history have typically relied on small, carefully selected data sets from well characterized terranes 3,11,12 . We take a different approach, reasoning that a sufficiently large, representative data set can be used to define global geochemical averages very precisely; and that changes in this average over time should reflect the influence of long-term geodynamic drivers such as mantle temperature.
To this end, we have constructed a database of ,70,000 major and trace element whole-rock analyses from pre-existing sources 3,4,11,13 that report sample location as well as crystallization age and uncertainty. These data are superimposed on existing geophysical models (Crust 2.0 14 and TC1 5 ) that estimate present-day crustal and lithospheric parameters. Results for global mean geochemistry with time are generated by Monte Carlo analysis with weighted bootstrap resampling (Methods; Supplementary Fig. 1). Resampling weights are inversely related to spatiotemporal sample density, to achieve a maximally uniform posterior sample density distribution and minimize selection bias ( Supplementary Fig. 2). The result is a best estimate of the average igneous geochemistry of exposed continental crust through time, reported as means with associated 2-standard-error (95% confidence interval) uncertainties of the mean for 100-Myr intervals between 0 and 3.8 Gyr. In order to distinguish between primary mantle melting and subsequent evolution and fractionation processes, we extract results from two silica ranges: 62-74% SiO 2 for intermediate to felsic samples, and 43-51% SiO 2 for less fractionated mafic samples.
The results for the least fractionated (43-51% SiO 2 ) mafic samples preserved in the continental igneous rock record reveal a systematic decrease in compatible element concentration and a corresponding increase in incompatible element concentration through time (Fig. 1). We can show that these trends reflect decreases in mantle partial melting through time, as a consequence of secular mantle cooling; alternative explanations such as a change in the degree of fractionation through time are incompatible with the observed data (Methods; Supplementary Figs 8,9). Although such low-silica samples do not generally represent primitive melts, compositional trends in these samples nonetheless reflect source processes, given a representative data set, allowing us to estimate minimum percent mantle partial melting through time (Methods). In order to derive such minimum estimates, we applied least-squares optimization to match thermodynamic batch melting simulations of primitive mantle from MELTS 15,16 with Monte Carlo simulations of basaltic major-element concentrations through time (Methods). The results (Fig. 1g) show a systematic decrease in mantle melt fraction from ,0.35 in the Archaean to ,0.1 today. These results agree well with studies of mantle xenoliths which suggest that Archaean subcontinental mantle lithosphere was produced by extraction of ,40% melt 17 .
In addition to the gradual trend expected from secular cooling, our results show a sudden, systematic decrease of ,10 percentage points (from ,27% to ,17%) in the average degree of mantle melting at ,2.5 Gyr (Fig. 1g). Due to the Earth's thermal inertia, it is unlikely that this event is the result of a rapid change in mantle potential temperature. Although a rapid change in the style of tectonics and crust formation at ,2.5 Gyr could provide one possible explanation for the rapid decrease in average mantle melting, an abrupt response in degree of melting to a gradual cooling path may provide a more parsimonious explanation; both experiments and thermodynamic models indicate substantially nonlinear mantle melting as a function of temperature, with changes in slope corresponding to clinopyroxene exhaustion 18 . Trace-element geochemistry provides some additional insight: simple trace-element modelling (Methods) suggests that the marked increases in La/Yb (Fig. 1f), La/Sm and Sm/Yb cannot be explained by the trend in degree of mantle melting in Fig. 1g alone, indicating a change in either trace-element composition of the source or effective bulk partition coefficients.
Geochemical trends in high-silica (62-74% SiO 2 ) rocks parallel many of those observed in basalts. Sodium, however, represents a notable exception. The high Na/K (Fig. 2a) and Na/(Ca 1 Na 1 K) of the TTG series have pervaded the discussion of Archaean crust formation 3,4,11,[19][20][21] . These geochemical signatures, along with high ratios of light to heavy rare-earth elements (for example, La/Yb; Fig. 2b), have been used to argue for the formation of TTG magmas by high-pressure (.40 km depth) partial melting of metabasalt at pressures above the Na-plagioclase stability field in the presence of restitic amphibole and/or garnet 3,4,11,[19][20][21] . Experimental data for partial melting of mafic rocks at mid-to lower-crustal pressures and temperatures have demonstrated a positive correlation between pressure and sodium concentration in the melt 20 . Trace-element results from our data set show correlations predicted by such a model: the secular evolution of Eu/Eu* and Sr (Fig. 2c,d) parallels that of Na/K in felsic rocks, arguing that deep crustal melting (and/or high-pressure fractional crystallization 22 ) played a larger role in forming preserved Archaean crust than post-Archaean crust.

LETTER RESEARCH
The two most frequently proposed settings for such deep crustal anatexis and TTG generation are the melting of subducted oceanic crust; and the melting of overthickened mafic lower crust, which derives its basalt from mantle melting and crustal underplating 19,21 . Of these, the latter provides a potential genetic link between TTG generation and observed high-degree mantle melting in the Archaean; extensive basalt production could lead to a thicker and/or more mafic lower continental crust, which would more often result in TTG generation. However, any process requiring a thick Archaean crust creates an apparent dilemma: preserved Archaean lithosphere is characterized by unusually thin crust 6 (Fig. 3a). A possible solution is crustal delamination, in which crustal thickening leads to a lower crust of dense garnet pyroxenite that can break off and founder into the mantle-triggering an influx of hot asthenosphere that leads to lowercrustal anatexis, crustal uplift and magmatism 23 . Overthickening of mafic lower crust thus leads to delamination, crustal thinning and further deep crustal magmatism. Indeed, previous workers have suggested that lower-crustal delamination was more prevalent in the Archaean, a process that would have been facilitated by high Archaean mantle potential temperatures 21,24 . Furthermore, because delamination involves a pressure-and temperature-dependent critical crustal thickness 24 , delamination may relate gradual secular cooling to rapid change in the geochemical record. Archaean delamination thus connects high-degree mantle melting, TTG abundance, dramatic geochemical shifts near the Archaean/Proterozoic boundary and thin preserved Archaean crust.
Although the timing of the observed geochemical transition has an uncertainty of about 6100 Myr, our results nonetheless reveal a striking temporal coincidence between an important discontinuity in the character of preserved continental lithosphere and the atmospheric great oxidation event (GOE) near the Archaean/Proterozoic boundary (Figs 1-3). It is generally agreed that the atmosphere was mostly anoxic before ,2.4 Gyr, as inferred from measurements of mass-independent sulphur isotope fractionation in pyritic shales 25 , occurrence of banded iron formations 26 , and the presence of reduced detrital grains (sulphides) in Archaean sediments 27 (Fig. 3e). But the cause of subsequent atmospheric oxidation remains a matter of debate; some indications suggest that the GOE did not occur for at least several hundred million years after the evolution of photosynthesis, leading some to argue that an additional geological stimulus was required to trigger widespread atmospheric oxidation 28,29 .
Proposed mechanisms that could link the deep Earth to atmospheric oxidation at the end of the Archaean include a change in magma degassing style (subaerial versus subaqueous) 28, or a change in mantle oxygen fugacity 29 , both of which would influence the atmosphere via the flux of reduced volcanic volatiles. The former process is possible, but it is not immediately clear how this mechanism would link to our observed geochemical trends. As for the latter process, studies of mantle redox proxies such as V/Sc, including our own analysis ( Supplementary Fig. 11), indicate that upper-mantle oxygen fugacity has remained roughly constant over the past 3.8 Gyr (ref. 12). However, despite providing useful constraints on mantle oxygen fugacity, redox tracers such as V/Sc best record changes in oxygen fugacity of the mantle/primitive melt system, not subsequent magma fractionation and evolution; consequently, constraints on the redox state of magmas during evolution and fractionation in the crust are largely absent 30 .
Magmatic redox fractionation is controlled by the ratio of mineral/ melt distribution coefficients, D Fe 31 /D Fe 21 , which are very poorly characterized. However, as shown by experimental iron speciation studies, eclogitic minerals such as garnet and sodic clinopyroxene can accommodate substantial Fe 31 , and have moderately elevated Fe 31 /Fe 21 ratios compared with other common restitic minerals such as olivine 31,32 . Consequently, melt extraction from an eclogitic or garnet pyroxenitic restite, as is expected during TTG genesis, could lead to lower melt Fe 31 /Fe 21 relative to melt extracted from a low-pressure restite rich in olivine (which contains substantial Fe but very low Fe 31 /Fe 21 ) 31,32 and plagioclase (which contains high Fe 31 /Fe 21 but very little total Fe) 32 . TTG magmatism may therefore leave more oxidizing power in eclogitic residue in the form of Fe 31 , while more reducing power is transported to the surface in the form of Fe 21 and equilibrated volatile species such as H 2 S, CH 4 and CO.
In this context, Archaean deep crustal magmatism provides a potentially significant influence on surficial oxygen fugacity: simple stoichiometric calculations show that a 15% change in restitic Fe 31 /SFe for a modest magmatic flux of 5 km 3 yr 21 would result in an increase in magma reducing capacity by approximately 4.0 3 10 12 mol e 2 yr 21 , similar in magnitude to the entire present-day volcanic volatile reducing flux (Methods). Although it is not simple to determine if the termination of this additional reducing flux at the end of the Archaean would be sufficient to trigger atmospheric oxidation, the magnitude of the flux may be constrained by further study of restitic iron partitioning.

RESEARCH LETTER
Though speculative, restitic redox partitioning provides a testable mechanism to explain the correlation of atmospheric oxidation with the ,2.5-Gyr geochemical transition revealed by our analysis. It is not yet clear whether these changes represent a fundamental shift in the style of crust formation and tectonics, or simply an abrupt response of melt production to gradual secular cooling. However, delamination and subsequent lower-crustal anatexis provide potential explanations for the observed geochemical trends in felsic rocks that are largely independent of the style of Archaean tectonics. Our analysis illustrates the opportunity presented by recent efforts to make large amounts of published geochemical data readily available and amenable to statistical analysis. These initial results substantiate the existence of a period of profound change in the Earth system near the Archaean/Proterozoic boundary, with implications for crustal differentiation, atmospheric oxidation and the habitability of the Earth for complex aerobic biota.

METHODS SUMMARY
We assembled a database of major-and trace-element analyses of ,70,000 igneous whole-rock samples from a number of pre-existing sources 3,4,11,13,20 . Each sample was paired with a geospatial location, supplemented by new estimates when necessary. Geochemical data for individual samples were then coupled with crustal and lithospheric thickness, and crustal seismic velocities (v P and v S ) and density, by projecting sample locations onto high-spatial-resolution global crustal and lithospheric models 5,14 . Sampling and preservation bias were mitigated by weighted bootstrap resampling, with sample weights inversely proportional to spatiotemporal sample density (see full methods). Estimates of all variables of interest were generated by Monte Carlo analysis as follows: (1) A subset of the data was randomly selected so that the probability of inclusion in the resampled subset is directly proportional to sample weight. (2) A synthetic data set for each sampled data point was drawn from a Gaussian distribution with a mean equal to the original value of the data point and standard deviation equal to the 1s uncertainty of the data point. (3) The resulting data were sorted into 100-Myr bins, with a mean and variance calculated for each variable. (4) Steps 1-3 were repeated 10,000 times. (5) A total mean and standard error of the mean were calculated for each variable in each bin. Mantle melt fraction (per cent melt) was estimated for each sample in the range 43-51% SiO 2 , in all 10,000 simulations, by numerical leastsquares optimization of the simulated major-element concentrations to those for slightly hydrous (0.1%) primitive mantle partial melts, as calculated by MELTS 15,16 for isobaric melting at a range of pressures. The resulting trends in per cent melt estimates agree regardless of melting pressure, suggesting that the assumption of an arbitrary constant pressure does not bias the results.
Full Methods and any associated references are available in the online version of the paper at www.nature.com/nature.

METHODS
Database compilation. We assembled a database of major-and trace-element analyses of ,70,000 igneous whole-rock samples from a number of pre-existing sources, including ,2,500 from Condie 3,11 , ,1,500 from Moyen 4,20 , and the remaining majority from the EarthChem repository (encompassing samples from Navdat, Georoc and the US Geological Survey) 13 . Each sample was paired with a geospatial location, with pre-existing location data supplemented by new estimates when necessary. Geochemical data for individual samples were then coupled with crustal and lithospheric thickness, and crustal v P , v S and density, by superimposing sample locations onto high-spatial-resolution global crustal and lithospheric models CRUST 2.0 14 and TC1 5 , respectively. To minimize the chance of introducing any unintentional selection bias, outliers were not removed from the database except in cases where the recorded data were demonstrably unphysical. The result was a data set including up to 98 defined variables (for example, sample age, SiO 2 , La, crustal thickness, etc.) for each of the ,70,000 samples, for a total of more than 3.9 million defined data points. Monte Carlo analysis. Estimates of all variables of interest were generated by Monte Carlo analysis with weighted bootstrap resampling, as illustrated for one case in Supplementary Fig. 1. Sample weights were assigned to be inversely dependent on spatiotemporal sample density, according to the relationship where n is the number of samples in the database, z is spatial location, t is age of the rock, and a and b are normalization coefficients of 1.8 arc degrees (200 km) and 38 Myr, respectively. After calculation of weight W i for each sample i in the database, bootstrap resampling was carried out in the following steps: (1) A subset of data was randomly selected so that the probability of inclusion in the resampled subset is directly proportional to sample weight. (2) A synthetic data set for each sampled data point was drawn from a Gaussian distribution with a mean equal to the original value of the data point and standard deviation equal to the estimated 1s uncertainty of the data point. (3) The resulting data were sorted into 100-Myr bins, with a mean and variance calculated for each variable. (4) Steps 1-3 were repeated 10,000 times. (5) A total mean and standard error of the mean were calculated for each variable in each bin. To distinguish between processes of mantle melting and subsequent evolution and fractionation, the entire procedure was conducted separately for samples from each of four silica ranges (43-51%, 51-62%, 62-74% and 74-80%).
Owing to the large quantity of random numbers required for Monte Carlo simulation, randomness in steps 1 and 2 was simulated with the Mersenne Twister pseudorandom number generator 34, as implemented by MATLAB. To assess the efficacy of the resampling procedure, prior and posterior density distributions were compared to the ideal case for latitude, longitude and time ( Supplementary Fig. 2). The results suggest substantial improvement in sampling quality, although tradeoffs exist between uniform sampling through time and uniform sampling in space. To evaluate the applicability of mean and standard error as appropriate statistics to characterize the results of Monte Carlo simulation, the distribution of sample means was examined and found to be approximately Gaussian, as illustrated for one variable in Supplementary Fig. 3. This is expected as a consequence of the Central Limit Theorem; regardless of the distribution of a data set, the distribution of the mean of the data set will naturally tend towards a Gaussian distribution. Finally, to mitigate the undesirable statistical properties of ratios with small denominators and minimize the effect of outliers, geochemical ratios such as La/Yb were calculated first as fractions (for example, La/(La1Yb)) during Monte Carlo analysis and averaging, and subsequently converted to ratios for presentation; the resulting 'fractional' averages are more robust than but not directly comparable to the average of a simple ratio.
In contrast to the widely discussed episodic nature of crust formation (for example, ref. 35), our analysis attempts to produce a uniform record by reducing the weight of abundantly sampled time periods or terranes. For instance, in accordance with equation (1), rock samples from a time period poorly represented in the database will be bootstrap-resampled more heavily than those from periods where much crust is preserved. However, as equation (1) attempts to balance optimal sampling in space and optimal sampling through time, periods where very few data are preserved (for example, 2.6-2.2 Gyr) will remain more poorly sampled, resulting in larger standard errors of the mean at these times. Additionally, spatiotemporal sample weighting in conjunction with Monte Carlo sampling of rocks with large reported age uncertainties result in a continuous geochemical record that smooths out stochastic and stepwise transitions; the slopes of any abrupt trends presented herein are therefore minimum estimates. The nature of the Monte Carlo sampling process has the additional consequence that there is no uniform or single-valued age uncertainty for the resulting trends. However, temporal uncertainty in the results is generally on the order of 100 Myr. Response of the Monte Carlo technique to the distribution of the uncertainty of the data. Although our Monte Carlo analysis does not assume any specific statistical distribution for the geochemical, geochronological or geophysical data, it is necessary to specify the distribution of the uncertainties associated with each individual data point in our database. As this uncertainty will be a product of the convolution of multiple independent error sources, it is reasonable to conclude that the distribution of the uncertainty will approach a Gaussian distribution, in accordance with the Central Limit Theorem. (A similar rationale explains why the distribution of the mean of the data approaches a Gaussian distribution (Supplementary Fig. 3), even though the distribution of the data is not Gaussian.) In the sources from which we compiled the database, geochronological data are implied to be Gaussian in accordance with standard practice; no uncertainty is reported for individual geochemical data, but uncertainties for the analytical methods used are typically on the order of a few per cent. Consequently, we have conducted Monte Carlo analysis using the reported geochronological uncertainty for each data point, and assuming a 62s uncertainty of 4% with a Gaussian distribution for the geochemical data, as described in 'Monte Carlo analysis', above. To ensure that the this choice of distribution does not bias the distribution of the results, two numerical experiments were carried out.
First, to demonstrate that the choice of error distribution used in Monte Carlo analysis does not bias the distribution of the data or its mean in the onedimensional case, three Monte Carlo experiments were conducted using: (1) a Gaussian distribution; (2) a log-normal distribution; and (3) a x 2 distribution with 1 degree of freedom. Using two test data sets (MgO content of mafic igneous samples between 400 and 410 Myr, and a synthetic log-normal distribution), a 10,000-run Monte Carlo simulation was then conducted using each distribution for each data set to simulate the effect of 4% (2s) analytical error-the same analytical error as is assumed for geochemical data in the Monte Carlo analysis of our 70,000-sample geochemical data set. As seen in Supplementary Fig. 4, for both the MgO data (a-d) and the simulated distribution (e-h), although the choice of error distribution may influence the small-scale smoothness of the resulting data distribution, there is no bias to the overall distribution (as particularly evident in e-h), and no bias to the mean of the distribution (red vertical lines).
To further demonstrate that the choice of statistical distribution to represent analytical error does not bias our Monte Carlo results, we conducted our full Monte Carlo analysis (as above) for three different combinations of error distributions. Supplementary Fig. 5 shows the results as a typical example for MgO in mafic samples (43-51% SiO 2 ) using: (a) a Gaussian distribution for both geochemical error and geochronological error (as in Figs 1-3); (b) a Gaussian distribution for geochronological error and a log-normal distribution for geochemical error; and (c) a log-normal distribution for both geochemical error and geochronological error. As seen in Supplementary Fig. 5, the resulting trend in MgO through time is nearly identical in each case. We conclude that the choice of error distribution in our Monte Carlo analysis has little to no effect on the resulting trends. Response of the Monte Carlo results to regions of missing data. Due to the well known paucity of early Proterozoic igneous rocks 36 , it may be argued that the position and/or shape of the ,2.5-Gyr geochemical transition that we observe may be influenced by a lack of data in this period. In our data set, this lack of data is particularly pronounced between 2.2 and 2.6 Gyr, a period during which the temporal density of geochemical data is only ,5% that at the pronounced data peaks at 2.2 and 2.7 Gyr. Consequently, we have conducted additional tests to examine the sensitivity of the Monte Carlo technique to periods of sparse data. First, a 10,000-point synthetic data set was generated using a log-normal distribution, with an abrupt shift in the mean and variance of the data at 2.5 Gyr. Then, varying amounts of data were deleted in the period from 2.2 to 2.6 Gyr, ranging from 100% of data present to 0% of data present (all deleted). The data set was then subjected to Monte Carlo analysis as above, using an average age uncertainty of 22 Myr and 4% analytical error. As seen in Supplementary Fig. 6, an abrupt transition, as seen in Fig. 1g and Fig. 2b-d, may be preserved to some extent down to the point of 5% data remaining; with 1% and 0% of the data preserved, only a smooth transition is observed. Similarly, deleting data from a continuous starting distribution will not result in any discontinuity in the resulting Monte Carlo output, only a large increase in uncertainty in the data-poor region ( Supplementary  Fig. 7). Thus, it appears unlikely that the observed abrupt transitions in secular geochemical evolution could be a product of the scarcity of early-Proterozoic data. More gradual transitions such as in Fig. 1f may be influenced by a lack of data; although such trends still provide convincing evidence for a major geochemical transition, the timing and rapidity of that transition is less certain. Although the dearth of data between ,2.2 and ,2.6 Gyr may contribute in cases such as Fig. 1f to the uncertainty in the timing of the observed geochemical transition, it remains the case that between the uncertainty in the timing of this transition and the RESEARCH LETTER uncertainty in the timing of the Great Oxidation Event, these two events easily overlap within uncertainty. Interpretation of major element trends. Assuming our data set provides a representative sample of the crust, four candidate mechanisms could explain the observed trends (Fig. 1) in compatible and incompatible element composition of mafic rocks: (1) similar degree melting from a mantle that is more enriched now than in the past; (2) higher degrees of crustal contamination of otherwise identical basalts in the present; (3) higher degrees of fractionation in the present; or (4) higher mantle melt fraction in the past, due to either higher temperature or higher water content. A modern-day mantle that is more enriched now than in the Archaean is inconsistent with well established evidence from Nd, Sr and Hf isotope studies of mid-ocean ridge basalts 37,38 , eliminating mechanism (1). Increasing crustal contamination (2) should be reflected in ratios such as Nb/U and Pb/ Ce 39, which are instead relatively constant through time (Supplementary Fig. 8). If increasing magma fractionation (3) were used to account for the entire observed trend, this would imply that modern basalts are on average created through fractionation of very-high-MgO magmas such as komatiites, which is clearly not the case. Indeed, if anything, a thick Archaean crust as suggested by crustal depth-of-melting/fractionation indicators (Fig. 2), as well as rapidly cooling high-Mg mafic magmas, suggest increased fractional crystallization in the Archaean, in opposition to the observed trend. To address the issue more quantitatively, we mimic the approach of ref. 40, by applying our bootstrap resampling approach to reconstruct global mean MgO-Na 2 O and MgO-CaO/Al 2 O 3 trends for time periods between the Archaean and the present day ( Supplementary Fig. 9). Both trends are consistent with globally averaged fractionation paths, illustrating that substantial fractionation occurs at any given time, but parallel trends are shifted as a function of age. Following Klein and Langmuir 40 , these parallel trends illustrate changes in mantle melt fraction. Given that decreasing average mantle melt fraction through time (4) is then the most compatible interpretation of the data, we are faced with two options: secular cooling or mantle dehydration. Evidence for a more hydrous Archaean mantle is equivocal 41 , and it has even been argued that the potentially undercompensated flux of water into the mantle at subduction zones suggests that the hydration state of the mantle is increasing through time, not decreasing 42,43 . Secular cooling, on the other hand, presents no such caveats and is widely expected on the basis of physical constraints.
The degree of mantle melting was estimated for each sample with 43-51% SiO 2 , in all 10,000 simulations, by numerical least squares optimization of the simulated major-element data to that of slightly hydrous (0.1%) primitive mantle partial melts as calculated by MELTS 15,16 for isobaric melting at a range of pressure conditions. Resulting per cent melt estimates were found to agree within error despite order of magnitude differences in melting pressure ( Supplementary Fig. 10), suggesting that the results are insensitive to changes in model pressure. Although MELTS does not extend to high enough pressures to directly simulate high-degree isentropic decompression melting, as long as the system is not far from equilibrium, per cent melt and melt composition will be roughly state functions of pressure and temperature-that is, path-independent. Consequently, at a given pressure (say, corresponding to the base of the crust) per cent melt and melt composition will both be state functions of temperature. In this case, then, there will be a simple relationship between per cent melt and melt composition, allowing for accurate estimation of per cent melt from melt composition regardless of the exact melting path. This approach will break down in the case of fractional melt extraction, in which the melt cannot be in equilibrium with the solid; the absolute value of the results may therefore be inaccurate in the case of predominantly fractional melting, though relative trends would be expected to remain robust.
Because our analysis focuses on the trends in melt percentage over Earth history rather than the absolute values of these estimates, a number of sources of systematic error can be disregarded. We have intentionally not corrected basalt major-element compositions for fractional crystallization to calculate primary melt compositions, as such an analysis has already been conducted on smaller data sets, and is only possible for a small proportion of the available data 44 . By neglecting the effect of fractional crystallization, our estimates of per cent melting through time are systematic underestimates. Nonetheless, in the absence of large changes in the degree of fractional crystallization through time, as is argued above, the results should accurately represent changes in per cent melt over time. Herzberg et al. 44 analysed a database of ,1,500 basalts through early Earth history, and identified those that could be corrected to primary mantle melts (n 5 33) to extract mantle potential temperature. The topology of a temperature curve corresponding to our degreemelting curve (Fig. 1f) agrees well with Herzberg et al.'s mantle potential temperature, although, as expected, they obtain higher mantle temperatures 44 . Moreover, our per cent melting results generally agree with independent estimates of modern and Archaean percent melting, suggesting that the degree of underestimation is not large 17,45 .
To determine if observed trends in trace-element ratios such as La/Yb in lowsilica samples could be generated solely by a single set of bulk partition coefficients given the estimated per cent melting curve, we used linearized least-squares regression to optimize bulk partition coefficients to fit per cent melting and traceelement curves with analytical non-modal batch melting equations. No single set of partition coefficients can reproduce the observed trace-element ratios in basalts by varying melt fraction alone, showing that other variables such as temperature and mineralogy may have influenced partition coefficients, or that significant changes in source trace-element composition (for example, for the example of an arc setting, change in the trace-element composition of the mantle wedge due to subduction erosion or sediment subduction) had an effect. Redox partitioning calculations. To assess the potential for magmatic redox partitioning to affect surficial redox state, a simple stoichiometric redox balance of mantle restite was calculated for a modest magmatic flux of 5 km 3 yr 21 (equivalent to estimates of modern arc magmatism 48 ), assuming a magma density of 2,500 kg m 23 , 6 wt % restitic Fe, and a ratio of 2 kg of restite to 1 kg of magma, as follows: Consequently, a 15% increase in Fe 31 /SFe of this restitic flux would liberate 4.02 3 10 12 mol e 2 yr 21 , or, stoichiometrically, consume 1.0 3 10 12 mol O 2 yr 21similar in magnitude to estimates of the entire present-day volcanic volatile flux of ,5.6 3 10 12 mol e 2 yr 21 (ref. 46). Notably, the Fe 31 /Fe 21 ratio of many magmas is low even today, apparently not leaving much room for variation through time. However, theory and experiment suggest that oxygen fugacity is proportional to approximately the fourth power of ferric/ferrous ratio 47 . Thus, a change in ferric/ferrous ratio from 0.1 to 0.01 would result in a decrease of roughly ten orders of magnitude in magmatic oxygen fugacity, with potentially important consequences for the production of volcanic volatiles.