Plate tectonics and continental basaltic geochemistry throughout Earth history

Abstract Basaltic magmas constitute the primary mass flux from Earth's mantle to its crust, carrying information about the conditions of mantle melting through which they were generated. As such, changes in the average basaltic geochemistry through time reflect changes in underlying parameters such as mantle potential temperature and the geodynamic setting of mantle melting. However, sampling bias, preservation bias, and geological heterogeneity complicate the calculation of representative average compositions. Here we use weighted bootstrap resampling to minimize sampling bias over the heterogeneous rock record and obtain maximally representative average basaltic compositions through time. Over the approximately 4 Ga of the continental rock record, the average composition of preserved continental basalts has evolved along a generally continuous trajectory, with decreasing compatible element concentrations and increasing incompatible element concentrations, punctuated by a comparatively rapid transition in some variables such as La/Yb ratios and Zr, Nb, and Ti abundances approximately 2.5 Ga ago. Geochemical modeling of mantle melting systematics and trace element partitioning suggests that these observations can be explained by discontinuous changes in the mineralogy of mantle partial melting driven by a gradual decrease in mantle potential temperature, without appealing to any change in tectonic process. This interpretation is supported by the geochemical record of slab fluid input to continental basalts, which indicates no long-term change in the global proportion of arc versus non-arc basaltic magmatism at any time in the preserved rock record.


Introduction
Systematic changes in crustal composition and magma chemistry between the Archean and the present have long been advanced on the basis of geochemical, petrological, and field observations. Differences between Archean and post-Archean crust -including the preponderance of granite-greenstone terranes, high-Na Trondhjemite-Tonalite-Granodiorite (TTG) plutons, and ultramafic komatiite lavas in the Archean -have led to much discussion about the style and rate of Archean tectonics (Green, 1975;Hargraves, 1976;Wilson et al., 1978;Sleep, 1992;de Wit, 1998;Hamilton, 1998;Bédard, 2006;Pease et al., 2008;Condie and Kröner, 2008;Stern, 2008). One interpretational paradigm invokes the abundance of TTGs, komatiites, and apparent lack of preserved blueschists and ophiolites in the Archean along with vertical tectonic interpretations of dome-and-keel granite-greenstone terranes as evi- To obtain meaningful temporal trends from the heterogeneous rock record, one productive approach has been to prioritize quality over quantity, for instance obtaining a well-constrained primitive magma composition by tracing a liquid line of descent from samples that have only fractionated olivine (Herzberg et al., 2010), or by focusing on the composition of carefully screened and geologically well-understood suites of samples Condie and O'Neill (2010); Condie et al. (2016). For instance, using 33 primary melt solutions culled from a dataset of 1500 "non-arc" basalt analyses, Herzberg (Herzberg et al., 2010) found primary magma temperatures consistent with a smooth trend in mantle potential temperature from ∼1530 • C in the Eoarchean to ∼1350 • C today, along a similar trend to the models of Korenaga (2008).
If, however, one obtains a sufficiently extensive, representative sample dataset, an alternative inclusive approach becomes possible. Given a perfectly representative sample set of the continental crust, the compositional average of this set has unique physical significance corresponding to the bulk composition of the crustthat is, the composition that would result from physically homogenizing the entire crust. Consequently, changes in the average composition of new igneous crust contain information about the processes of crust formation: at long time scales, changes in the mean reflect secular variation in fundamental terrestrial parameters such as mantle potential temperature, crustal thickness, and tectonic mechanism of crust formation. A continuous temporal record of average crustal geochemistry would thus allow for the critical examination of temporal variations in such Earth system parameters and, consequently, the extrapolation of modern geologic processes back through Earth history in a conditional application of uniformitarianism. However, the quantification of systematic changes in mean crustal composition has been impeded by the difficulty of obtaining a representative sample population given the vast heterogeneity produced by magmatic processes, along with persistent questions about sampling and preservation bias Rudnick and Gao (2014).
Recently, increasing availability of large geochemical datasets, such as EarthChem, Georoc, and Navdat, has provided a new opportunity to quantify crustal composition by employing a greatly enlarged sample set. Our previous work (Keller and Schoene, 2012) has taken advantage of these freely-available datasets to construct and analyze a ∼70,000-sample whole-rock geochemical record of secular crustal evolution over the past 4 Ga, using weighted bootstrap resampling to obtain maximally representative temporal trends. First-order results revealed smooth decreases in compatible element abundance and increases in incompatible element abun-dance in basaltic samples ( Fig. 1), along with some indication of more discontinuous processes at the end of the Archean. While many indications of geochemical variation across the Archean-Proterozoic boundary were specific to the felsic record, changes in the estimated degree of average mantle melting and in basaltic La/Yb ratio circa 2.5 Ga suggested an analogous transition in the mafic record as well (Keller and Schoene, 2012).
Given the petrological complexities involved, previous work has not presented any definitive conclusions on the origin and significance of this potential transition in the mafic record. However, many additional geochemical constraints are available due to the breadth of the underlying dataset. In order to better understand the implications of the new record of crustal evolution provided by statistical geochemistry, we focus here on the ∼30,000-sample mafic subset of our database (defined as 43-51 wt.% SiO 2 , generally including komatiites but excluding many ultramafic cumulates), here termed basaltic sensu lato. As the product of 4 Ga of partial mantle melting, this mafic record preserves the most direct account of mantle conditions through time, and provides the necessary foundation for any interpretation of subsequent processes of crustal differentiation. Our analysis of this record combines weighted bootstrap resampling with major-and trace-element geochemical modeling. This analysis allows us to study global variations in the process and setting of basaltic magma production over geologic time, in turn providing insight into the style and existence of plate tectonics throughout Earth history.

Weighted bootstrap resampling
In order to minimize sampling bias and accurately propagate temporal and compositional uncertainties to obtain an optimal estimate of the compositional distribution of Earth's exposed crust, we followed the weighted bootstrap resampling method of Keller and Schoene (2012), summarized in Appendix A and Supplementary Figs. 1-6. As in our previous work, this technique was applied to a ∼70,000-sample record of crustal whole-rock major and trace element geochemistry, including ∼67,000 samples from the EarthChem repository, ∼2,500 samples from Condie and O'Neill (2010), and ∼1500 from Moyen (2010) including age estimates and geospatial sample locations for each sample. Although outliers were not rejected categorically, the dataset was manually screened to identify suspicious data patterns and exclude any physically implausible data (e.g. concentrations <0 or >100%). We consider only samples collected from the present-day continental crust (i.e., excluding MORB and OIB except where naturally ob-ducted), focusing unless otherwise noted on samples with between 43-51 wt.% SiO 2 . This silica range is chosen to include a majority of mantle-derived mafic magmas (including komatiites), but exclude many ultramafic cumulates; however the exact basaltic silica range chosen has relatively little influence on the results (Supplementary Fig. 6). To minimize filtering-induced interpretational circularity, no additional geochemical screening was conducted. The full dataset and all computational source code is available at https :/ /github .com /brenhinkeller /StatisticalGeochemistry under an MIT open-source license.

Gradual secular variation in continental basaltic composition over Earth history
The results of this statistical analysis (Figs. 1-6) reveal a diverse array of trends in average basaltic composition well beyond the two-standard-error level, ranging from smooth and continuous over 4 Gyr to abrupt and discontinuous at the scale of our 100 Myr binning resolution. Considering first the continuous trends, we observe what may be summarized as a consistent increase in incompatible element abundance accompanied by a decrease in compatible element abundance for average basalts and gabbros of a constant silica range (43-51% SiO 2 ) over the last 4 Ga. In addition to the observations of Keller and Schoene ( Fig. 1), such trends persist for a wide range of compatible and incompatible elements ( Supplementary Fig. 5), hinting at an igneous origin for these firstorder variations.

Sampling, preservation, and alteration
In principle, systematic changes in average continental basalt composition may represent igneous variations in average crustal basaltic composition, secondary variations in the composition of preserved basalts due to increasing alteration or metamorphism of older samples, or illusory variations driven only by sampling or preservation bias. Evidence that variations in basaltic compatible and incompatible element abundances since 4 Ga represent magmatic signatures, and are not overly influenced by alteration or weathering, is provided by the matching variations in element pairs with similar olivine/melt partition coefficients but with very distinct chemical properties. Fig. 2 illustrates a number of such element pair ratios: U and La, both highly incompatible but with dramatically different redox chemistry; K and Sm, both incompatible but very distinct in abundance and aqueous fluid solubility; and Sc and Co, both compatible but with distinct size, charge, and redox chemistry. Differences in ionic radius, charge, and redox chemistry will greatly influence elemental mobility and solubility during alteration and metasomatism, yet all ratios remain nearly invariant over the last 4 Ga.
Further, while sampling and preservation bias are difficult to eliminate even with our weighted bootstrap resampling methods (Appendix A), there is no known a priori reason that either sampling or preservation bias should favor more mafic compositions in the past and more felsic compositions towards the present. While difficult to quantify given its dependence on the collective sampling behavior of the geological community, the expected distribution of sampling bias in silica space is likely to disproportionately favor less common rock types -acting counter to any underlying changes in lithology. Meanwhile, preservation bias should favor (if anything) more felsic samples over time, given the increased chemical weathering rate and greater density of mafic rocks, which make mafic crust more vulnerable both to erosional loss at the surface and convective loss at the base of the crust. Here and throughout, ratios are calculated using the method of Keller and Schoene (2012), starting with fractions a/(a + b) (to avoid numerical instability) for each whole-rock analysis with both elements defined, then resampled and averaged before converting back ratio space a/b.

Variation in the extent of mantle melting
Given an igneous interpretation, the most straightforward hypothesis to produce the systematic temporal trends in basaltic compatible and incompatible element abundance observed by Keller and Schoene (2012) is that secular mantle cooling has decreased the average extent of mantle partial melting through time, resulting in increasing incompatible element concentrations and decreasing compatible element concentrations as seen in Fig. 1. However, a range of confounding factors could bias these trends, including changes in mantle depletion, basaltic magma differentiation, and magma water content.

Mantle depletion
Since Earth's mantle forms the depleted complement to the silicic continental crust, the mantle has become depleted in incompatible elements as a result of crust formation (Hofmann, 1997). While terrestrial crustal growth rates remain actively debated, nearly all models agree that the crust has either grown or remained constant in volume since the end of the Hadean Harrison (2009). Constant or increasing crustal volume would result in either constant or increasingly depleted mantle composition. Any increase in mantle depletion would not, however, explain the observed trends: an increasingly depleted mantle would result in decreasing incompatible element concentrations in basalts over time, opposite the trends observed from our dataset.

Differentiation
A more subtle point concerns the possibility for changes in the extent or character of basaltic magma differentiation. For many elements, increasing magma differentiation by fractional crystallization and decreasing degree of primary mantle melting produce Fig. 3. Deconvolving the effects of mantle melting and fractional crystallization. Observed trends throughout this figure were produced by binning the full dataset of Keller and Schoene (2012) -by SiO 2 in a-b and by MgO in c-f. Independent variable uncertainties are smaller than the symbol size, with binning ranges chosen to illustrate a region of petrologic interest in each subfigure; full methodology is described in the text. (a) MgO decreases as a function of silica content during fractional crystallization, converging at high SiO 2 . Reductions in mantle melting over the last 4 Ga decrease MgO at nearly constant silica, producing offsets between different age differentiation trends at high angle to the differentiation trend. (b) FeO*/MgO ratio increases with silica during fractional crystallization, but decreases at nearly constant silica during partial melting, producing systematic offsets between fractionation trends of different ages that cannot be explained by changes in basalt differentiation. FeO* represents total Fe 2+ + Fe 3+ normalized as FeO. (c) Relative effects of partial melting and fractional crystallization on FeO* as a function of MgO. (d) Changes in mantle partial melting and basalt fractional crystallization have distinct slopes in CaO/Al 2 O 3 -MgO space, producing resolvable offsets in differentiation trends over time. (e) The systematic offset between Archean and modern MgO-TiO 2 differentiation trends is not well explained by changes in either extent of mantle melting or fractional crystallization, but may instead reflect systematic changes in depth of mantle melting or reequilibration. (f) Similar offsets between Archean and later MgO-Na 2 O differentiation trends again suggest systematic variation in melting and/or equilibration pressure over the past 4 Ga. Notably, Archean basalts have lower Na 2 O for a given MgO, despite the well-known high Na signatures of Archean granitoids. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.) similar changes in elemental abundance. As a result, differentiating the contributions of these two processes can be problematic. In order to resolve such processes, we turn to a variety of differentiation diagrams featuring elements or ratios that behave differently during mantle partial melting as opposed to basalt fractional crystallization (Fig. 3).
It is important to note that the differentiation trends shown in Fig. 3 are not liquid lines of descent, but rather the average of many discrete lines of descent from different magmas under different conditions. As such, the interpretation is distinct from a traditional liquid line of descent. Nonetheless, changes in the average extent of mantle melting or average trend of differentiation will remain clearly visible. Due to the heterogeneity of mantle-derived basalt compositions, a single "primitive" origin is not shown, but in general mantle-derived basalts will originate on the left and differentiate to the right in Fig. 3a-b, and originate on the right and differentiate to the left in 3c-f, as illustrated by the schematic "fractional crystallization" arrows. The directions of the accompanying mantle "melting extent" and mantle "melting pressure" were established by pMELTS-mode alphaMELTS simulations (Ghiorso et al., 2002;Asimow et al., 2004;Smith and Asimow, 2005) of isobaric batch mantle melting from solidus to liquidus at a range of pressures between 1 and 3 GPa using the primitive mantle composition of McDonough (McDonough and Sun, 1995) with varying bulk mantle water contents between 0.1% and 1% H 2 O (Appendix B; Supplementary Fig. 7). Variation in water content over this range resulted in substantial differences in primitive magma composition for some elements, but neither varying water content nor varying oxygen fugacity substantially changes the overall vector of mantle melting extent and melting pressure at the scale of the schematic effect vectors.
As seen in Fig. 3, differentiation trends for successive 1-Gyr periods of Earth history define parallel arrays, indicating changes in mantle partial melt composition perpendicular (in elemental crossplot space) to the prevailing differentiation trend (c.f. Klein and Langmuir, 1987). In particular, changing degree of mantle melt appears to have little effect on primary magma silica content relative to other major elements such as MgO (e.g., Supplementary Fig. 7), such that binning basaltic magmas by silica serves as a relatively robust way to extract the effects of mantle partial melting for magmas which have undergone similar degrees of differentiation (Fig. 3).

Water content
Ruling out changes in mantle depletion or basaltic differentiation, a final caveat stems from the ability of water to decrease mantle solidus temperature (e.g. Mysen and Boettcher, 1975;Grove et al., 2006). Decreasing mantle water content could in principle lead to decreasing mantle melting extent for a constant potential temperature, with solidus temperature increasing in response to ongoing mantle degassing. However, changing mantle water content also has major predictions for subsequent magma differentiation. Notably, increased magmatic water activity promotes magnetite saturation relative to the major rock-forming silicate min- Fig. 4. Estimated average extent of mantle melting (a) and mantle potential temperature (b) through time, determined by least-squares comparison of observed major element concentrations for each basaltic sample in the database (43-54% SiO 2 ) to those of calculated basalt compositions corresponding to 0 to 100% partial melting of primitive mantle. Mantle melting was simulated by pMELTS-mode alphaMELTS isobaric equilibrium batch melting calculations (Ghiorso et al., 2002;Asimow et al., 2004;Smith and Asimow, 2005) using the primitive mantle composition of McDonough (McDonough and Sun, 1995) with 0.15 wt.% initial H 2 O at 2.2, 2.7 and 3.2 GPa. Compared to a record derived from oceanic basalts, this continental basaltic record is less directly representative of the freely convecting mantle, influenced by, e.g., conditions in the arc mantle wedge. erals, changing magma Fe/Mg ratios during differentiation (Sisson and Grove, 1993;Botcharnikov et al., 2008). As demonstrated in Fig. 3b-c, neither the FeO*/MgO* slope during differentiation nor the relative occurrence of magnetite saturation in FeO*-MgO space is significantly altered over the past 4 Ga. Consequently, declining extent of mantle melting in response to secular mantle cooling remains the simplest explanation for the observed first-order variations in average basaltic geochemistry through time.

Melting extent
If it is possible to calculate the approximate extent of mantle melting for each basaltic sample in the database, we may use this information to produce a temporal record of global average mantle melting extent. Along these lines, one first-order methodology is to compare the major element composition of each whole-rock sample to pMELTS simulations of primitive mantle melting and find the (least-squares) best-fit extent of mantle melting for each basaltic sample. A representative average record through time is then produced by the usual weighted bootstrap resampling technique (Keller and Schoene, 2012). This approach naturally produces a corresponding estimate of mantle melting temperature for each sample, allowing for the calculation of apparent mantle potential temperature through time. This first-order approach makes no attempt to either exclude basaltic samples that have undergone small degrees of fractional crystallization or crystal accumulation, or to numerically correct back to a primitive magma composition. Since most basaltic magmas have undergone some degree of fractional crystallization, this approach likely results in an underestimation of mantle melting extent. However, in the absence of changes in degree of fractional crystallization through time, this underestimation results in only a systematic offset with little influence on relative temporal trends.

Mantle potential temperature
The results (Fig. 4) are recalculated using this approach for isobaric batch melting at pressures between 2.2 and 3.2 GPa using the primitive mantle composition of McDonough (McDonough and Sun, 1995) with 0.15 wt.% initial H 2 O. Percent melt estimates display a clear decrease from nearly 40% mantle melting in the Archean to ∼10% today, in agreement with independent estimates of modern and Archean mantle melting (Kelley et al., 2010;Portnyagin et al., 2007;Walter, 1998;Bernstein et al., 2007) despite the lack of a fractionation correction. Over the same time period, the corresponding mantle potential temperature record cooled from approximately 1500 • C in the Archean to 1330 • C at present, following a trend intermediate between the models of Korenaga (Korenaga, 2008) and Davies (Davies, 2009), which are based on scaling laws for the mantle heat budget through time. Although average basaltic water content appears not to have changed substantially over the past 4 Ga (no change in FeO*/MgO differentiation slope in Fig. 3b-c, discussed further in section 4.1), the same cannot be said for mantle melting pressure; systematic offsets in MgO-TiO 2 and MgO-Na 2 O differentiation trends ( Fig. 3e-f) likely require systematic variation in melting, equilibration or differentiation pressure; a comparison to MELTS simulations would suggest increasing average mantle melting pressure through time, though this inference is not definitive (Appendix B). If present, a systematic increase in mantle melting pressure would allow a smaller range of mantle cooling to explain the observed decrease in mantle melting extent. Nonetheless, the total magnitude of average mantle cooling (∼170 • C) is generally intermediate between those suggested by Grove and Parman (2004) and Herzberg et al. (2010), and broadly consistent with previous estimates despite contrasting methodology (Herzberg et al., 2010;Condie et al., 2016).
By examining basalts preserved in the continental crust, we are largely excluding oceanic basalts such as MORB and OIB. Instead, considering the results of Section 4.1, the temperature trends of Fig. 4 are arguably most comparable to the hydrated mantle basalt estimates of Condie et al. (2016). However, in many other ways the results display only limited sensitivity to model assumptions. Varying the sample set by rejecting samples with komatiitic MgO contents only slightly diminishes estimated temperatures and melting extents ( Supplementary Fig. 8), since excluding komatiites has little effect on most major elements other than MgO. Excluding samples with high least-squares residuals similarly has little effect on the results (Supplementary Fig. 9). An equivalent model employing continuous polybaric fractional melting in place of isobaric batch melting ( Supplementary Fig. 10) results in higher estimated potential temperature due to rapid removal of any initial water during fractional melting, but has almost no impact on estimated melting extent. Overall, relative trends are broadly consistent for a wide range of melting models, and while absolute mantle potential temperature estimates are necessarily model-dependent, estimated mantle melting extents are much less so (Supplementary Fig. 11; Appendix B).

Rapid geochemical variation at the end of the Archean
While the estimated decreases in melting temperature generally follow a smooth curve over the past 4 Ga, a marked offset (at the two-standard-error level) is observed in average apparent mantle percent melt across the Archean/Proterozoic boundary at 2.5 Ga (Fig. 4a). The absence of any similarly rapid shift in the corresponding temperature record (consistent with the thermal inertia of the mantle) suggests that, if representative of a real phenomenon, this offset does not stem directly from a change in mantle temperature alone, but rather from a systematic threshold in process, mineralogy, or rate.
In addition to the mantle melting extent record in Fig. 4a, a number of basaltic trace element abundances and ratios also exhibit shifts at or near the Archean-Proterozoic boundary. This effect was first observed for La/Yb ratio, where trace element calculations indicate that the increase in average basaltic La/Yb at the end of the Archean cannot be explained by changes in percent mantle melting alone without either a change in partition coefficients or source composition (Keller and Schoene, 2012). As seen in Fig. 5a-c, broadly step-like increases in rare earth fractionation (for both light/middle and middle/heavy rare earth elements) are demonstrated by average basaltic La/Yb, La/Sm, and Sm/Yb ratios. Similarly rapid (<1 Gyr) transitions are observed in the abundance of a number of incompatible high-field-strength elements such as Zr (Fig. 5d) and Ti.

Geochemical signatures of subduction
One potential model for the observed 2.5 Ga transition involves the initiation of modern-style subduction at the end of the Archean, which could increase average basaltic La/Yb, La/Sm, and Sm/Yb by addition of high La/Yb slab fluids and subducted sediment. It has been suggested that such a scenario could also produce a decrease in apparent degree of mantle melting due to the introduction of cold subducted slabs into the mantle (Condie and O'Neill, 2010). If this is the case, the same basaltic samples should also record a dramatic change in subduction signatures which record the effect of slab fluids due to the relative insolubility of high field strength elements such as Nb and Ta compared to large ion lithophiles like K and Rb.
In particular, Condie and O'Neill (2010) suggested the use of the Nb/Th ratio, which unlike many other high field strength/large ion lithophile ratios is relatively insensitive to changing percent melt due to the similar incompatibility of Nb and Th in the mantle. In particular, Nb/Th ratios below 5 are generally associated with a subduction setting, where slab fluids are abundant, while Nb/Th ratios above 7 are associated with melts produced in the absence of slab fluids. However, as seen in Fig. 6, both the average Nb/Th and the percentage of high-and low-Nb/Th continental basalts are essentially constant throughout recorded Earth history. Alternative subduction signatures such as K/Hf and Rb/Nb -also with matching peridotite/melt partition coefficients but with even higher solubility contrast (Tatsumi et al., 1986) -similarly show no major discontinuity across the Archean-Proterozoic boundary and no clear trend over geologic time (Fig. 6b-c).
Further insight into the prevalence of subduction throughout geologic time is provided by the record of calc-alkaline vs. tholeiitic magma series. More hydrous and more oxidized magmas saturate Fe-Mg oxides such as magnetite and ilmenite sooner relative to the major rock-forming silicate minerals, driving the remaining magma composition away from the Fe axis and resulting in a calc-alkaline differentiation trend (Sisson and Grove, 1993;Botcharnikov et al., 2008). Consequently, calc-alkaline differentiation is observed to correlate strongly with magmatic water content (Zimmer et al., 2010). While tholeiitic magmas may occur in multiple tectonic settings, the calc-alkaline series is restricted almost exclusively to subduction zones, where it is a dominant differentiation trend. Consequently, the initiation of subduction in the rock record should coincide with an increase in the abundance of the calc-alkaline series.
Quantitatively, Zimmer et al. (2010) have proposed the "tholeiitic index" as a numerical measure for the degree of tholeiitic vs. calc-alkaline differentiation in a given magma series, defined as total Fe at 4% MgO over total Fe at 8% MgO for magmas of the same series (e.g. as in Fig. 3c). Extending this metric to consider the average tholeiitic index of the continental crust -comparing the FeO content of magmas with 4 wt.% MgO to those with 8 wt.% MgO for each 100 Ma time bin -should allow us to discern systematic changes in the abundance of calc-alkaline versus tholeiitic differentiation. Remarkably, as seen in Fig. 6d, the average tholeiitic index of juvenile basalts have remained constant over the last ∼3.8 Ga, strongly suggesting that average continental basaltic water content has remained relatively constant for most of Earth history. These results, including an average THI of ∼0.9 are consistent with the predominantly calc-alkaline average differentiation trends observed in Fig. 3b-c.

Trace element modeling
If tectonic setting and magma water activity are constant over the last ∼4 Ga, one possible remaining source of discontinuous threshold processes in the basaltic record involves melting systematics -the modal mineralogy of the residual mantle as a function of melt extraction. Although our earlier results suggested that the observed La/Yb trends could not be fit with a constant set of partition coefficients (Keller and Schoene, 2012), bulk peridotite/melt partition coefficients are not constant during mantle melting, nor is melt productivity linear as a function of mantle temperature. Instead, as is well known from both experimental and theoretical studies, the relationship of mantle percent melt to melting temperature is markedly non-linear, with discontinuities in slope corresponding to the exhaustion of peridotite minerals upon progressive melting, most notably at the point of clinopyroxene exhaustion around 25-30% melt (Herzberg and Zhang, 1996;Hirschmann et al., 1998).
In order to investigate the effects of mantle melting systematics on melt productivity and partition coefficients, we paired the pMELTS mantle melting simulations from Section 3.2 with trace element partitioning calculations (Appendix B). To this end, we first extracted the relative mass of melt and residue present for each temperature step, as well as the proportion of each rock-forming mineral in the corresponding mantle residue. Trace element concentrations in the melt were then calculated using the equilibrium batch melting equation along with mineral/melt partition coefficients derived from the GERM partition coefficient database (GERM, 2013) interpolated as a function of SiO 2 as in Keller et al. (2015) by bootstrap resampling. To allow for processes such as slab fluid input and magmatic differentiation which influence the trace element budget more dramatically than the major element budget, the initial primitive mantle trace element composition was inverted by scaling the model to fit the mean of the basaltic record. This trace element budget was subsequently held constant for all time, providing no influence on the model trend over time. Compared to standard primitive mantle trace element estimates (e.g. McDonough and Sun, 1995), this inversion resulted in higher effective mantle trace element abundances, particularly of fluid mobile elements -consistent with an enriched slab fluid input.
The results reveal that at pressures typical for mantle wedge flux melting, garnet and clinopyroxene are both exhausted in the residuum at ∼30% melt (Fig. 7a), resulting in a discontinuous break in slope for both melt productivity (Fig. 7b) and peridotite/melt partition coefficients (Fig. 7c). The nonlinearity and discontinuous slope of the percent melt-temperature relationship in Fig. 7b are responsible for reconciling the discontinuous percent melttime curve in Fig. 4a with the smooth temperature-time curve in Fig. 4b. Applying the partition coefficients from Fig. 7c to the percent melt-time trend from Fig. 5a accurately reproduces observed trace element trends, as seen in Fig. 8, with no fitting or adjustment of partition coefficients. As with the mantle melting  Fig. 7, GERM partition coefficients, and a constant mantle composition. estimates of Section 3.2, calculated trace element trends are not substantially altered by substituting continuous fractional polybaric melting in place of isobaric batch melting ( Supplementary Fig. 10).

Constraints on the prevalence of subduction through Earth history
The initiation of subduction and modern-style plate tectonics would directly precipitate an increase in average basaltic water content due to the addition of wet basaltic magmas produced by flux melting with ∼4.0% H 2 O (e.g. Plank et al., 2013) to the preexisting population of dry rift and plume basalts produced by decompression melting with only ∼0.2% H 2 O (e.g. Sobolev and Chaussidon, 1996), corresponding to a decrease in tholeiitic index from 1.8 at 0.2% H 2 O towards 0.8 at 4.0% H 2 O (Zimmer et al., 2010). However, no such trend is observed in our results (Fig. 6d). Instead, average tholeiitic index has held roughly constant at THI 0.9 since the earliest Archean, indicating a consistent predominance of calc-alkaline differentiation (Fig. 6d). This conclusion is strongly supported by a wide range of fluid-sensitive trace element signatures including Nb/Th, K/Hf, and Rb/Nb (Fig. 6a-c, e), all of which show no systematic variation in the relative extent of slab fluid input since the earliest Archean. Indeed, since arc basalts today display higher degrees of melting than coeval rift and plume basalts (Watson and McKenzie, 1991;Kelley et al., 2010) due to the increased water content of arc basalts, the sudden addition of hydrous arc basalts to the geologic record would arguably result in a rapid increase in average mantle melting extent, which is not observed at any time in the preserved rock record (Fig. 4).
Objections to Archean plate tectonics have been raised on a wide variety of grounds, notably among them the absence from the preserved Archean record of certain rock types considered characteristic of plate tectonics (Hamilton, 1998;Stern, 2008). For example, blueschists and ultra-high pressure (UHP) metamorphic terranes, considered to form exclusively at subduction zones, first appear in the early Neoproterozoic, broadly concomitant with an increase in the occurrence of preserved ophiolites (Stern, 2008;Brown, 2007;Harlow et al., 2015). However, a Neoproterozoic start to plate tectonics would conflict with an abundance of geologic evidence thought to reflect modern-style horizontal tectonics well before this time, including the formation of the well-established supercontinent Rodinia (e.g. Hoffman, 1988;Li et al., 2008). Such conflicts may be resolved by incomplete preservation in the geologic record: in light of preservation bias, the first appearances of plate tectonic hallmarks supply only one-sided constraints -providing evidence of plate tectonic processes at a given date, without necessarily ruling out earlier operation. This effect of preservation bias is broadly applicable to a wide range of records. For example, Shirey (Shirey and Richardson, 2011) provide a global record of lithospheric inclusions in diamond, identifying eclogitic inclusions as old as 3.0 Ga, but only peridotitic inclusions prior to that time. Due to Os isotopic enrichment, the eclogitic inclusions provide evidence for subduction and underplating of basaltic crust (Richardson et al., 2001). However, the absence of observed eclogitic inclusions older than 3 Ga does not exclude subduction before this date; unlike depleted lithospheric peridotite, eclogite is denser than the underlying asthenosphere and thus subject to delamination -a process which is only expedited by higher mantle temperatures (Jull and Kelemen, 2001) which likely prevailed in the Archean.
Jadeitites, blueschists, and other very high pressure, low temperature metamorphic rocks are thermodynamically fragile -far out of thermodynamic equilibrium in any cratonic geotherm -and are difficult to successfully exhume and preserve even in modern arcs (Ernst, 2003); such assemblages may not have formed in ancient subduction zones with higher mantle temperatures and more Mg-rich oceanic crust (Palin and White, 2015). Ophiolites, meanwhile, remain recognizable after more extensive metamorphism and alteration, but are susceptible to erosion due to both their chemical weatherability and their emplacement by obduction atop existing continental crust. At the other extreme, geochemical subduction signatures and isotopic tracers of crustal recycling can be preserved in principle as far back in time as the earliest surviving rocks and minerals. If the pervasive diachroneity of these data is the result of preservation bias, we should then be led to conclude that the oldest preserved constraint is the most accurate.
Meanwhile, the geochemical constraints (Figs. 3, 6) are difficult to escape. The effect of water on the extent of calc-alkaline differentiation is well established, and the major elements by which this effect is quantified are not readily biased (e.g. Sisson and Grove, 1993;Zimmer et al., 2010). The effect of slab fluids to produce relative enrichment of fluid-mobile elements and depletion of non-mobile elements is similarly well established (Perfit et al., 1980;Tatsumi et al., 1986;McCulloch and Gamble, 1991;Pearce and Peate, 1995;Schmidt and Poli, 2014), with no clear mechanism available to mask the absence of such arc signatures from the geologic record in the absence of subduction. If anything, the progressive leaching of fluid-soluble elements such as Rb from older rocks would only enhance such an effect. Together, we are left to conclude that the global long-term average proportion of flux and dehydration melting has not changed systematically at any time in the preserved rock record, and that subduction of oceanic crust has likely been in operation since at least the earliest Archean.
Nonetheless, our results need not preclude the operation of novel, nonuniformitarian processes in the Archean -e.g., widespread intracrustal diapirism (Van kranendonk and Smithies, 2007) -as long as such processes are compatible with the constant proportion of flux and decompression melting preserved in the continental basaltic record. Some of the most intriguing features of Earth's igneous rock record are endemic to the Archean, including granite-greenstone terranes, trondhjemites and komatiites. While our results suggest that basalts are generated by largely the same processes today as they were in the in the Archean, rapid (<200 Myr) transitions in granitic trace element compositions 2.5 Ga ago (Keller and Schoene, 2012) appear to suggest a change in the process of crustal differentiation, which remains coincident with widespread atmospheric oxidation at the end of the Archean.

Implications of ongoing subduction throughout the preserved rock record
The observed TiO 2 -MgO and Na 2 O-MgO systematics (Fig. 3e-f) tentatively suggest a continual, progressive increase in the average depth of mantle melting over Earth history -opposite the expectation from secular mantle cooling in a mantle decompression melting scenario (Iwamori et al., 1995). In a plate tectonic flux melting regime, however, the depth of flux melting is generally thought to decrease with increasing mantle potential temperature due to dehydration of the slab at lower pressures if all else is held constant (Peacock, 1990;Schmidt and Poli, 2014). One complication to this model arises from the fact that the composition of subducting slabs has not been constant through Earth history. Recently, Palin and White (2015) have suggested that more Mg-rich Archean slabs contain more stable hydrous minerals that can carry more water to greater depths, acting counter to the effect of increased mantle potential temperature. Mantle lithosphere, which may control the depth of last magma equilibration with the convecting mantle, provides another potential complication. Given that preserved continental lithosphere is roughly coeval with its overlying crust (Pearson, 1999), and increases in thickness towards the Archean (Artemieva, 2006), one might expect deeper depths of last magma equilibration in the Archean -also contrary to the observed pressure trends. However, this complication applies in detail only to decompression melting below preexisting continental crust with an extant mantle lithosphere; the depth of last equilibration in modern arc settings appears to be set instead by crustal thickness (Plank and Langmuir, 1988), suggesting that a highly rigid mantle lithosphere is not yet present in arc settings. While there is much room for further research, the trend in melting pressure suggested by Fig. 3e-f may well be most consistent with a predominance of flux melting in the continental basaltic record, as is suggested by the consistently calc-alkaline (THI<1) average observed tholeiitic index in Fig. 6d.
Other features of the observed basaltic record, such as the rapid (<1 Gyr) shift in trace element ratios (La/Yb, etc.), are similarly internally consistent with a plate tectonic setting, given that mineral stability thresholds in an arc mantle wedge setting readily explain the observed trends (Fig. 8). While we have not attempted to quantify the extent of the apparent increase in average mantle melting pressure, one notable consequence would be to shift the site of flux melting from the garnet peridotite stability field towards the spinel peridotite stability field in early Earth history -though potentially to little effect given the apparent extent of mantle melting beyond the point of garnet and clinopyroxene exhaustion in the Archean (Figs. 4, 7).
Even within a plate tectonic regime, the constant average continental basaltic water content and slab fluid input since ∼4 Ga is notable, suggesting a roughly constant proportion of arc and nonarc magmatism in the preserved continental record. If accurate, these results demonstrate that a constant ratio of flux to decompression melting is a persistent feature of terrestrial plate tectonics on Ga timescales.

Conclusions
The integration of large, freely available geochemical datasets with computational statistical methods has allowed for the quantification of systematic changes in average basaltic composition over time (Figs. 1-6), including gradual (Fig. 1), abrupt (Fig. 5), and invariant ( Fig. 6) trends at the 4-Ga scale. Although it is not possible to demonstrate complete independence of these results from sampling and preservation bias, weighted bootstrap resampling provides a robust method to obtain maximally representative information from the heterogeneous rock record. Internal consistency of the resulting trends, along with observed variations contrary to expected weathering and preservation biases, provides some confirmation that the results likely reflect secular variation in underlying drivers of the Earth system. These trends represent a significant advance in our ability to quantify changes in the composition of Earth's crust over time, allowing us to build our understanding of geological rates, states, and processes back through Earth history from the known state of the present.
Due to the lack of any systematic variation in abundance ratio for element pairs with similar igneous partitioning behavior but contrasting mobility during alteration and weathering (Fig. 2,  5a-b), we consider that these trends reflect igneous variability attributable to systematic temporal changes in magmatic process driven by mantle thermal state. An igneous origin for observed basaltic major and trace element geochemical trends is consistent with the interpretation of Keller and Schoene (2012) that the smoothly decreasing compatible and increasing incompatible element trends of average basalt over the past 4 Ga (Fig. 1) reflect decreasing average extent of mantle melting driven by secular mantle cooling. This interpretation is supported by the distinct observed signature of decreased partial melting relative to increased magma differentiation (Fig. 2) and the lack of apparent variation in average basaltic magma water content since 4 Ga (Fig. 3b-c,  6). Moreover, when properly accounting for the effects of mineral stability thresholds on trace element partitioning and melt productivity, (Figs. 7), it is observed that even the most abrupt observed transitions in trace element geochemistry are consistent with smooth, continuous mantle cooling as observed in Fig. 8.
Consequently, we propose that the geochemistry of the preserved basaltic crustal record can be most parsimoniously explained by plate tectonic processes modified by secular mantle cooling. In this interpretation, rapid variations in the geochemical record at the end of the Archean are effectively explained by the discontinuous effects of mineral stability thresholds on the petrologic processes of magma generation and differentiation, specifically garnet and clinopyroxene exhaustion during mantle partial melting. In addition to the consistent reproduction of observed trace element geochemistry (Fig. 8), support for this model comes from its reconciliation of a discontinuous (at the 2-standard-error level) percent melt trend with a smooth and continuous mantle melting temperature trend as observed in Fig. 4. This model is reinforced by the invariance of geochemical slab fluid input signatures through geologic time (Fig. 6a-c, d) which indicate a roughly constant proportion of crust produced in subduction zones over the entirety of the preserved rock record. Together, this evidence favors a horizontal, largely uniformitarian model of Archean tectonics. Given the profound geologic consequences of a planetary transition from a stagnant-lid or vertical tectonic regime to horizontal plate tectonics, the invariance of slab fluid and flux melting signatures (Fig. 6) and the consistency of the basaltic record with simple secular mantle cooling (Fig. 8) both suggest that subduction and associated horizontal tectonics had begun on Earth by at least 3.8 Ga and have operated consistently since.

Author contributions
C.B.K. conducted the computations and prepared the manuscript. Both authors contributed to the interpretation of the results.

Data availability
The full dataset and all computational source code is available under an MIT open-source license at https :/ /github .com / brenhinkeller /StatisticalGeochemistry. Acknowledgements C.B.K. was supported by the U.S. Department of Energy Computational Science Graduate Fellowship Program under contract DE-FG02-97ER25308. Reviews by Kent Condie and an anonymous referee substantially improved the manuscript.

Appendix A. The weighted bootstrap resampling method
To obtain an optimally even spatiotemporal sample weighting, each sample is assigned a weight inversely proportional to its proximity to other samples in space and time. Subsequently, the data are resampled (with replacement) in accordance with this weight such that data points from sparsely represented regions of the crust are resampled more frequently, while those from densely sampled regions are selected less frequently. Meanwhile, analytical uncertainty is accurately represented at each resampling step by the addition of Gaussian noise with appropriate variance for each variable in the dataset. The result is a continuous temporal record of the composition of preserved, exposed crust, represented as mean and two-standard-error uncertainty of the mean for 100 Ma time intervals over the last 4 Ga (e.g., Supplementary Fig. 1). Spatial and temporal prior and posterior sample distributions for the mafic sample subset examined in this paper are shown in Supplementary Fig. 2.
This bootstrap resampling methodology was introduced and discussed by Keller and Schoene (2012), but several points may bear further emphasis. In particular, the proximity of any sample i to all other samples in the dataset is calculated as a proximity sum of the form: where n is the number of samples in the dataset, z is spatial location, t is sample age, and a and b are constants where a = 1.8 arc degrees and b = 38 Myr. This equation and the associated constants are derived as follows.
We wish to obtain a more uniform sample distribution in space and time by assigning a greater weight (or resampling probability) to any samples that are far away from other samples in space and time. However calculating sample weight on the basis a naive distance sum will obviously not suffice, as this will be dominated only by far-field samples. A sum of inverse distances may be more appropriate, but this will introduce a singularity in the case of samples with identical positions. Instead, we desire to sum over a variable that should equal 1 for two samples with identical position, and 0 for two samples infinitely far apart. Consequently, we define proximity for some position variable x as: 1 ( x/s) p + 1 where s and p are constants. Any positive real values for s and p will satisfy our requirement that proximity is one when x = 0 and approaches zero as x approaches infinity. From inspection, it is clear that greater values of p diminish the proximity value for sample pairs with x > s and increase the proximity value for sample pairs with x < s.
When combining proximities for two variables that do not share a scale (e.g., space and time), we have a choice of scale constant for each variable. For geologic time and distance on the surface of the Earth, we have a clear maximum scale for each: 180 arc degrees for space, and ∼3800 Myr for the preserved continental rock record; this helps us set the desired proportionality between a and b. We must also consider the relative importance of near-field versus far-field samples. In considering a single geologic sample, the proximity of that sample to any number of samples on the other side of the Earth should not greatly influence the resampling probability. Instead, we wish to emphasize spatial and temporal differences on the scale of the fundamental unit of the igneous rock record: a single pluton or a single magmatic episode. As higher values of p diminish the proximity value of far-field samples and increase that of near-field samples (relative to scales a and b), we would like to set p > 1 along with a and b near the scale of a pluton or magmatic event. Within these constraints, variations in the exact value of p have relatively little impact on the resultant temporal trends (Supplementary Fig. 3).
A pluton or batholith may have a spatial scale up to tens or hundreds of km and a lifetime up to ones or tens of million years.
Setting a = 1.8 arc degrees (200 km) and b = 3.8 Myr satisfies these constraints while maintaining the desired proportionality between spatial scale and temporal scale. However, different weights may be chosen if one wishes to study a process with a different fundamental scale.
Having now a full expression for the proximity norm of any two samples, we may calculate the proximity sum k i for each sample i relative to all other samples j = 1 to n. However, as defined, k i will vary as a function of the total number of samples in the dataset, as well as the choice of constants a, b, and p. We thus require a normalization constant c, which may be adjusted to set the mean or median resampling probability (typically ∼0.2: lower mean resampling probability would allow for greater improvements in sample evenness, but at the cost of sample richness). We then may calculate sample weight: where an infinite weight corresponds to a 100% probability of resampling, or equivalently: Given the resampling probability for each sample, the full dataset may be resampled (with replacement) as many times as desired, binned by the independent variable of choice (typically age, in this study), and means and two-standard-error uncertainties calculated for each age bin. Convergence may occur within as few as ∼50-100 iterations, but even the conservative default of ∼10,000 iterations is not computationally prohibitive. Notably, where the dependent variable is a ratio of two elemental analyses, as in (Keller and Schoene, 2012) we represent each element pair as a fraction a/(a + b) prior to binning (for each analysis with both elements a and b defined), and subsequently return the binned results to a/b ratio form; this fractional average avoids the numerical instability inherent in ratios where the denominator may approach zero. Finally, although geological age uncertainties may be correlated, for e.g., samples derived from the same pluton or outcrop, the binned results are conveniently found to be relatively insensitive to the presence or absence of even quite extreme correlations in independent variable uncertainty (4).
While such weighted resampling reduces sampling bias, it cannot eliminate it completely (e.g., rocks that were never sampled cannot be included in the dataset). Similarly, rocks that have been lost from the crustal record or are otherwise not presently exposed on Earth's surface are necessarily absent from our record. The resulting preservation and exposure bias cannot be corrected without quantitative knowledge of the relationship between preservation / exposure potential and composition. The displayed twostandard-error uncertainties throughout this paper show only how well we know the average composition for basaltic rocks of a given age assuming a random, unbiased sampling of the crust, and thus must be treated as minimum uncertainties. Consequently, sampling, preservation, and exposure bias should be considered when interpreting temporal geochemical trends.
These limitations notwithstanding, the quantification of average compositions for preserved, surficially exposed continental crust through time provides an invaluable window into the deeper workings of the Earth system due to the ability of averages to isolate meaningful signals of secular variation from the heterogeneous rock record. Reducing a geochemical variable to a single global average scalar necessarily discards information. However, this reduction in complexity mitigates the effect of many confounding processes (e.g. local variability in structural/tectonic setting, mantle heterogeneity, partitioning within the sample set, etc.) such that the resulting average value is a more representative estimator of underlying secular variation. This follows directly from the significance of the arithmetic mean as an analogue for physical homogenization, and can be considered to reflect the inverse relationship between sample size and signal/noise ratio. This above effect is crucial for our purposes, since temporal variations in crustal composition are typically dwarfed by the instantaneous geochemical variability of the igneous rock record at any given time. Even though our resampled averages are likely not perfectly representative, they are the most representative data available at present, and are observed to display properties consistent with first-principles expectations, e.g, for secular temporal variation in mantle potential temperature. Observed variations follow clear systematics as expected on the basis of elemental compatibility during mantle melting (Supplementary Fig. 5) and are broadly robust for a range of SiO 2 and MgO windows which may be considered to encompass the majority of mantle-derived magmas ( Supplementary Fig. 6).
When interpreting the obtained results, it is worth bearing in mind that neither mid-ocean ridge basalts nor mantle plume basalts are common in the continental record. Since we include only samples geographically originating from present-day continental crust, non-continental basalts such as mid-ocean ridge basalts and ocean island basalts are excluded from our record except where they have been incorporated into the continental crust by natural geologic processes such as obduction. This exclusion is fully intentional, given that we are attempting to understand in particular the formation of the continental crust. Consequently, no subsequent geochemical interpretations should be considered representative of non-continental basalts such as MORB or OIB, nor the source mantle thereof.
The results demonstrate that weighted bootstrap resampling of a sufficiently large global dataset reveals a consistent set of clear temporal trends in average basaltic composition over Earth history. These trends are resolvable only thanks to the statistical properties of large sample sets, demonstrated in the general case by the law of large numbers. Although individual basalt compositions may vary widely from the true population mean due to geological (and to a lesser degree analytical) sources of variance, the effect of this variance on the mean of the observed sample set diminishes in proportion to the number of independent data in the set; equivalently, the standard error of the observed mean decreases as the square root of the number of independent data. Consequently, while it is impractical to exclude all low-quality data from a large dataset, even highly uncertain data provide new information about the population mean as long as individual sample compositions are independent and accurate to within their stated uncertainty.
Since both geological and analytical uncertainty are the result of successive convolutions of many independent sources of error, the resulting distribution of the uncertainty will (following the Central Limit Theorem) will approach a Gaussian distribution as the number of sources of error increases. Moreover, even drawing from a highly non-Gaussian distribution of geochemical data, the probability distribution of the mean will also approach a Gaussian distribution in accordance with the central limit theorem, such that variance and standard deviation are appropriate parameters with which to characterize the variability of the mean (Keller and Schoene, 2012). The result is a statistically robust record of secular geochemical evolution on which to build interpretations of geologic process.
Additional discussion of the weighted bootstrap resampling method and the standard error of the mean can be found in the full methods of Keller and Schoene (2012). In particular, we note that any rate of change derived from our Monte Carlo method is a minimum: geochronological uncertainty acts to smooth any rapid discontinuities into a continual transition. This effect applies even in the case of a gap in the underlying data, as discussed by Keller and Schoene (2012). Nonetheless, we note that our coarse binning resolution of 100 Myr implies that we cannot conclusively distinguish an instantaneous discontinuity from a protracted transition occurring over as much as 200 Myr (i.e., due to a Nyquist frequency of 1/200 Myr −1 ).
All data and source code is available at https :/ /github .com / brenhinkeller /StatisticalGeochemistry under an MIT open-source license.

Appendix B. alphaMELTS calculations and trace element modeling
To obtain thermodynamic insight into the geochemistry and equilibrium mineralogy of mantle partial melting, we conduct a wide range of pMELTS simulations (Ghiorso et al., 2002;Asimow et al., 2004) via the alphaMELTS command-line batch interface (Smith and Asimow, 2005). Such pMELTS simulations are used to (1) asses the sensitivity of major element melting differentiation trends to changes in melting extent and melting pressure, (2) invert for average mantle melting extent from the observed whole-rock major element record, and (3) estimate modal mineralogy of the mantle residue during partial melting to inform trace element partitioning calculations.
Deconvolving fractional crystallization, partial melting, and melting pressure In order to differentiate between fractional crystallization and mantle partial melting in the continental magmatic record, we use a range of bivariate differentiation plots where the opposing effects of mantle melting extent and fractional crystallization extent are oblique rather than subparallel in cross-plot space (Fig. 3). To inform the schematic effect direction arrows shown in Fig. 3, pMELTS simulations are conducted and plotted for each diagram over a range of pressures and melting extents, using a Matlab script to automate interaction with the alphaMELTS command-line interface, as illustrated in Supplementary Fig. 7. Since each observed differentiation line shares a common global mantle potential temperature (and global average mantle melting extent), fractional crystallization dominates the differentiation trend within a single time interval, while changes over time perpendicular to this trend are driven by either changes in mantle melting extent or (in some cases) melting pressure.
Specifically, the effects of melting pressure and melting extent were quantified by pMELTS simulations for equilibrium batch melting of primitive (McDonough and Sun, 1995) mantle, with 0.1 to 1 wt.% H 2 O and pressures ranging from 1 to 3 GPa in 250 MPa increments, and with a fixed oxygen fugacity set by the fayalitemagnetite-quartz buffer. The results of changing melting extent or water content are typically straightforward, while changing melting pressure may have more indirect implications. For instance, increasing pressure both increases solidus temperature and increases melt MgO for a given melting extent, resulting in higher concentrations of mantle-incompatible elements such as Na or Ti at a given MgO content (or equivalently, higher MgO at a given Na 2 O).
These results are consistent a range with experimental observations. The observation of increasing melt MgO as a function of increasing pressure at constant melt fraction is well known (e.g., Klein and Langmuir, 1987). Since we do not attempt to isolate or reconstruct primitive unfractionated magma compositions, we cannot observe this effect directly, but instead observe it indirectly through incompatible element systematics at constant MgO. For instance, Herzberg and Zhang (1996) illustrate that, at constant mantle melting extent, increasing mantle melting pressure decreases Na 2 O by nearly a factor of three (e.g. from ∼5.7 to 2 wt.%), but increases melt MgO by more than a factor of five (from ∼5 to 27 wt.%). However, equilibrium lines for variation in Na 2 O with variable extent of melting or differentiation at constant pressure are much steeper in Na 2 O-MgO space, such that at a 2 GPa magma with 10 wt.% MgO has nearly twice the Na 2 O content at equilibrium of a 1 GPa magma at 10 wt.% MgO. Analogous effects should apply for any mantle-incompatible element that increases steeply with decreasing MgO.
Increased Na partitioning into clinopyroxene with increasing pressure introduces an additional complicating factor when considering the above pressure systematics; however, even at high pressure (∼6 GPa), D Cpx/Melt Na does not exceed one, and thus does not strongly sequester Na (Blundy et al., 1995). Indeed, at intermediate pressures appropriate to arc mantle melting (2-3 GPa) Na remains Instead, Na 2 O and TiO 2 show equivalent systematics in Fig. 3, consistent with the hypothesis that both Na and Ti are both acting only as passive incompatible markers recording the underlying increase in primary magma MgO with pressure.

Estimating mantle melting extent
Quantitative estimates of mantle melting extent are produced by inversion of the resampled basaltic whole rock dataset as discussed in Section 3.2.4. In this procedure, each basaltic major element whole rock analysis in each of the (10,000) weighted bootstrap sample sets is assigned a mantle melting extent by least-squares comparison to model mantle melt compositions for a range of melting extents from 0 to 100% melt at temperatures between the solidus and the liquidus. After assigning a best-fit mantle melting extent and temperature for each (resampled) analysis, these variables are reduced in the same manner as any of the other resampled geochemical variables, to yield a mean and standard error for 100 Myr age bins between 0 and 3900 Myr. The required model compositions are produced by pMELTS-mode al-phaMELTS simulations for primitive mantle melting under a range of model melting conditions and melting modes including isobaric batch melting and continuous fractional polybaric melting. Fig. 4 illustrates the results in the case of isobaric batch (equilibrium) melting of McDonough and Sun (1995) primitive mantle ("pyrolite") with oxygen fugacity fixed by the fayalite-magnetitequartz buffer and 0.15 wt.% initial H 2 O at 2.2, 2.7, and 3.2 GPa, calculated between the liquidus and the solidus in 10 • C increments. Since mantle melting extent is estimated for each individual (resampled) whole-rock analysis for each set of pMELTS results, each (resampled) whole-rock analysis is individually assigned a estimated melting extent and estimated melting temperature. While estimated melting extents and subsequent averages are expected to be relatively insensitive to variations in the melting model, the corresponding estimated melting temperatures are obviously highly model dependent. Consequently, to asses the robustness of these results to data quality and model assumptions, a wide range of sensitivity tests have been conducted.
While the results appear to fall within an a-priori reasonable range (as discussed in the text), their accuracy is in part dependent on the degree of correlation among the error terms introduced for each separate major-element record by each potential confounding igneous process (fractional crystallization, phenocryst accumulation, assimilation, slab fluid input, etc.). In the simplest case of fully uncorrelated errors for N major elements, error in the end result (e.g., mantle melting extent) would decrease by a factor of √ N. Fully correlated errors, which would arise if the effects of confounding processes such as fractional crystallization were fully parallel to the mantle melting path in N-dimensional major element space, would allow for no reduction in the confounding bias. In reality, each of the above confounding processes acts neither parallel nor perpendicular to the mantle melting path, with precise effects varying as a function of numerous geological variables for each sample. Consequently, the resulting mantle potential temperature and melting extent estimates should not be taken as absolute. Potential igneous biases are schematically illustrated in Fig. 4, and should be kept in mind throughout.
Firstly, to test sensitivity of the results to the presence of high-MgO lithologies such as komatiites or cumulates, the mantle melting extent calculations are repeated excluding all analyses with more than 18 wt.% MgO, shown in Supplementary  Fig. 8. Due to their exceptional petrology and general rarity outside the Archean, komatiites are sometimes considered separately from other mantle-derived magmas. While excluding komatiites from the mafic magma set substantially decreases average Archean MgO abundances, it has little impact on most other major elements trends such as CaO, Al 2 O 3 , TiO 2 , Na 2 O or K 2 O, and thus only modestly decreases estimated average potential temperatures and melting extents. While excluding komatiites would likely be prudent if attempting to estimate, e.g., Archean oceanic basalt compositions or mid-ocean ridge melting conditions, Archean komatiites nonetheless represent real mantle-derived magmas preserved in the continental record, and thus should not be excluded from our continental sample set. However, as shown in Supplementary  Fig. 8, even if we were to exclude komatiites, the impact on estimated mantle potential temperature and melting extent would be unexceptional.
Similarly, reasonable variations in the compositional (leastsquares) comparison procedure also have little impact on the results, as illustrated in Supplementary Fig. 9. While in many cases samples with lower residuals may be considered more reliable, this is not necessarily the case here, given the substantial nonuniqueness of model melting curves introduced by uncertainty regarding the P, T, H 2 O, etc. conditions of mantle melting (c.f. Supplementary Fig. 7). Given the cost to sample richness and evenness of excluding half the potential sample population, it is not clear that excluding high-residual samples should be expected to yield more accurate results. In either case, the results are not dramatically altered when high-residual samples are excluded ( Supplementary  Fig. 9).
More complicated considerations involve the choice of mantle melting style. While isobaric equilibrium melting provides a simple, computationally stable, and reversible model, real mantle melting is rarely isobaric, and only asymptotically approaches equilibrium. Since equilibrium melting is path-independent, polybaric equilibrium melting would be equivalent to isobaric melting at the pressure of last equilibration. In the case of fractional melting, however, polybaric and isobaric melting may yield distinct results.
As an alternative endmember model to (isobaric) equilibrium melting then, we may consider continuous polybaric fractional melting, where new increments of melt are removed and completely isolated from the solid once they exceed an abundance threshold of, e.g. 0.001 wt.%. Such an alternative model is explored in Supplementary Fig. 10.
Given the arc genesis expected for many of the basalts preserved in the continental record, the polybaric fractional melting model used in Supplementary Fig. 10 follows a P-T path of the form that might be experienced by a packet of upwelling mantle in a subduction zone mantle wedge, beginning at the hydrous solidus just above at slab at 3.3 GPa and increasing in temperature with decreasing pressure until 2 GPa (c.f. Grove et al. (2012)). Continuous melting fractional melting is enabled using the ALPHAMELTS_CONTINUOUS_MELTING option, with the AL-PHAMELTS_MASSIN environment variable set to 0.001 g for a 100 g system to establish a 0.001 wt.% melt extraction threshold. While the vertical structure of the arc mantle melting region (and thus melting pressure) is constrained by the geometry of the mantle wedge, variations in global mantle potential temperature will be reflected by variations in potential temperature at the core of the mantle wedge. Consequently, the melting simulation is repeated for 82 distinct P-T paths with maximum temperatures ranging from 1265 • C to 1675 • C in 5 • C increments.
As seen in Supplementary Fig. 10a, this radical change in melting model has little impact on the global average mantle melting extent estimated over the past 4 Gyr. However, the estimated mantle potential temperature inverted using for a continuous polybaric fractional melting model markedly exceed those calculated for the isobaric equilibrium case (Supplementary Fig. 10b). This latter deviation is attributed to tendency of fractional melting to relegate all water to the isolated melt after the first few fractional melting steps, leaving the residue dry and thus driving up the liquidus temperature, as seen in Supplementary Fig. 10c. Such an effect is not considered realistic, deriving from the arbitrarily efficient melt segregation of the continuous melting model. However, a more accurate full reaction-infiltration model is not yet feasible in alphaMELTS, and would more closely approach the equilibrium melting model.
Finally, to more empirically validate the mantle melting model, we have conduct a range of inversions for modern mantle melts including MORB, OIB, and arc, and continental rift basalts using both isobaric equilibrium melting and polybaric fractional melting models. As seen in Supplementary Fig. 11, variations in estimated melting extent between different melting models are easily resolvable, but are small in absolute terms relative to the variations observed over the past 4 Gyr (∼3 wt.% vs ∼30 wt.%). Model variation in absolute melting temperature is much more dramatic, given that mantle melting temperature is directly influenced by mantle water content and melting pressure. In mantle potential temperature terms, continuous fractional melting behaves almost as anhydrous melting at similar pressures, since in a strict fractional melting model any initial water is rapidly extracted in the earliest melting steps, and thus is unable to depress melting temperature at any but the very lowest melting extents.
The results of Supplementary Fig. 11 highlight the disconnect between mantle melting extent and mantle melting temperature that is introduced by variable P and H 2 O conditions of mantle melting. For instance, while modern plume basalts are generated at higher mantle potential temperatures than modern MOR basalts, the shallowest extent of mantle decompression is limited by impingement on the rigid mantle lithosphere in the case of OIBs, whereas decompression and melting may extend to the surface in modern MORBs. Consequently, the difference in mantle potential temperature between OIB and MORB source mantle does not translate to a corresponding difference in melting extent. Arbitrarily inverting MORB and OIB compositions at the same final equilibration pressures yields unrealistically similar mantle potential temperatures for both lithologies; the true mantle potential difference between MORB and OIB is only evident when each is inverted under geodynamically accurate P and H 2 O conditions (denoted with asterisks in Supplementary Fig. 11b).
Unsurprisingly, MORBs yield reasonable temperatures when inverted assuming shallow and dry final equilibration conditions, OIBs when assuming deep and dry final equilibration, and arcs when assuming deep and hydrous final mantle equilibration (Supplementary Fig. 11b). In the case of arcs, melting is enhanced by the input of slab fluids, but hampered by the chilled temperatures of the mantle wedge and the limited extent of decompression due to thick overlying arc crust and/or lithosphere. Together, these observations suggest that while our absolute mantle potential temperature estimates are unavoidably model-dependent, relative temporal trends in temperature are more robust to model variation than are absolute values. As a corollary, we may note that the absolute temperature estimates of Fig. 4 which assume generally deep, consistently hydrous mantle/melt equilibration will be most appropriate if our determination of predominant arc origin and consistent proportion of flux melting (c.f. Section 5.1) is correct.

Trace element modeling
To better understand the observed continental basaltic trace element record, we compare observed trace element trends over the past 4 Gyr to expected trace element trends calculated on the basis of secular mantle cooling. In this model, expected trace element compositions are calculated using GERM mineral/melt trace element partition coefficients along with modal mineral abundances for the mantle residue calculated using pMELTS for the case of equilibrium melting of primitive mantle as above. This melting extent -trace element abundance relationship may then be used to map from mantle melting extent (as determined in Fig. 4) to av-erage basaltic trace element evolution over the past 4 Gyr (as in Fig. 8).
As in previous simulations to determine average mantle melting extent, modal mineralogy of the mantle residue was determined by pMELTS-mode alphaMELTS simulations of isobaric batch melting at 2.7 GPa with a fayalite-magnetite-quartz oxygen fugacity buffer, calculated for 10 • C increments between the liquidus and the solidus. Residual mineralogy was correlated by temperature for each step, and normalized to yield modal mineralogy for the purposes of bulk solid/liquid partition coefficient calculation.
Since mineral/melt partition coefficients of most elements increase with increasing melt polymerization, GERM partition coefficient data were fit by bootstrap resampling and binned averaging (in log-space) as a function of melt SiO 2 , as in Keller et al. (2015). Subsequently, appropriate partition coefficients were applied for each mineral in the mantle residue for each melting step in accordance with the appropriate melt SiO 2 , a bulk solid/liquid partition coefficient calculated, and equilibrium partitioning applied via the partitioning equation i is the concentration of element i in the liquid at for at melt extent F given initial concentration C 0 i solid/liquid partition coefficient D s/l .
While mantle major element composition may be satisfactorily approximated as primitive mantle during mantle partial melting, the same is not true for trace elements in a subduction zone setting, where trace element budgets are substantially modified by the input of slab fluids. Consequently, the initial trace element budget is inverted by scaling the model mantle trace element curve against the observed trend. However, since this inverted mantle trace element budget is constant for all time, the model is only expected to produce a temporal trend in agreement with the observed trend if (1) the mantle mineralogy and trace element partitioning model accurately reproduces the trace element evolution of the mantle during mantle partial melting, and (2) the melting extents inverted from the major element record alone in Fig. 4 record a similar mantle melting history to the trace element record alone. The close agreement of model and observed trace element temporal trends (Fig. 8) suggests that both conditions are met.