Volcanic–plutonic parity and the differentiation of the continental crust

The continental crust is central to the biological and geological history of Earth. However, crustal heterogeneity has prevented a thorough geochemical comparison of its primary igneous building blocks—volcanic and plutonic rocks—and the processes by which they differentiate to felsic compositions. Our analysis of a comprehensive global data set of volcanic and plutonic whole-rock geochemistry shows that differentiation trends from primitive basaltic to felsic compositions for volcanic versus plutonic samples are generally indistinguishable in subduction-zone settings, but are divergent in continental rifts. Offsets in major- and trace-element differentiation patterns in rift settings suggest higher water content in plutonic magmas and reduced eruptibility of hydrous silicate magmas relative to dry rift volcanics. In both tectonic settings, our results indicate that fractional crystallization, rather than crustal melting, is predominantly responsible for the production of intermediate and felsic magmas, emphasizing the role of mafic cumulates as a residue of crustal differentiation.

The continental crust is central to the biological and geological history of Earth. However, crustal heterogeneity has prevented a thorough geochemical comparison of its primary igneous building blocks-volcanic and plutonic rocks-and the processes by which they differentiate to felsic compositions. Our analysis of a comprehensive global data set of volcanic and plutonic whole-rock geochemistry shows that differentiation trends from primitive basaltic to felsic compositions for volcanic versus plutonic samples are generally indistinguishable in subduction-zone settings, but are divergent in continental rifts. Offsets in major-and trace-element differentiation patterns in rift settings suggest higher water content in plutonic magmas and reduced eruptibility of hydrous silicate magmas relative to dry rift volcanics. In both tectonic settings, our results indicate that fractional crystallization, rather than crustal melting, is predominantly responsible for the production of intermediate and felsic magmas, emphasizing the role of mafic cumulates as a residue of crustal differentiation.
Earth is the only planet in our Solar System with a high-silica continental crust, which is complemented by a low-silica oceanic crust similar to that of the Moon, Mars and Venus 1,2 . The resulting bimodal topography 3 of oceans and continents is essential for regulating Earth's climate 4,5 and providing elemental nutrients to the biosphere 6 . Igneous geochemistry drives our understanding of the formation and evolution of the continental crust, along with its depleted complement, the mantle 2,7 . Detailed investigations of specific igneous systems have revealed that mantle-derived basaltic melts differentiate to higher silica contents through processes such as fractional crystallization, metamorphism, and remelting of crustal material. Though these processes are generally agreed to result in a geochemically-stratified lithosphere which increases in silica abundance from the crust-mantle boundary to the surface 2,8 , the compositional and geological heterogeneity of the continents-scaling from kilometres to micrometres-has made it difficult to assess the relative importance of specific endmember processes. Solving this problem requires an understanding of why magmas (1) ascend through the crust and erupt to form volcanic rocks, or (2) stall and freeze at depth as plutons. Some models suggest that magma composition and volatile content dictate eruptibility 9,10 , while others emphasize physical processes such as crystal/melt segregation in magma chambers 11 resulting in complementary volcanic and plutonic compositions. Alternatively, eruptibility may be controlled by thermal constraints whereby higher magma fluxes result in more magma reaching the surface 12 , in which case systematic geochemical differences between volcanic and plutonic rocks are not required. Here we test such models by taking advantage of increased availability of whole-rock compositional data to perform a comprehensive statistical analysis of the geochemical relationship between volcanic and plutonic rocks. We posit that this comparison provides insight not only into the predominant mechanisms of igneous differentiation and the origins of crustal stratification, but also-since magmas that reach the Earth's surface are more efficiently recycled back into the mantle via erosion and subduction-into crust-mantle geochemical evolution over billion-year timescales.

Quantifying volcanic-plutonic parity
To determine if continental plutonic and volcanic rocks are compositionally equivalent, we have applied computational statistical techniques to a data set of major-and trace-element geochemical analyses of 122,751 plutonic and 171,690 volcanic whole-rock samples of the continental crust, extracted from the EarthChem database 13 . After filtering spurious data (Methods), each data set was subjected to weighted bootstrap Monte Carlo analysis 14 , with sample weights inversely proportional to spatial sample density in order to obtain a uniform, maximally representative posterior sample distribution on the globe (Methods, Extended Data Fig. 1). This method maximizes the accuracy of statistical parameters such as the mean and median, and minimizes the effect of oversampling in regions that have been subjected to more intensive field study. Sampling bias is further minimized in our interpretations by focusing on intensive compositional variables, binned by differentiation proxies such as SiO 2 and MgO (Methods, Figs 1,2). Results are reported as the mean and 2 standard error (s.e.) of the mean for 2 weight percent (wt%) SiO 2 bins.
The geochemistry of igneous rocks is known to correlate with tectonic setting 15 , arising in part from different styles of mantle melting: water-driven flux melting in subduction zones (volcanic arcs) versus decompression melting in rifts [16][17][18][19] . These contrasting processes form the basis of the classically defined calc-alkaline versus tholeiitic differentiation series, as observed in both major-and traceelement systematics 16,[19][20][21] . Since the proportion of volcanic versus plutonic magmatism may vary with tectonic setting, we have examined differentiation trends separately for arc and rift samples (based on geographical location for samples predominantly less than 200 Myr old; Methods, Extended Data Figs 2, 3) to prevent such an effect from biasing our geochemical analysis. The results (Figs 1, 2) are shown along with mid-ocean ridge (MOR) differentiation trends as an endmember tholeiitic case for comparison. Samples whose original tectonic environment could not be discerned (grey in Extended Data Fig. 2b) were excluded from plots of arc versus rift geochemistry (Figs 1-3), but included in subsequent differentiation figures.
Major and trace elements from both tectonic settings follow curvilinear differentiation trends (Figs 1, 2), with rift volcanics falling closer to the tholeiitic MOR data set. Elemental abundances consistently diverge for compositions more mafic than average primitive basalt (that is, ,50% SiO 2 ), with plutonic rocks enriched in MgO and depleted in nominally incompatible major elements (Fig. 1), consistent with the exposure of mafic plutonic cumulates. Enrichments in volcanic TiO 2 and K 2 O at low silica may also suggest the influence of highly eruptible peralkaline volcanics with few intrusive equivalents. Above 55% SiO 2 , elemental abundances for arc plutonic and volcanic uncertainties of the mean in 2% SiO 2 bins for plutonic (blue) and volcanic (red) samples from arc (a-g; left column) and rift (h-n; right column) tectonic settings. Mean and 2 s.e. confidence intervals for average MORB and felsic MOR differentiates are shown in grey for comparison. Throughout the figures, major elements (given as oxides) are reported as wt%, while trace elements (given as element only) are reported as p.p.m. Europium anomaly Eu/Eu* represents depletion (Eu/Eu* , 1) or enrichment (Eu/Eu* . 1) of Eu 21 relative to adjacent trivalent rare earth elements.

RESEARCH ARTICLE
rocks follow overlapping trends at the 2 s.e. level, while rift volcanics reveal subtle enrichments in FeO and depletions in CaO and Al 2 O 3 relative to rift plutonic samples. More dramatic differences between rift volcanic and plutonic rocks are observed in the trace-element data, with offsets by up to a factor of two across a wide range of SiO 2 . Ba and Sr are significantly enriched in plutonic relative to volcanic rift rocks above 50% SiO 2 (Fig. 2), while conversely, felsic rift plutonic rocks show depletions in Zr, Hf, and the heavy rare earth elements (HREEs). Divergences in rift samples remain when sampling only rocks with known ages ,100 Myr (Extended Data Fig. 5), such that plutonic compositions are not biased towards older, basement continental crust. In addition, while rift samples are generally associated with high heat fluxes 22 that may permit partial melting of pre-existing crust to produce incompatible-element-enriched melts, average rift differentiation trends in Fig. 1 do not show a strong tendency towards S-type 23 signatures, for example, with lower Al at high silica than in arcs. As Zr and Hf are nominally incompatible in all major rock-forming minerals, the inflection towards decreasing Zr and Hf with increasing silica reflects saturation of zircon (ZrSiO 4 ) 24 . Yttrium and the HREEs are similarly influenced by zircon fractionation, as well as by fractionation of amphibole and garnet 25,26 . The large divalent cations Ba and Sr, meanwhile, are primarily compatible in feldspar, particularly high-Na plagioclase 26 . Consequently, decreased feldspar fractionation and increased zircon fractionation in rift plutonic relative to rift volcanic magmas could produce the observed offsets in Zr, Hf, HREEs, and Sr at a given SiO 2 , consistent with a less dramatic but systematic depletion of CaO and Al 2 O 3 in rift volcanics compared to plutonics (Fig. 1).

The role of water in magma eruptibility
In one scenario for mineral fractionation between volcanic and plutonic rocks, eruptible felsic melt is extracted from a crystal-rich magma shortly before eruption, leaving behind a lower-SiO 2 , cumulate-like plutonic residue 11 . This mechanism may be expected to enrich the plutonic residue, as is observed, in Ba and Sr due to feldspar accumulation. However, it is not obvious why such a process would be active in rift environments but completely absent in arcs.
Differences between rift and arc magmas indicate a higher proportion of anhydrous decompression melting in rifts compared to dominantly fluid driven flux melting beneath arcs [17][18][19] . As seen in Fig. 1, rift volcanics consistently plot closer to the dry, tholeiitic MOR differentiation trend than other continental igneous rocks, with increased FeO* and TiO 2 and decreased CaO and Al 2 O 3 at a given silica content 19,21 (FeO* is defined in Fig. 3 legend). The tholeiitic character of rift volcanics relative to all other continental samples is illustrated most prominently in Fig. 3, in which rift volcanic rocks show initial FeO enrichment with decreasing MgO that typifies the tholeiitic differentiation trend 19,21 . In contrast, rift plutonic and arc rocks show lower FeO* at a given MgO and no FeO enrichment with decreasing MgO. Given the strong correlation between water content and the degree of tholeiitic versus calc-alkaline differentiation 20,21,27 , this observation is consistent with a systematically more hydrous composition for plutonic rift relative to volcanic rift samples.
Dissolved water in silicate melts has the effect of lowering the liquidus and solidus temperatures such that most silicate minerals become saturated in the magma at lower temperatures 20,27-30 , with a particularly strong effect for sodic plagioclase 20,28 . In contrast, experiments reveal little to no effect of water content on the saturation temperature of minerals such as magnetite 31 and zircon 24 . Thus, at high water activity, crystallization of common silicates will be suppressed relative to these phases at a given temperature: zircon and magnetite are therefore expected to saturate at lower magma SiO 2 contents during hydrous fractionation. Increased fractional crystallization of magnetite relative to plagioclase and other silicates in hydrous magmas is responsible for FeO depletion in the calc-alkaline trend 20,21,27,31 . We argue that given the dominant control of zircon on magma Zr and Hf, and combined with decreased partitioning of Ba and Sr in anorthite compared to albite 26 , calc-alkaline differentiation also explains enrichment of Ba and Sr and depletion of Zr, Hf, and HREEs in the melt phase relative to differentiation at lower water contents.
If Ba and Sr contents are elevated in hydrous magmas owing to suppressed plagioclase stability, the Eu anomaly (Eu/Eu*; Fig. 2), typically viewed as an indicator of plagioclase fractionation during differentiation, may be expected to follow a less negative trend for hydrous differentiation. In apparent contradiction to the plagioclase fractionation hypothesis, the data sets show that Eu/Eu* is identical in all plutonic and volcanic rocks above 52% SiO 2 regardless of tectonic setting. Surprisingly, however, arc Eu/Eu* is observed to follow equally or even more negative differentiation trends than MOR samples (Fig. 2). The relative lack of negative Eu/Eu* in dry MORB tholeiites has been observed in other data compilations 32 , and may suggest that water availability in a differentiating magma has little effect on resulting Eu/Eu*. Indeed, while plagioclase is often credited as the primary driver of magma Eu/Eu*, a compilation of published mineral/melt partition coefficients 26 illustrates that, for instance, orthoclase may often have a higher Eu 21 /Eu 31 partition coefficient ratio than plagioclase (Extended Data Table 1). Consequently Eu/Eu* is not directly diagnostic of plagioclase fractionation, and the lack of a more pronounced Eu anomaly in MOR and rift volcanic samples does not rule out increased plagioclase crystallization in these dry lithologies relative to arc and plutonic samples if coupled with complementary H 2 O-induced changes in the fractionation of other Eu/Eu*-active minerals.
Given the high fluid solubility of Ba and Sr 16,33 , source enrichment by aqueous fluids provides an alternative origin of enhanced Ba and Sr concentrations in rift plutonics. Such metasomatic fluids may be released along with fluid-mobile elements such as Ba, Sr, Pb, and Rb 16,33 by rift-associated heating of mantle lithosphere or hydrous continental crust. However, enrichment in fluid-soluble elements cannot solely explain Ba and Sr offsets in rift plutonics given the lack of equivalent offsets in other fluid-soluble elements like Rb and Pb (Fig. 2), which are observed to be strongly elevated in known cases of fluid-mobile element enrichment in metasomatic veins and aureoles 16 , and in experimental measurements of aqueous fluids produced by dehydration of oceanic crust 33 . Melts enriched by such fluids would be driven towards calc-alkaline differentiation trends  Shown are mean and 2 s.e. uncertainties of FeO* (total iron normalized as Fe 21 oxide) and FeO*/MgO ratio as functions of MgO and SiO 2 abundance, respectively. a, FeO* binned as a function of MgO in 1 wt% intervals. Increasing FeO* during differentiation from 8% MgO to 4% MgO indicates a tholeiitic differentiation trend, while decreasing FeO* indicates a calc-alkaline trend. b, FeO*/MgO binned as a function of SiO 2 in 2 wt% intervals. Greater FeO*/MgO at a given silica content indicates more tholeiitic differentiation.

ARTICLE RESEARCH
by the increased water content, with suppressed plagioclase fractionation and increased zircon fractionation generating a majority of observed volcanic-plutonic offsets. Enrichment by metasomatic fluids may provide an initial water source for hydrous plutonic magmas in rifts, but this cannot yet be conclusively confirmed or rejected with the present data set. Examination of the hydrous and anhydrous solidi for granitic compositions predicts a plutonic fate for hydrated granitoid magmas. Given the negative slope of the water-saturated solidus, hydrous magmas approach the solidus and may freeze upon decompression due to devolatilization, while decompression of dry magmas moves them away from the solidus towards greater degrees of partial melt ( Fig. 4; Supplementary Discussion) 9,10,29,30 . Consequently, water content may correlate with emplacement depth in settings where magmas with highly variable water contents are observed. In arcs, where magmas are almost ubiquitously hydrated 34 and often water-saturated by upper crustal depths for felsic compositions, the probability of eruption may instead be dictated by variables such as magma flux, crustal stress and thermal state. These quantities have already been suggested to control eruption on the basis of geochronological, thermal, and structural considerations 12,35,36 and arguably must do so in arc settings given identical arc volcanic and plutonic magma compositions. Therefore, while dry tholeiitic magmas have been documented in arc settings 19 , our analysis indicates that such melts do not affect average arc differentiation trends.

Generation of felsic continental crust
Despite the critical role of the continental crust on Earth, the relative contribution of two endmember processes responsible for crustal differentiation-fractional crystallization of mantle-derived basaltic magma and partial melting of older crust-remains uncertain [37][38][39][40] . The distribution of silica in crustal rocks, both globally ( Fig. 5) and regionally, reveals a pronounced gap in the abundance of intermediate compositions, generally referred to as the 'Daly Gap' [37][38][39]41,42 . This bimodality initially suggests partial melting as the dominant mechanism of differentiation 37 : low-degree melting of mantle-derived basalt with ,50% SiO 2 can produce granodioritic and trondhjemitic melts with ,70% SiO 2 or more-corresponding to the two peaks at ,50 and ,70% SiO 2 in the observed silica distribution 27,43 . In this model, intermediate magmas are produced primarily by magma mixing between basaltic and granitic endmembers 44 , as suggested by the observation of crystal-rich andesites 42 . Comparatively, fractional crystallization is less clearly predisposed to produce a bimodal silica distribution. Mechanisms have nonetheless been suggested by which compositional gaps may result from the systematics of crystal fractionation 38 or the dynamics of liquid-crystal separation 39 .
Fractional crystallization and partial melting models may be distinguished in our data set, however, since the two models result in contrasting major-and trace-element differentiation trends at intermediate compositions. Given that intermediate magmas are produced by mixing between granitic and basaltic endmembers in the partial , while the high-silica peak predominates in plutonic samples 30 . Since the volcanic record is more influenced by preservation and exposure bias than the plutonic record, the spatiotemporally weighted record of ref. 14 is used in ref. 30 to obtain a better estimate of average time-integrated continental crust, despite the smaller size of this data set. melting model, differentiation trends in this model must follow straight mixing lines, as physical mixing between two compositional endmembers invariably produces straight lines in element-element space 45 . Fractional crystallization, in contrast, produces curved trends through continuous differentiation [46][47][48][49] . Consistent with a fractional crystallization model, the observed volcanic and plutonic differentiation trends follow distinctly curved paths, even at intermediate compositions (Figs 1, 2). These curved paths resemble liquid lines of descent for fractional crystallization, overlapping with experimental batch fractional crystallization results (Fig. 6a) [47][48][49] . Comparison of observed trends to the differentiation trends of primitive basalt subject to varying degrees of magma mixing suggest a maximum 40% contribution of magma mixing to the production of intermediate magmas relative to idealized continuous fractional crystallization (Fig. 6a inset), decreasing to almost no allowed mixing relative to the experimental batch results. Abundant field and petrological evidence for magma mixing 44,50 in turn suggests an underlying (pre-mixing) differentiation trend with greater curvature than the observed record, and thus a fractionation mechanism intermediate in efficiency (on average) between batch and continuous fractional crystallization.
One first-order prediction of the fractional crystallization model is the production of abundant crystal cumulates as a by-product of magma differentiation. Such cumulates will vary in mineralogy and composition as a function of the pressure, temperature, and extent of melt extraction 46 , but must be compositionally complementary to the extracted melt, resulting in divergent cumulate and melt compositions 51 . As a result, the presence of cumulates in the plutonic record will produce observable compositional offsets between volcanic and plutonic averages for any nonlinear differentiation trends (see, for example, Extended Data Fig. 4). Clear compositional divergence attributable to plutonic cumulates is observed below 50% SiO 2 in both arc and rift settings (for example, MgO, Eu/Eu*). However, compositional offsets between felsic (.62% SiO 2 ) volcanic and plutonic differentiation trends are absent in arcs (Figs 1, 2) and are parsimoniously explained in rifts by the petrological effects of water. This lack of a felsic cumulate signature suggests that any cumulates with more than 60% SiO 2 are either small in volume, minimal in degree of melt extraction, or extensively remobilized into the volcanic record-and thus of limited significance for crustal differentiation. Consequently, we propose that the production of even felsic magmas by fractional crystallization results in a predominantly mafic cumulate residue with ,62% SiO 2 .
To test this prediction, we conducted geochemical modelling in which cumulate compositions, volatile contents, and pressure-temperature (P-T) paths representative of average crustal differentiation are estimated by inversion of 1.36 million pMELTS 52,53 simulations to fit the observed major-element differentiation paths given in Figs 1 and 2 (Methods). As displayed in Fig. 6b-e, felsic magmas with up to 80% SiO 2 are produced in equilibrium with a cumulate of no more than 60% SiO 2 . In addition to anticipated outputs of MgO-enriched ultramafic olivine-clinopyroxene-garnet-spinel cumulates, MELTS modelling elucidates the origin of several previously unexplained features of our data set (for example, subtle plutonic K 2 O enrichment between 52% and 62% SiO 2 due to orthoclase accumulation; Fig. 6e) while reproducing observed differentiation trends.
The production of intermediate and felsic crust by fractional crystallization is therefore consistent with both the curvature of observed differentiation trends and the existence of mafic cumulates in the plutonic record. However, questions remain regarding crustal silica distributions in a fractional crystallization scenario. As shown in Fig. 5, the compositional gap in crustal igneous rocks is largely driven by the contrasting silica distributions of volcanic and plutonic rocks, with a high-silica peak in plutonic samples and a low-silica peak in volcanic samples. This relationship is the opposite of that predicted by proposed mechanisms of compositional gap formation by mineral-melt separation 38,39 , and may indicate that the global whole-rock compositional gap is distinct from that observed in some individual volcanic suites 39,41 . The global Daly Gap as reflected in exposed igneous rocks appears to originate instead from

ARTICLE RESEARCH
superimposing a mantle-derived primitive basaltic distribution with a peak at ,52% SiO 2 onto a differentiated plutonic silica distribution skewed towards high silica.
While mafic cumulates clearly exist in the plutonic record, the higher average silica of plutonic samples is significant in the context of differentiation by fractional crystallization, and suggests that the majority of the mafic residue from crustal differentiation is lost from or otherwise not represented in the plutonic record. This is one expression of a widely-acknowledged problem in crustal mass balance: namely, the production of continental crust with ,61% SiO 2 from basaltic mantle inputs with ,50% SiO 2 (refs 2, 54). This issue is not unique to a fractional crystallization origin of the felsic crust, but can be averted in some crustal melting models by mechanisms where the mafic residue never leaves the mantle, for example, slab melting 43 .
In a fractional crystallization model (and most crustal melting models), however, a mechanism by which cumulates are physically removed from the continental crust is required-for example, lower crustal delamination 54 .
Further implications for the evolution of the crust-mantle system arise from the observed differences in volcanic and plutonic differentiation in rift settings. Because volcanic rocks are more easily eroded and removed from the crust, the tendency towards higher water in rift plutonic rocks should bias the mean crustal composition towards calc-alkaline, thus downplaying the importance of intracontinental magmatism to crustal growth. Meanwhile, the weatherable volcanics are more readily recycled and incorporated into surface biogeochemical cycles-for instance, elevated phosphorus in mafic volcanics (Fig. 1) will necessarily contribute to increased phosphate availability over geological timescales. Conversely, the preferential sequestration of elements enriched in plutonic rocks in the continental crust, such as K, Ba, and Sr (Figs 1, 2), when compounded over billions of years, influences the crust-mantle elemental budget and radiogenic isotope mass balance, and provides additional constraint for quantitative models of crustal growth from the Archaean to the present.
Although our results emphasize the distinction between arc and rift magmatism, the curvature of differentiation trends underscores the importance of fractional crystallization in both environments. Given the exclusively Phanerozoic age of surveyed arc and rift samples, our results suggest that it is possible to produce a bulk andesitic crust by processes of fractional crystallization and subsequent loss of mafic cumulates in modern tectonic settings, consistent with production of the continental crust by tectonic processes similar to those now in operation.
Online Content Methods, along with any additional Extended Data display items and Source Data, are available in the online version of the paper; references unique to these sections appear only in the online paper. Supplementary Information is available in the online version of the paper.

METHODS
In order to elucidate the geochemical relationship between volcanic and plutonic rocks, we have prepared and analysed a data set of major-and trace-element compositional data for 171,690 volcanic and 122,751 plutonic whole-rock samples together with location and lithological information (Supplementary Data 1). Data were obtained from the freely-accessible EarthChem Portal, including samples from NavDat, Georoc PetDB, and the US Geological Survey 13 . Oceanic samples were extracted to calculate tholeiitic MOR differentiation trends, while continental samples were separated into volcanic and plutonic data sets based on reported lithology. Volcanic and plutonic samples were further classified by cross-correlating sample location with present-day tectonic environment, using a geospatially referenced tectonic setting map modified from the USGS Global Crustal Database 55 (Extended Data Fig. 2). As illustrated in Extended Data Fig. 2, only samples from active oceanic and continental arcs are included in the 'arc' category, while samples from active rifts along with some well-established failed rifts (aulacogens) and plume-associated continental large igneous provinces were included in the 'rift' category. Consequently, rift samples represent those whose parental basaltic magma was produced primarily by decompression melting of the mantle, while arc samples are characterized by the contributions of flux melting in the arc mantle wedge. Continental samples whose original tectonic environment could not be discerned (grey area in Extended Data Fig. 2) were excluded from plots of arc versus rift geochemistry, but included in subsequent figures.
In contrast to the data set of Keller and Schoene 14 , a substantial proportion (,40%) of the whole-rock samples in our larger volcanic-plutonic data set are not associated with age constraints. Tectonic setting assignments based on present-day location are, naturally, only valid for samples produced in their current setting, limiting the accuracy of assignments for samples older than a few hundred million years. If numerous, mis-assigned samples would probably blur the distinction between arc and rift subsets. However, based on the distribution of known sample dates (Extended Data Fig. 3), the proportion of older 'inherited' samples (.300 Ma) in modern arc and rift settings is minor in our data set (,10%). This dominantly Phanerozoic age distribution has the additional benefit of increasing the probability that the process of tectonics and crust formation has not changed within the period of our analysis. Accordingly, differentiation diagrams including only samples with known ages in the last 100 Myr (Extended Data Fig. 5) display the same trends and offsets as those including all available data, only with somewhat larger standard errors of the mean due to smaller sample size.
Owing to the large size of the data set, spurious data was filtered both manually and algorithmically (Supplementary Data 4). Although outlier removal may introduce bias by rejecting valid but exceptional data, the presence of extreme outliers (for example, several orders of magnitude farther from the median than the 25% CI and probably not representing real bulk rock compositions) necessitated a method of outlier rejection to produce optimally accurate distributions of crustal composition. Consequently, conservative outlier rejection bounds were set manually in severe cases. Data representing physically impossible conditions (for example, compositional normalization ?100%) was not considered meaningful, and thus deleted categorically.
Data were then subjected to a modified version of the statistical approach of Keller and Schoene 14 . Since geochronological constraints are not available for many samples in the volcanic-plutonic data set, resampling weights were based only on sample location, rather than location and age as in Keller and Schoene 14 , to obtain an optimally uniform posterior spatial distribution (Extended Data Fig. 1). Although elemental distributions in the data set are often not Gaussian (Extended Data Figs 6,  7), the arithmetic mean of a set of sample compositions nonetheless has unique physical significance, representing the bulk composition that would be produced by physically homogenizing that sample set. Further, in contrast to the distribution of the elemental data, the distribution of the mean for each element is nearly Gaussian, due to the effects of the central limit theorem. Consequently, Figs 1, 2 show the arithmetic mean of binned volcanic and plutonic samples, along with the two-sigma standard error of the mean; equivalent plots based instead on the median are shown in Extended Data Fig. 6, with congruent results.
As the single dominant oxide component of most igneous rocks, wt% SiO 2 is simpler to interpret and easier to correlate with observed rock compositions than most alternative differentiation indicators. Nonetheless, similar volcanic-plutonic trends to those in the Harker diagrams (Figs 1, 2) may be observed as well with wt% MgO as the independent variable (Extended Data Fig. 8), with decreasing MgO corresponding to increasing differentiation. MELTS simulations. To test our predictions regarding the influence of fractional crystallization on average differentiation trends, the characteristics of cumulates in the plutonic record, and the effect of water on differentiation, we have conducted trace-and major-element geochemical modelling using highperformance computational resources at the DOE's National Energy Resource Supercomputing Center and the Princeton Institute for Computational Science and Engineering. The massive parameter space and the nonlinear, discontinuous equations of igneous differentiation processes such as fractional crystallization lead to highly nonunique solutions which cannot be obtained analytically except in the simplest cases. Consequently, we have approached the problem through a bruteforce Monte Carlo approach, minimizing 1.36 million MELTS fractionation paths with differing starting compositions and P-T paths to fit the observed crustal differentiation trends. These simulations were conducted using the alphaMELTS 53 v1.4.1 command-line version of the MELTS and pMELTS thermodynamic modelling software 52 along with parallel scripting codes (Supplementary Data 4, Supplementary Data 2). As each forward MELTS simulation is independent, the computation is inherently highly parallel, with minimal communication requirements between processes. However, due to the design of MELTS and alphaMELTS, numerous small configuration, input, and output files must be written for each simulation, such that networked filesystem performance may be a limiting factor for scalability. The average primitive oceanic basalt of Kelemen 56 was used as a starting composition, with varying initial H 2 O, CO 2 contents and P-T conditions. Initial H 2 O and CO 2 contents were drawn from uniform distributions between 0-4% by weight for H 2 O and between 0-1% for CO 2 . Initial oxygen fugacity was equilibrated at the fayalite-magnetite-quartz buffer (FMQ) but subsequently allowed to evolve in response to mineral fractionation.
Pressure ranges of differentiation were set following the assumption that basaltic magmas must be generated in the mantle (conservatively, at depths greater than ,33 km) and solidify in the crust. Maximum pressure is limited in a typical arc setting by the depth of the subducted slab (,100 km) but was constrained slightly shallower in practice so as not to greatly exceed the calibration range of pMELTS. Consequently, initial and final pressures were drawn from uniform random distributions between 1-2.5 and 0-1 GPa, respectively, with initial temperature set by the primitive basalt liquidus and a final temperature of 700 uC. Pseudorandom P-T paths were then generated by selecting five points from a uniform distribution for both pressure and temperature between the previously generated minima and maxima for each simulation, subject to the constraint that subsequent points must be non-increasing in both pressure and temperature. Each P-T path was then constructed by linearly interpolating between these seven points. The production simulations were then run using alphaMELTS v1.2 with a fractionation threshold of 0.005 and a 10 uC cooling step size. To explore the effect of oxygen fugacity, an equivalent set of simulations was conducted with an initial oxygen fugacity one log unit above the FMQ buffer (FMQ 1 1). The initial volatile compositions, average volatile compositions, and P-T paths of the best-fitting 350 simulations out of 1.36 3 10 6 simulations for initial oxygen fugacities of FMQ and FMQ 1 1 are shown in Extended Data Fig. 8, along with a comparison of best-fitting major-element paths for these two different initial redox states. Results elsewhere throughout the text show the FMQ version.
Given the large compositional and physical parameter space sampled, some simulations will inevitably fall in regions where MELTS is poorly calibrated. In addition, saturation of hydrous mafic phases (largely amphibole) is probably underestimated in general, while at high silica MELTS may underestimate the saturation of quartz and overestimate that of potassium feldspar; rhyolite-MELTS, which addresses the latter two problems 57 (not available in an alphaMELTS scripted version at time of our simulations) may result in different inversion results. Consequently, our simulations should be reproduced with multiple versions of MELTS before assigning particular significance to second-order results such as the inverted optimal P-T paths or water contents shown in Extended Data Fig. 8. Nonetheless, the MELTS results do unequivocally demonstrate that fractionation of the calculated cumulate compositions can accurately reproduce the observed average major-element differentiation trends. Trace-element calculations. In order to additionally simulate trace-element fractionation, we calculated representative average mineral/melt partitioning using partition coefficients from the GERM partition coefficient database 26 (Supplementary Data 3). Owing to the strong compositional dependence of most trace-element partition coefficients (mineral/melt partition coefficients tend to increase with increasing magma silicate network strength), partition coefficients were binned and interpolated as a function of magma SiO 2 content, propagating errors through unweighted Monte Carlo bootstrap resampling (Supplementary Data 4).
Trace-element differentiation paths were calculated for ideal fractional crystallization using mineral phase proportions from the best-fit MELTS simulations. As seen in Extended Data Fig. 9, calculated fractionation paths conform to the expectation of increasing Sr, and to a lesser degree Ba with increasing magma water content. As seen in Extended Data Table 1, the compiled partition RESEARCH ARTICLE coefficient data reveals high K d Eu/Eu* in orthoclase as well as plagioclase feldspar, resulting in no clear correlation of magma water content with Eu anomaly during differentiation (Extended Data Fig. 9). Solidus temperature of hydrous and dry magmas. While the opposing slopes of the hydrous versus dry solidi are unambiguously determined from experimental data for a wide range of silicate systems 30,58 , connecting these solidi to the preferential stalling of hydrous magmas upon decompression also requires lower average temperatures for hydrous magmas. The assumption of colder temperatures in hydrous magmas is somewhat justified since these magmas inevitably derive from systems with depressed solidi and liquidi. As an example, let us consider an arbitrary silicate protolith to which a finite quantity of heat is added, sufficient to produce an extractable proportion of silicate magma over a geologically relevant temporal and spatial scale. Two factors suggest that, in this scenario, final magma temperature will be linked to the solidus and liquidus temperatures of the protolith.
First, for processes within the solid earth, rates of heat transport are competitive with those of heat addition. Once the solidus is reached in one region of the protolith, heat can diffuse and advect to adjacent regions of unmelted protolith, until the available heat is consumed as latent heat of fusion. This is demonstrated by the extreme rarity of terrestrial supersolidus magmas, which are thought to be produced only in cases of extremely rapid heat addition such as large bolide impacts 59 .
Second, the molar heat capacity C p (in J mol 21 K 21 ) of most silicates is much lower (by two to three orders of magnitude) than their molar latent heat of fusion DH f (in J mol 21 ) 60 . In this sense, a difference in solidus temperature is equivalent in terms of heat consumed to a comparatively small change in the extent of protolith melting. As long as the liquidus is not exceeded, then, the final magma temperature will be a consistent offset above the solidus over a range of solidus temperatures.
In other words, once the solidus is reached, the remaining available heat is primarily used to produce more magma, not raise magma temperature. We may then infer that the addition of heat to a hydrous (low-T solidus) versus dry (high-T solidus) protolith will result in a relatively similar extent of melting-and thus a similar position relative to the appropriate solidus-but with higher temperature for the dry magma due to its higher solidus temperature.
The same reasoning holds for the addition of water to an existing magma: the addition of water will promote further melting (of entrained crystals, if any, or of the surrounding wall rock), decreasing magma temperature as the melting reaction consumes the available excess enthalpy. Causality of the plutonic water correlation. If our observed correlation between plutonic emplacement and geochemical signatures of higher water content is not fortuitous, there are at least two plausible causal mechanisms we must consider, with opposite directions of causality: that plutonic magmas are hydrous because they are not degassed near the surface, or that their higher water content forces them to stall at depth.
Given that the solubility of water in silicate magma increases with increasing pressure 61 , there is an inherent bias towards higher water content in magmas differentiated and emplaced at higher pressure. However, our geochemical evidence for higher water in plutonic magmas requires increased water content during the period of differentiation, rather than that of final emplacement. The apparent dearth of cumulates in the upper crust (Figs 1, 2) may indicate that proportionately little differentiation occurs in shallow subvolcanic magma chambers, partially decoupling depth of emplacement from depth of differentiation. If, despite these considerations, greater emplacement depth does contribute significantly to higher water content in plutonic magmas during differentiation, this effect and the associated geochemical signatures of higher water content in plutonic magmas during differentiation should be observed equally in both arc and rift environments, which is not what our data show. In contrast, higher water content leading to an increased likelihood of plutonic emplacement invokes no such contradiction if arc magmas are (as expected) uniformly hydrous at crustal depths. Such considerations, combined with a thermodynamic expectation for the stalling of hydrous magmas as discussed in the section on solidus temperature lead us to favour the latter direction (plutonic because hydrous) of causality. Eu anomalies. Although negative Eu anomalies in magmas are typically attributed to plagioclase fractionation 62,63 , plagioclase is not the only mineral with substantial potential effect on Eu/Eu*. Instead, as observed in Extended Data Table 1, minerals such as orthoclase, sphene, apatite, allanite, clinopyroxene, and garnet may also substantially both the Eu budget and Eu/Eu* ratio 26 (Supplementary Data 3). One potential explanation for the apparent lack of a systematic difference in Eu anomaly between dry and hydrous magmas hinges on the relative abundance of orthoclase and plagioclase feldspar. While water suppresses plagioclase stability (other factors held constant), it also increases the relative proportion-and possibly the absolute proportion-of orthoclase feldspar. Given the substantially greater K d Eu/Eu* and K d Eu of orthoclase relative to anorthite (Extended Data Table 1), the increased partitioning of Eu into orthoclase may offset the decrease in Eu fractionation by plagioclase, resulting in an equally negative Eu anomaly of the magma even in the case of decreased plagioclase fractionation. Such an explanation is not without caveats: for instance, considering only the feldspars, this mechanism would predict that any waterrelated Ba offset produced by plagioclase suppression would be similarly negated by orthoclase, which is not observed. However, the lack of any difference in Eu anomaly between arc and rift (including MOR) samples for a given silica range demonstrates that, whatever the its cumulative effects on mineral stability, water activity evidently does not substantially alter magma Eu anomalies (Fig. 2). Melt extraction from crystal-rich magmas. One potentially significant mechanism for the production of compositional offsets between volcanic and plutonic compositions is the extraction of mobile (and potentially eruptible) melt from a magma, leaving a cumulate plutonic residue 11,51 . At low silica, plutonic cumulates would be expected to be enriched in minerals such as olivine and pyroxene relative to non-cumulate volcanic magmas, resulting in enrichments of compatible elements such as Mg and Ni. Meanwhile, accumulation of plagioclase and pyroxene, as in (similarly mafic) anorthosites and cumulate gabbros would be expected to result in positive plutonic Eu anomalies. Evidence of such mafic cumulate processes is evident in the geological record 64-66 as well as from our analysis of geochemistry in both arc and rift environments (Figs 1, 2).
At higher silica, a similar process of melt extraction would be expected to leave plutonic cumulates enriched in feldspar along with potentially lesser proportions of amphibole, biotite, and quartz. This feldspar accumulation would enrich the felsic plutonic residue in Ba and Sr 51 , as is observed in rift settings. However, felsic plutonic cumulates resulting from shallow melt extraction would not readily account for the occurrence of plutonic Ba and Sr enrichments only in rift (and not arc) settings, nor the additional volcanic enrichments of Zr and Hf in rift settings. While plagioclase saturation may be suppressed in arcs relative to rifts due to higher overall water content, plagioclase remains a dominant rock-forming mineral in arcs. Accordingly, suppression of plagioclase in arc settings is evidently not so dramatic as to prevent or even markedly reduce the formation of anorthosite or leucogabbro-like mafic cumulates with positive Eu anomalies in arcs (Fig. 2), and would not be expected to prevent or reduce the formation of feldspathic-cumulates at higher silica when magma temperatures have fallen farther into the feldspar stability range. As such, melt extraction in should be expected to result in Ba and Sr enrichment of the residue in arc settings as much as in rifts. Consequently, although shallow volcanic melt extraction leaving a plutonic residue has been argued to occur in some scenarios 67 (and our results should not be taken to indicate otherwise), it does not appear to be a volumetrically significant mechanism for the production and differentiation of felsic crust.  (FMQ 1 1, right). a-f, Initial and average volatile distributions of 350 best-fitting simulations (a, b); P-T and P-SiO 2 paths of 350 best-fitting MELTS simulations (c, d); and major-element melt and cumulate compositions for 200 best-fitting MELTS simulations (e, f).

ARTICLE RESEARCH
Extended Data Figure 9 | Trace-element modelling based on MELTS simulation results. The figure shows simulated Sr (a), Ba (b) and Eu (c) differentiation paths (blue lines) and corresponding cumulate compositions (green points) calculated by applying GERM partition coefficients 26 to the mineral percentages derived from 200 MELTS simulations with best-fitting major-element trends. As in Extended Data Fig. 8, best-fitting MELTS simulations were selected out of a total of 1.36 3 10 6 MELTS simulations with varying initial H 2 O, CO 2 , and P-T paths (Methods) and initial oxygen fugacity set by the FMQ buffer. The diameter of each cumulate point is proportional to the volume of cumulates produced at each simulation step.

RESEARCH ARTICLE
G2015 Macmillan Publishers Limited. All rights reserved