TY - JOUR
AU1 - Grose, Christopher, J
AU2 - Afonso, Juan, C
AB - Abstract We examine REE (Rare-Earth Element) and isotopic (Sr–Hf–Nd–Pb) signatures in OIB (Ocean Island Basalts) as a function of lithospheric thickness and show that the data can be divided into thin- (<12 Ma) and thick-plate (>12 Ma) sub-sets. Comparison to geophysically constrained thermal plate models indicates that the demarcation age (∼12 Ma) corresponds to a lithospheric thickness of about 50 km. Thick-plate OIB show incompatible element and isotopic enrichments, whereas thin-plate lavas show MORB-like or slightly enriched values. We argue that enriched signatures in thick-plate OIB originate from low-degree melting at depths below the dry solidus, while depleted signatures in MORB and thin-plate OIB are indicative of higher-degree melting. We tested quantitative explanations of REE systematics using melting models for homogeneous fertile peridotite. Using experimental partition coefficients for major upper mantle minerals, our equilibrium melting models are not able to explain the data. However, using a new grain-scale disequilibrium melting model for the same homogeneous lithology the data can be explained. Disequilibrium models are able to explain the data by reducing the amount of incompatible element partitioning into low degree melts. To explore new levels of detail in disequilibrium phenomena, we employ the Monte-Carlo Potts model to characterize the textural evolution of a microstructure undergoing coarsening and phase transformation processes simultaneous with the diffusive partitioning of trace elements among solid phases and melt in decompressing mantle. We further employ inverse methods to study the thermochemical properties required for models to explain the OIB data. Both data and theory show that OIB erupted on spreading ridges contain signatures close to MORB values, although E-MORB provides the best fit. This indicates that MORB and OIB are produced by compositionally indistinguishable sources, although the isotopic data indicate that the source is heterogeneous. Also, a posteriori distributions are found for the temperature of the thermomechanical lithosphere-asthenosphere boundary (TLAB), the temperature in the source of OIB (Tp, oib) and the extent of equilibrium during melting (i.e. grain size). TLAB has been constrained to 1200–1300°C and Tp, oib is constrained to be <1400°C. However, we consider the constraints on Tp, oib as a description of all OIB to be provisional, because it is a statistical inference from the global dataset. Exceptional islands or island groups may exist, such as the classical ‘hotspots’ (Hawaii, Reunion, etc) and these islands may originate from hot sources. On the other hand, by the same statistical arguments their origins may be anomalously hydrated or enriched instead. Mean grain size in the source of OIB is about 1–5 mm, although this is also provisional due to a strong dependence on knowledge of partition coefficients, ascent rate and the melting function. We also perform an inversion in which partition coefficients were allowed to vary from their experimental values. In these inversions TLAB and Tp, oib are unchanged, but realizations close to equilibrium can be found when partition coefficients differ substantially from their experimental values. We also investigated bulk compositions in the source of OIB constrained by our inverse models. Corrections for crystallization effects provided ambiguous confirmations of previously proposed mantle compositions, with depleted mantle providing the poorest fits. We did not include isotopes in our models, but we briefly evaluate the lithospheric thickness effect on isotopes. Although REE data do not require a lithologically heterogeneous source, isotopes indicate that a minor enriched component disproportionately contributes to thick-plate OIB, but is diluted by high-degree melting in the generation of thin-plate OIB and MORB. INTRODUCTION The oceanic crust consists of geochemically and geophysically distinct regions: a roughly uniform layer of mid-ocean ridge basalt (MORB), pockmarked by ocean-island basalts (OIB). MORB originates from partial melting of upwelling mantle at mid-ocean ridges, while OIB are thought to originate from deeper mantle sources by an active upwelling impinging on the lithosphere (Morgan, 1972; McKenzie & Bickle, 1988). The characteristics of these two settings should, therefore, provide important constraints on the processes of mantle convection and the properties of the mantle. One such characteristic is the progressive enrichment of incompatible elements in OIB compared to MORB (Sun & McDonough, 1989; Hofmann, 2006; Niu et al., 2012). Although it is known that elevated incompatible element concentrations can be attributed to low extents of melting, it is widely believed that melting phenomena alone cannot explain the major features of the data (Gast, 1968; Hirschmann et al., 2003; Hofmann, 2006; Armienti & Gasperini, 2007; Humphreys & Niu, 2009; Niu et al., 2012). Thus, significant chemical variations observed among and between OIB and MORB have been attributed to compositional heterogeneity of the mantle (e.g. Stracke et al., 2003; Willbold & Stracke, 2006). In addition, an important characteristic of the data encompasses correlations between the chemistry of intra-plate lavas and the thickness of the oceanic lithosphere on which they erupt. As the OIB source ascends in the mantle, it must eventually encounter the base of the lithosphere. Beneath spreading centers the source may ascend to shallow levels, but elsewhere the final depth of melting is mechanically limited by the thickness of cool rigid lithosphere (Fig. 1). Correlations between lithospheric thickness and chemical signatures have been documented for trace and major elements as well as isotopes (e.g. Park, 1990; McKenzie & O’Nions, 1991; Ellam, 1992; Haase, 1996; Ito & Mahoney, 2005a, 2005b; Humphreys & Niu, 2009; Dasgupta et al., 2010; Niu et al., 2011). These correlations are important because, if the thickness of the oceanic lithosphere is known, then the final depth of melting can also be estimated. Thus, by sampling intra-plate lavas erupted on lithosphere for a range of seafloor ages, the compositional evolution of aggregate primary liquids may be estimated over a corresponding range of depths. Since seafloor age is generally well constrained (e.g. Müller et al., 2008) and the thermal structure of oceanic lithosphere has been well characterized from decades of modelling studies (e.g. Sclater & Francheteau, 1970; Davis & Lister, 1974; Stein & Stein, 1992; Doin & Fleitout, 1996; McKenzie et al., 2005; Afonso et al., 2008a, 2008b; Sleep, 2011; Grose, 2012; Grose & Afonso, 2013) with various approaches to the interpretation and analysis of global and regional geophysical datasets (e.g. Parsons & Sclater, 1977; Marty & Cazenave, 1989; Carlson & Johnson, 1994; Hillier & Watts, 2005; Crosby et al., 2006; Zhong et al., 2007; Korenaga & Korenaga, 2008; Crosby & McKenzie, 2009; Goutorbe & Hillier, 2013; Hasterok, 2013; Grose & Afonso, 2015), the lithospheric thickness and final depth of melting could be estimated with some confidence if the age of lavas which form the crust is also known. Fig. 1. View largeDownload slide Illustration of the cooling and thickening of the oceanic lithosphere with distance and age from the ridge axis. Isotherms are to scale based on model G13R developed by Grose & Afonso (2013) for mantle potential temperature Tp = 1300 °C. Also shown are plume-like conduits which ascend in the mantle until interrupted by mechanical interaction with the base of the lithosphere (LAB). The three small images illustrate examples of microstructure and microchemistry which our melting model predicts to exist in the ascending OIB source. Fig. 1. View largeDownload slide Illustration of the cooling and thickening of the oceanic lithosphere with distance and age from the ridge axis. Isotherms are to scale based on model G13R developed by Grose & Afonso (2013) for mantle potential temperature Tp = 1300 °C. Also shown are plume-like conduits which ascend in the mantle until interrupted by mechanical interaction with the base of the lithosphere (LAB). The three small images illustrate examples of microstructure and microchemistry which our melting model predicts to exist in the ascending OIB source. Trace element compositions in lavas and mantle rocks are typically interpreted with equilibrium melting models such as batch, near-fractional, or perfectly fractional melting (e.g. Shaw, 2006). In these models, melt instantaneously equilibrates with the solid residue as it is produced, but the segregation of melt from the source occurs differently. Batch melting assumes that no melt is extracted from the source until it reaches a final melting depth and, therefore, all melt remains in equilibrium with the solid at all times. Near-fractional (also referred to as dynamic, continuous, or critical) melting assumes that melt remains in equilibrium with the source until a critical porosity is reached, at which point the porosity is held constant and additional melting results in melt removal. Finally, perfect fractional melting assumes that the critical porosity is infinitesimal. Because of their simplicity these melt extraction models are widely employed in the interpretation of trace elements in melts and residues. The main limitation of these models is the equilibrium assumption, which requires that the rate of chemical exchange is much faster than the rate of melt production. Recognizing that volume diffusion rates for trace elements in mantle minerals may be too low to facilitate fast chemical exchange, previous studies have presented disequilibrium melting models which include diffusive exchange between grain interiors and surrounding melt (e.g. Qin, 1992; Spiegelman & Kenyon, 1992; Iwamori, 1992; 1993; Spiegelman & Elliot, 1993; Van Orman et al., 1998, 2002a, 2002b; Kogiso et al., 2004; Peate & Hawkesworth, 2005; Liang & Liu, 2016). Although our understanding of the role of disequilibrium in mantle melting processes was advanced by these previous studies, a significant caveat to previous models is that the grain geometry and phase distribution (i.e. microstructure) of a melting polycrystalline assemblage are not taken into account. Because the microstructure contains information such as the statistics of inter-phase grain contacts, grain geometry, variations in grain size distributions for different phases and complex microstructural changes associated with solid-state phase transformations and melting, a complete description of melting systems far from equilibrium requires a microstructural approach. In fact, a primary motivation of the present study is to work towards the development of a general model framework for non-equilibrium chemical thermodynamics within the framework of classical irreversible thermodynamics, which requires a grain-scale description of kinetics and energy. Since a true non-equilibrium thermodynamics approach adds additional layers of complexity, we use a numerically simpler model of microstructural evolution and focus on the predictions of trace element microchemistry during decompression melting of fertile mantle peridotite. In our model, diffusive partitioning of trace elements is coupled to a digital polycrystalline microstructure containing multiple solid phases. We employ the Monte-Carlo Potts model to simulate the microstructural evolution of grain boundaries, solid-state phase transformations and melting. Throughout the simulations, both the microstructure and the trace element microchemistry in the residue continuously evolve by non-equilibrium processes. We will explain the nature of the model and discuss general features of its predictions in melting systems. In the application of our model to mantle melting, we were motivated by our observation (discussed later) that while equilibrium melting models for realistic assemblages predict excessively enriched signatures in OIB, disequilibrium can suppress incompatible element enrichment to values which fit the data. An interesting consequence of this is that while the mantle is known to be isotopically heterogeneous, we find that enriched components are not required to explain major trends in Rare-Earth Elements (REE) in OIB. Also, by employing a Markov Chain Monte-Carlo inversion method we demonstrate that it is possible to constrain the extent of equilibrium in the source as well as important thermomechanical properties of the upper mantle, including the mantle potential temperature in the OIB source and the temperature at the thermomechanical lithosphere-asthenosphere boundary (LAB). REE SYSTEMATICS AND THE LITHOSPHERIC THICKNESS EFFECT OIB data compiled by Haase (1996) and Ito & Mahoney (2005a) showed that REE variations with seafloor age consist, generally, of MORB-like (if not slightly enriched) signatures in eruptions on thin lithosphere (<15 Ma) and incompatible element enriched signatures in lavas erupted on thicker lithosphere. More recently, Humphreys & Niu (2009) and Niu et al. (2011) used a larger dataset and argued that the variation of chemical signatures is more monotonic with age (or thickness). However, we have found that a reappraisal of seafloor ages using the global area–age grids of Müller et al. (2008) and geographic coordinates from Humphreys & Niu (2009) results in significantly different ages for their database of 115 islands (see Supplementary DataFigs S1 and S2, available for downloading at http://www.petrology.oxfordjournals.org). We plot their compiled data for La/Sm and Sm/Yb using our seafloor ages in Fig. 2. The basic features of the data agree with that recognized by Ito & Mahoney (2005a): MORB-like values in eruptions on young seafloor (<15 Ma) and higher values on older seafloor. We find no evidence of continued sensitivity of REE signatures to lithospheric thickness on seafloor older than about 15 Ma. Most lavas erupted on >12 Ma seafloor appear to be randomly distributed about a consistently high value (La/Sm∼5 and Sm/Yb∼4.5). We will refer to intraplate lavas erupted on <12 Ma seafloor as thin-plate OIB and lavas erupted on older seafloor as thick-plate OIB. Because the lithospheric thickness is a first-order control on REE signatures in lavas, characterizing OIB as a super-set including thin-plate and thick-plate variants should help provide constraints on thermochemical properties in the source as well as melting and crystallization processes which depend little on the final depth of melting. Fig. 2. View largeDownload slide Plots of (a) La/Sm and(b) Sm/Yb vs seafloor age, with no normalizations, in OIB, using the chemical data compiled by Humphreys & Niu (2009) and our seafloor ages for each island. The vertical error bars are standard deviations of the data for each island and horizontal error bars are age errors averaged over a 1°×1° box centered on the island’s geographic coordinates using the gridded seafloor area-age maps reported by Müller et al. (2008). Fig. 2. View largeDownload slide Plots of (a) La/Sm and(b) Sm/Yb vs seafloor age, with no normalizations, in OIB, using the chemical data compiled by Humphreys & Niu (2009) and our seafloor ages for each island. The vertical error bars are standard deviations of the data for each island and horizontal error bars are age errors averaged over a 1°×1° box centered on the island’s geographic coordinates using the gridded seafloor area-age maps reported by Müller et al. (2008). In our assessment, the demarcation age between these sub-sets likely originates from a marked jump in the rate of melt production (dF/dz, where F is extent of melting and z is depth) associated with the dry solidus. Melts accumulated below lithosphere thicker than the dry solidus depth primarily tap low-degree melts in the presence of highly incompatible elements such as volatiles (forming thick-plate OIB), whereas melts accumulated below thinner lithosphere progressively dilute the enriched signatures of these low-degree melts (forming thin-plate OIB). The thin- and thick-plate OIB subsets and their demarcation A statistic showing that the data do not follow a simple linear correlation with lithospheric thickness (or t1/2) is a comparison of slopes (dCk/dt1/2, where Ck is the concentration of element k) in best fit regression lines for the thick-plate data and for the full dataset. Results are shown in Table 1. Although the thin-plate dataset is only 20% of the total, slopes in the total dataset (for incompatible elements and Light/Heavy REE ratios) are about a factor of 2 to 3 higher than for the thick-plate only dataset. A more detailed view is provided in Supplementary DataFig. S3, available for downloading at http://www.petrology.oxfordjournals.org, showing the slope of the thick-plate OIB data for all possible values of the demarcation age. A jump in the data around 15 Ma is clear in both La/Sm and Sm/Yb, although it is strongest in Sm/Yb. We use 12 Ma as the demarcation age although the data may be best viewed as transitioning over a larger range (perhaps 9–16 Ma). Table 1: Slopes of regression lines through thin- and thick-plate OIB subsets and the full dataset for each element (ppm/My1/2) and ratio (My−1/2) Slope of Least Squares Regression line (TLAB = 1225°C), dCk/dt1/2 La/Sm Sm/Yb La Sm Yb Thin dataseta: 0·5194 0·2769 4·7974 1·1494 0·2425 Thick dataset: 0·0806 0·1244 1·0821 0·1201 0·0169 Full dataset: 0·2339 0·2601 3·078 0·3937 −0·0238 Full/thick datasetb: 2·902 2·091 2·845 3·278 −1·408 Slope of Least Squares Regression line (TLAB = 1225°C), dCk/dt1/2 La/Sm Sm/Yb La Sm Yb Thin dataseta: 0·5194 0·2769 4·7974 1·1494 0·2425 Thick dataset: 0·0806 0·1244 1·0821 0·1201 0·0169 Full dataset: 0·2339 0·2601 3·078 0·3937 −0·0238 Full/thick datasetb: 2·902 2·091 2·845 3·278 −1·408 a Regression of the thin-plate dataset is not interpreted but is shown for completeness. b The last row also shows the increase in slope of the full dataset compared to the thick-plate subset as a ratio. View Large Table 1: Slopes of regression lines through thin- and thick-plate OIB subsets and the full dataset for each element (ppm/My1/2) and ratio (My−1/2) Slope of Least Squares Regression line (TLAB = 1225°C), dCk/dt1/2 La/Sm Sm/Yb La Sm Yb Thin dataseta: 0·5194 0·2769 4·7974 1·1494 0·2425 Thick dataset: 0·0806 0·1244 1·0821 0·1201 0·0169 Full dataset: 0·2339 0·2601 3·078 0·3937 −0·0238 Full/thick datasetb: 2·902 2·091 2·845 3·278 −1·408 Slope of Least Squares Regression line (TLAB = 1225°C), dCk/dt1/2 La/Sm Sm/Yb La Sm Yb Thin dataseta: 0·5194 0·2769 4·7974 1·1494 0·2425 Thick dataset: 0·0806 0·1244 1·0821 0·1201 0·0169 Full dataset: 0·2339 0·2601 3·078 0·3937 −0·0238 Full/thick datasetb: 2·902 2·091 2·845 3·278 −1·408 a Regression of the thin-plate dataset is not interpreted but is shown for completeness. b The last row also shows the increase in slope of the full dataset compared to the thick-plate subset as a ratio. View Large In addition, a statistic showing that the thin- and thick-plate OIB are characterized by different signatures can be expressed by comparing their distributions with the similarity function ψk= ∫0∞Pk1Pk2dCk, (1) where Pk1 and Pk2 are normalized (integral is unity) density functions for the thin- and thick-plate sub-sets, respectively, for each element or ratio k. Distributions for these sub-sets for La/Sm and Sm/Yb are shown at the far right in each panel of Fig. 2 (also see Supplementary DataFig. S4, available for downloading at http://www.petrology.oxfordjournals.org). Similarity coefficients φk for these distributions are provided in Table 2. Similarity coefficients decrease for progressively more incompatible elements and Sm/Yb provides the lowest value as their distributions have the least overlap. Table 2: Mean values for thin- and thick-plate subsets of the OIB database, similarity coefficients for the distributions of each subset and the thick/thin ratio for each value Similaritya Thin-plate mean Thick-plate mean Thick/thin La/Sm 0·608 2·89 4·96 1·71 Sm/Yb 0·322 1·94 4·15 2·14 La 0·502 21·14 48·67 2·30 Sm 0·652 6·36 9·26 1·46 Yb 0·738 3·23 2·44 0·760 Similaritya Thin-plate mean Thick-plate mean Thick/thin La/Sm 0·608 2·89 4·96 1·71 Sm/Yb 0·322 1·94 4·15 2·14 La 0·502 21·14 48·67 2·30 Sm 0·652 6·36 9·26 1·46 Yb 0·738 3·23 2·44 0·760 a Similarity coefficients for thin-plate and thick-plate OIB distribution functions. View Large Table 2: Mean values for thin- and thick-plate subsets of the OIB database, similarity coefficients for the distributions of each subset and the thick/thin ratio for each value Similaritya Thin-plate mean Thick-plate mean Thick/thin La/Sm 0·608 2·89 4·96 1·71 Sm/Yb 0·322 1·94 4·15 2·14 La 0·502 21·14 48·67 2·30 Sm 0·652 6·36 9·26 1·46 Yb 0·738 3·23 2·44 0·760 Similaritya Thin-plate mean Thick-plate mean Thick/thin La/Sm 0·608 2·89 4·96 1·71 Sm/Yb 0·322 1·94 4·15 2·14 La 0·502 21·14 48·67 2·30 Sm 0·652 6·36 9·26 1·46 Yb 0·738 3·23 2·44 0·760 a Similarity coefficients for thin-plate and thick-plate OIB distribution functions. View Large Comparison of MORB and thick-plate OIB distributions One of the most insightful features of the data is seen in the compositional distributions of thick-plate OIB compared to MORB, shown in Fig. 3a and b. The distributions normalized to mean concentrations are also shown in Fig. 3 c and d. The purpose of comparing MORB to thick-plate OIB, instead of the full OIB dataset, is to isolate signatures in lavas with a roughly uniform degree of melting. Thus we remove OIB lavas which are likely to originate from some intermediate degree of melting. Notice that the normalized distributions for Sm, Yb and Sm/Yb in MORB and thick-plate OIB are remarkably similar. This similarity suggests that although the compatible to moderately incompatible elements are affected by processes which adjust all samples uniformly, processes which would cause skewing of the distributions are not different in these settings. For example, both the mean and distribution of Sm and Yb in the melting source must be similar in the sources of MORB and OIB. In contrast, La and La/Sm show skewing to lower values in MORB compared to thick-plate OIB. Although this could indicate a compositional effect of the source, this is likely a result of differences in shallow-level magmatic processes between these settings (O’Neill & Jenner, 2012; Coogan & O’Hara, 2015; O’Neill, 2016). Fig. 3. View largeDownload slide (a) Density functions for REE in MORB based on the dataset of Gale et al. (2013). Closed circles are the means of the distributions. (b) Density functions for REE in thick-plate OIB data based on the dataset of Humphreys & Niu (2009) with our corrected seafloor ages. (c) Density functions for REE in MORB normalized to the mean value. (d) Density functions for REE in thick-plate OIB normalized to the mean value. Fig. 3. View largeDownload slide (a) Density functions for REE in MORB based on the dataset of Gale et al. (2013). Closed circles are the means of the distributions. (b) Density functions for REE in thick-plate OIB data based on the dataset of Humphreys & Niu (2009) with our corrected seafloor ages. (c) Density functions for REE in MORB normalized to the mean value. (d) Density functions for REE in thick-plate OIB normalized to the mean value. Lithospheric thickness and the thermomechanical LAB Figure 1 shows isotherms from the thermal plate model of Grose & Afonso (2013) which we used to estimate lithospheric thickness. Using the 1100 °C and 1300 °C isotherms from this model the lithospheric thickness beneath each oceanic island is calculated and the chemical data are plotted in Fig. 4. Notice that for this amount of uncertainty in the temperature defining the thermomechanical LAB, the range of possible lithospheric thickness is large. For example, for TLAB = 1100 °C, plate thickness (as defined by this isotherm) at ∼12 Ma (the demarcation age for thin-plate and thick-plate OIB) is 36 km. However, for TLAB = 1300 °C, plate thickness is 63 km. So, a 200 °C increase in TLAB nearly doubles the estimate of lithospheric thickness at this age. Later we will show how this strong dependence of lithospheric thickness on the choice of TLAB allows TLAB to be precisely constrained. Fig. 4. View largeDownload slide Plots of La/Sm (top) and Sm/Yb (bottom) as a function of final melting depth (or lithospheric thickness) for LAB temperatures of 1100 °C (left) and 1300 °C (right). LAB temperatures are isotherms from the thermal plate model of Grose & Afonso (2013). Fig. 4. View largeDownload slide Plots of La/Sm (top) and Sm/Yb (bottom) as a function of final melting depth (or lithospheric thickness) for LAB temperatures of 1100 °C (left) and 1300 °C (right). LAB temperatures are isotherms from the thermal plate model of Grose & Afonso (2013). EQUILIBRIUM MELTING To test whether equilibrium melting can explain the data, we first test the predictions of equilibrium near-fractional and batch melting models. For both equilibrium and disequilibrium models evaluated in the present study we use three sources of information: (1) experimental partition coefficients for olivine (ol), orthopyroxene (opx), clinopyroxene (cpx) and garnet (grt) from Salters et al. (2002) and Adam & Green (2006); (2) a parameterization of melt productivity and major element composition in the residue of fertile peridotite based on the model of Herzberg (2004); and (3) calculations of phase abundances at thermodynamic equilibrium using the composition, temperature and pressure of the solid residue. The parameterizations for melt production and solid residue composition and the calculation of phase fractions, are described in Supplementary Data Appendix A, available for downloading at http://www.petrology.oxfordjournals.org. Total melt fraction as well as equilibrium assemblages in the residue are shown in Fig. 5 as a function of depth and mantle potential temperature. Overall, from great depth in the upper mantle, garnet gradually decomposes during decompression producing pyroxene phases which eventually begin to dissolve into the melt at a slow rate above the wet solidus depth and rapidly after the dry solidus (e.g. Solomatov & Reese, 2008). For partition coefficients, we prefer showing results using the values of Salters et al. (2002) due to the higher temperature and lower water content of their experimental system. However, because they did not report values for La we estimate them by extrapolation using the more complete dataset of Adam & Green (2006), ionic radius and Salters et al. (2002) values for Sm and Yb. Fig. 5. View largeDownload slide Phase fractions obtained from phase equilibrium calculations using the composition of the residue in our melting model for Tp = 1400 ± 100 °C as a function of depth. For line plots of the same data, see the Supplementary Data, available for downloading at http://www.petrology.oxfordjournals.org. Note that liquid fractions are parameterized (Supplementary Data Appendix A, available for downloading at http://www.petrology.oxfordjournals.org), rather than obtained from thermodynamic calculations. Fig. 5. View largeDownload slide Phase fractions obtained from phase equilibrium calculations using the composition of the residue in our melting model for Tp = 1400 ± 100 °C as a function of depth. For line plots of the same data, see the Supplementary Data, available for downloading at http://www.petrology.oxfordjournals.org. Note that liquid fractions are parameterized (Supplementary Data Appendix A, available for downloading at http://www.petrology.oxfordjournals.org), rather than obtained from thermodynamic calculations. Our predictions of La/Sm and Sm/Yb as a function of final melting depth for equilibrium batch melting and near-fractional melting are shown in Fig. 6. Equilibrium batch melting predictions are also shown using partition coefficients from both Salters et al. (2002) and Adam & Green (2006), whereas predictions for near-fractional melting use only Salters et al. (2002). For equilibrium near-fractional melting we show predictions for critical porosities between 0·05 and 1·0 %. All models assume that source potential temperature Tp is 1300 °C and values shown are normalized to the bulk source concentration (La/Sm = Sm/Yb = 1 in the source). Finally, we note that ‘potential concentrations’ of liquids are also shown for depths below the wet solidus (∼150 km in our models at this Tp). These are the concentrations in liquid predicted for an infinitesimal melt, although melting has not yet occurred. Fig. 6. View largeDownload slide (a) Equilibrium batch melting and (b) near-fractional melting models using our thermodynamically constrained equilibrium assemblage in the residue (Tp = 1300 °C) and the partition coefficient databases of Adam & Green (2006; dashed lines) and Salters et al. (2002; solid lines). The red lines are predictions of La/Sm and black lines are for Sm/Yb. Near-fractional melting models use the Salters et al. (2002) database only. Results for 20 near-fractional melting models are shown with the amount of residual melt varied between 0·05 and 1·0 %. The depths at which wet melting and dry melting occur in our models are indicated. On top of the predictions we also show the binned La/Sm and Sm/Yb data as connected filled circles. Note that model values are for a bulk source concentration of 1 (uncorrected for fractional crystallization). Thus, for a proper comparison of the data and models, one must imagine that the model values can be linearly shifted to higher or lower values to potentially fit the data. Fig. 6. View largeDownload slide (a) Equilibrium batch melting and (b) near-fractional melting models using our thermodynamically constrained equilibrium assemblage in the residue (Tp = 1300 °C) and the partition coefficient databases of Adam & Green (2006; dashed lines) and Salters et al. (2002; solid lines). The red lines are predictions of La/Sm and black lines are for Sm/Yb. Near-fractional melting models use the Salters et al. (2002) database only. Results for 20 near-fractional melting models are shown with the amount of residual melt varied between 0·05 and 1·0 %. The depths at which wet melting and dry melting occur in our models are indicated. On top of the predictions we also show the binned La/Sm and Sm/Yb data as connected filled circles. Note that model values are for a bulk source concentration of 1 (uncorrected for fractional crystallization). Thus, for a proper comparison of the data and models, one must imagine that the model values can be linearly shifted to higher or lower values to potentially fit the data. Comparing equilibrium model predictions (Fig. 6) with the data (Fig. 4) reveals significant similarities. High La/Sm and Sm/Yb values are predicted in lavas erupted on thick lithosphere (where the final melting depth is deep) and low values are predicted over thin lithosphere (where the final melting depth is shallow). Moreover, the predictions show that REE ratios as a function of final melting depth are non-linear and step-like, particularly for fractional melting. In the case of batch melting (Fig. 6a), Sm/Yb gradually increases with depth by a factor of 2–3 over the 0–60 km depth range, followed by a step-like increase to values up to 15 times that in the source at greater depth. La/Sm in the batch melting case behaves similarly, except that the gradual decrease over the 0–60 km depth range is smaller and instead of a pronounced step at the dry solidus depth (60 km), La/Sm gradually increases with depth. In the case of near-fractional melting (Fig. 6b) both La/Sm and Sm/Yb are more elevated in deep low-degree melts generated below the dry solidus. While the general agreement between equilibrium melting models and data is revealing, at least two important differences remain. Firstly, concentrations predicted in deep low degree melts are much higher than observed in the data. Secondly, although the observations show little to no variation in concentration with lithospheric thickness after a seafloor age of about 12 Ma, predicted values continue to increase with depth. These misfits may be attributed to one or more of the following: (1) incorrect partition coefficients; (2) lithological heterogeneity; (3) low abundance phases; (4) the form of the function for melt production with pressure; (5) reactive melt transport; or (6) chemical disequilibrium. Although we do not dismiss the role of other factors, here we focus mainly on the effect of disequilibrium and show that, even in the case of a homogeneous fertile peridotite, disequilibrium melting can fit the REE data without ad hoc adjustments to physical coefficients. DISEQUILIBRIUM MELTING MODEL In the following we describe the major characteristics of our disequilibrium melting model. Additional information can be found in Appendices A–C. In summary view, we model solid-state phase transformations and melting jointly with trace element diffusion within a 2 D microstructure as the rock decompresses in the mantle (Fig. 1). Explorative 3 D models were also investigated, but we found that the important dimensional and geometric aspects of 3 D microstructures are preserved in 2 D models (Supplementary Data, available for downloading at http://www.petrology.oxfordjournals.org). The diffusion model calculates the trace element composition of liquids in the microstructure and we predict lava compositions by gradually extracting and aggregating the melt according to near-fractional melt segregation. The initial phase fractions are obtained by Gibbs free-energy minimization (Fig. 5). Equilibrium thermodynamic calculations are also used to inform the evolution of phase transformations within the microstructure (as discussed later). The Monte-Carlo Potts model is used to build an initial multi-phase polycrystal to represent the initial condition of the melting model at depth (200–300 km) in the mantle (Fig. 7) and to evolve the microstructure during decompression. At the start of the melting model, trace element concentrations in the microstructure are distributed according to equilibrium mineral/melt partition coefficients at each lattice site. However, phase transformations change the local mineral/melt partition coefficients of the lattice sites where the transformations take place. Since trace element concentration is the same before and after transformation, the new state must be far from equilibrium, requiring diffusion throughout the entire system until a new equilibrium is achieved. Fig. 7. View largeDownload slide Illustration showing the main features of the Monte-Carlo Potts Model. The top row of panels show the phase distribution, the second row shows the spin distribution (where the ‘spin’ is a numerical id differentiating contiguous crystals), the third row shows the functional for the spatial distribution of interfacial energy and the bottom row shows an example of microchemistry for a trace element highly compatible in garnet and moderately compatible in cpx. Warm colours indicate high concentration, although quantitative details of the illustration are not discussed. The left column shows that the Potts model is initiated with a fine grained random distribution of phases (of known volume fractions) and maximum interfacial energy (interfacial energy is everywhere high). The middle column shows the result of running the Potts model (holding phase fractions constant) until grain sizes are large enough to serve as the initial condition for the melting model (where mean grain radius is about Lg/13 and Lg is the model edge length). Chemical components are initially distributed at equilibrium in the microstructure (constant values in same phase crystals). The right column of panels shows a typical result of running the Potts model with the phase transformations (solid-state and melting) predicted by decompression, concurrently with the model for diffusive partitioning of trace elements. Fig. 7. View largeDownload slide Illustration showing the main features of the Monte-Carlo Potts Model. The top row of panels show the phase distribution, the second row shows the spin distribution (where the ‘spin’ is a numerical id differentiating contiguous crystals), the third row shows the functional for the spatial distribution of interfacial energy and the bottom row shows an example of microchemistry for a trace element highly compatible in garnet and moderately compatible in cpx. Warm colours indicate high concentration, although quantitative details of the illustration are not discussed. The left column shows that the Potts model is initiated with a fine grained random distribution of phases (of known volume fractions) and maximum interfacial energy (interfacial energy is everywhere high). The middle column shows the result of running the Potts model (holding phase fractions constant) until grain sizes are large enough to serve as the initial condition for the melting model (where mean grain radius is about Lg/13 and Lg is the model edge length). Chemical components are initially distributed at equilibrium in the microstructure (constant values in same phase crystals). The right column of panels shows a typical result of running the Potts model with the phase transformations (solid-state and melting) predicted by decompression, concurrently with the model for diffusive partitioning of trace elements. Microstructural simulation with the Monte Carlo Potts model The main ingredient of our model is a numerical technique for time-dependent microstructural evolution (i.e. evolution of grain size, grain geometry and the distribution of phases). For this we employ the kinetic Monte-Carlo Potts model (hereafter, Potts model), a powerful computational approach to simulating microstructural phenomena in solids and liquids. The method has been widely used in computational metallurgy, ceramics and foam physics (cf., Janssens et al., 2007) and has occasionally been used in the context of geological systems (e.g. Solomatov et al., 2002; Piazolo & Jessell, 2008). Microstructures produced by the Potts model, including grain geometries, grain size distributions, phase distributions and dihedral angles between multi-phase triple junctions strongly resemble microstructures observed in solid and liquid-bearing solid materials (e.g. Anderson et al., 1989; Yang et al., 2000; Hui et al., 2003). Forward evolution in the Potts model is governed by the minimization of interfacial energy. Using an energy function to describe the energy density of a digitized microstructure and a probabilistic determination of local changes of state, forward simulation evolves a fine-grained high-energy system toward an ordered coarse low-energy state. This coarsening phenomenon is depicted in Fig. 7, in which the left column panels show the initial Potts model state and the middle column shows the result after a large number of time steps. Phases and grains are initially randomly distributed and the energy is high everywhere. After a number of state changes which lower system energy occur, large crystals form and high energy only occurs along grain boundaries. Models have periodic boundary conditions. The Potts model schemes presented below are basically similar to previously developed models (e.g. Tikare, 1995; Tikare et al., 1998; Janssens et al., 2007), but are significantly modified to encompass multiple phases, phase transformations and better resolved surfaces. The model begins by defining an N × M matrix where the elements represent a spatial distribution of lattice sites. All sites are then populated with properties such as the stable phase, P and the so-called spin, S (i.e. S defines a particular state of the lattice, see below). The phases are initially distributed randomly (Fig. 7), although with known modal abundances (e.g. 50% olivine, 10% opx, 5% melt, etc.). The lattice spin numerically differentiates individual grains and their boundaries and in a more physically explicative model it might be taken to represent the crystal orientation. Typically the number of possible spin states Q is a taken to be a small integer on the order of 102. In our formulation, the number of spin states is unlimited (infinite Q-state Potts model). The energy density Eα is computed from the energy function Eα=12∑βn=36eαβ, (2) where the subscripts α and β refer to the site and its neighbours, respectively, n is the number of interacting neighbours and eαβ is the energy of the interaction. The interaction energies depend on the spins and phases of the interacting sites such that eαβ=0ifeαβ=γssifS(α)≠Sβ,S(αs)≠Sβs,eαβ=γslifeαβ=γsl ifS(αl)≠Sβs,S(αs)≠Sβl, where S(α) is the spin of site α, subscripts s and l indicate that the sites are solid and liquid, respectively, γss is the solid-solid interaction energy and γsl is the solid-liquid interaction energy. The energy of the entire system is E=12∑αijN∑βn=36eαβ, (3) where N is the total number of lattice sites. While typical Potts models use five or nine point stencils (1st and 2nd nearest neighbours), here we use a much larger 37-point stencil with n = 36 nearest neighbouring elements (7th nearest neighbours). This stencil for the energy functional is shown in the small subpanels of Fig. 7. The main reasons for the large stencil are to inhibit lattice pinning effects (e.g. Holm et al., 1991) and to allow better resolution of boundary curvature (Mason et al., 2015), thus promoting realistic Zener pinning effects around interstitial grains. While Potts model microstructures normally generate jagged grain edges, a large energy stencil results in rounded features which appear similar to the petrography of mantle rocks and the predictions of phase-field models (compare our results to those of Tikare et al., 1998). The interaction energy employed in the Potts model is related to the extensive thermodynamic quantity of interfacial free energy, which is also related to the concept of surface tension (c.f., Porter & Easterling, 1981). In the presence of ≥2 phases whose boundaries are characterized by different free energies per unit area (or length in our 2 D models), the interfacial configuration of the system will evolve toward a structure which preferentially minimizes the surface area of interfaces with high free energy. Due to the low free energy of liquid-solid interfaces compared to solid-solid interfaces (γss/γsl=2·5 in our models), the system will increase the area of solid-liquid interfaces by drawing excess liquid away from triple junctions and into inter-crystalline boundaries. The resulting migration of liquid sites is a type of capillary action which causes wetting of grain boundaries by the liquid. It is also known that interfacial energy depends on misorientation angle, but due to limited data and for model simplicity we assume that interfacial energy depends only on the interacting phases. The Monte-Carlo procedure in the Potts model is as follows: Calculate energy density configuration using periodic boundary conditions. Randomly select a lattice site and one of its eight neighbouring sites. If the spin state of the site and its selected neighbour are equal, go back to step 2. Otherwise, continue to next step. If the phase of the site and its neighbour are unequal, continue to step 5. Otherwise, the site adopts the spin of the neighbour and continue to step 7. This step facilitates grain boundary migration processes. If the site is a solid phase, it adopts a new random spin (liquid sites always have equal spin). This step facilitates Ostwald ripening processes. The site and its neighbour exchange places. Spin and phase are exchanged. The new energy configuration is calculated. Use a probability transition function to determine whether to accept or reject the change. We use the common Metropolis function where the probability h of accepting a change is (Janssens et al., 2007) h=1 if ΔE≤0h=exp  (-ΔE/KT) if ΔE>0 (4) where ΔE is the energy change and KT is the thermal energy, or noise, of the system. The above algorithm is repeated once for all lattice sites. Once all lattice sites have been tested for state changes, the first Monte-Carlo step (MCS) has passed. For our models, the normal value of KT is taken to be 1·0 and the normal value for time scaling is 1 MCS = 1 year. Disequilibrium petrology The petrology of the mantle in disequilibrium melting models uses the equilibrium model previously described with features added to account for phase transformations in a 2 D microstructure. During decompression, pressure and temperature decrease, resulting in a different equilibrium phase assemblage. The microstructure must thus transform step-wise from the existing assemblage toward the equilibrium assemblage. A phase transformation model must describe how phase transformation takes place in the numerical lattice, where they may take place and the forces controlling the rates of phase transformation and which phases may transform into other phases. In our model, phase transformations occur by changing the phase ID of a lattice site, are only allowed to occur near the edges of grains (heterogeneous nucleation), the ‘force’ governing phase transformation is the difference between the existing and equilibrium phase fractions, the rate of transformation requires a rate coefficient in the probability equation and we require that phases may only appear from a restricted set of parent phases (e.g. olivine is not allowed to transform into olivine). Based on these rules, a probability that the lattice will transform into any different phase is calculated for every lattice site in the model and whether or not the transformation takes place is determined by sampling a transition function. The details of phase transformation are described in Supplementary Data Appendix B, available for downloading at http://www.petrology.oxfordjournals.org. Perhaps the most important simplification of our phase transformation model is the force controlling phase transformation. In our model the force is proportional to the deficit or excess in the abundance of a phase, whereas in real systems the force is the bulk free energy change associated with the transformation. Because this simplification compromises the kinetics of phase transformation and major element diffusion as a thermodynamically coupled process, we employ a rate coefficient which is high enough to ensure that the system is always close to phase equilibrium. In future work we will generalize our model for major element diffusion in a self-consistent non-equilibrium thermodynamics framework. Such improvements will be necessary to confirm the textural evolution observed in our models and to study the impact of phase transformation and major element kinetics on trace element disequilibria. A snapshot of a typical model near the garnet-out transition depth is shown in the right column of panels in Fig. 7. Because the growth of pyroxenes from garnet breakdown involves the growth of a new population of small interstitial grains, size distributions are bimodal (Supplementary Data Appendix K, available for downloading at http://www.petrology.oxfordjournals.org). Such unstable textures are characteristic of phase transformations and would not normally arise during normal coarsening processes. Equilibrium modal abundances of olivine, orthopyroxene, clinopyroxene and garnet computed as a function of total melt fraction and depth are shown in Fig. 5 for Tp = 1375 ± 75 °C. Our phase equilibrium calculations show good agreement with estimates of mantle petrology based on abyssal-peridotite data (Niu et al., 1997). Diffusive partitioning At the start of the melting model, we assume equilibrium distributions of trace elements (Fig. 7). Because phases are characterized by unique partition coefficients, phase transformations result in a state of local chemical disequilibrium, thereby requiring that diffusion occurs throughout the entire system until a new equilibrium is achieved. However, as diffusive exchange may be slow for a given rate of phase transformation and grain size, a highly complex microchemical heterogeneity typically arises (Fig. 7). Diffusive partitioning of trace elements is described with an explicit finite-difference approximation of the advection-diffusion equation ∂Ck∂x+v⋅∇Ck=∇⋅Dk∇C^k, (5) where Ck is the concentration of the trace component k, Dk is the diffusivity of the component and the hat on Ĉk indicates that the concentration gradients require a special treatment when there exists a discontinuous variation in the partition coefficients (interphase partitioning). Note that although we show an advection term, advection is treated exclusively by lattice site motions in the Potts model and tends to occur only along grain boundaries. Diffusion across phase boundaries (grain boundaries of unequal phase) is described by implementing one simple operation. To evaluate the diffusive flux between two phase boundary nodes (or lattice sites in the Potts model) with different partition coefficients, the concentration of the node with the lower partition coefficient is changed by the amount δK= Kn/K, (6) where K is the partition coefficient of the site and Kn is the partition coefficient of the neighbouring site. This procedure is illustrated in Supplementary DataFig. S14, available for downloading at http://www.petrology.oxfordjournals.org. We may give a difference equation for volume diffusion throughout the domain as: ∇⋅Dk∇C^k= ∑n41Δx2DeCnθi-Cθk, (7) where Dek= 2DnDDn+Dk (8) is the effective chemical diffusivity (harmonic mixing rule) between a central node and its neighbour, the subscript k indicates the component, subscript n indicates one of four (for 2 D models) orthogonal neighbouring nodes, Δx is the grid spacing and θ is a concentration modifier of the central node and θi for the neighbouring node (with respect to the central node), which depends on the partition coefficients. Specifically, If δK<1 then θi=δK−1and θ=1, (9) otherwise θi = θ  =  1 (grain interiors). Finally, we note that in Eq. 7, θi has been used rather than θn, to denote that θi is not a property of the neighbouring node. Instead, it is a property of the central site in relation to the neighbouring node. Also note that a subscript outside of brackets indicates that all quantities in brackets are for the component k. Benchmarks of numerical accuracy in our diffusive partitioning simulations are discussed in Supplementary Data Appendix C, available for downloading at http://www.petrology.oxfordjournals.org. 2D Near-fractional melting Our predictions of aggregate primary melt compositions use the assumption of near-fractional extraction, where melt remains in contact with the residue until a porosity threshold is exceeded, after which all additional melt is mixed with existing pore liquid and then extracted. Because we model diffusive partitioning of trace elements in a 2 D domain at the sub-grain scale and melt may be continuously expelled as it is produced in the case of near-fractional melting, an important problem arises. How do we model the expulsion of melt? It may seem straightforward to remove a melt site by having it migrate toward an edge and then removing it, or by transforming it into a ‘void space’ and allowing the surrounding media to compact. However, such implementations are not simple. The total amount of melt removed is on the order of 25%, or larger for Tp > 1300 °C, resulting in substantial compaction. This compaction will result in the motion of grains with respect to each other, but our Potts model is not capable of advecting whole grains as solid units. Therefore we approximate the effects of melt extraction by assuming that all ‘melt’ is retained in the model domain, but the representative volume of melt decreases with extraction. This means that the Potts model behaves ‘mechanically’ as a batch melting domain, whereas the transport behavior of trace components is consistent with near-fractional melting. So, if a total of 4% melt is produced, but near-fractional melting predicts that residual melt is ϕr =1%, then 4% of the model lattice sites will be melt, but the ‘volume density’ of each melt site is 1/4. This may be introduced into the diffusion equation (Eq. 7) as ∇⋅Dk∇C^k= 1ρv∑n41Δx2DeCnθi-Cθk, (10) where ρv is the volume density (1 for solid, ≤ 1 for liquid). An additional consequence of compaction which we do include is the reduction in the ascent velocity of the residual solid vs as melt is extracted (details in Supplementary Data Appendix B, available for downloading at http://www.petrology.oxfordjournals.org). The reduced ascent velocity is approximately vs, 0exp(-F), where vs, 0 is the initial ascent velocity, although not exactly since there is a finite porosity and a critical threshold for extraction. Erupted compositions If the instantaneously extracted melt volume is Fe= dFdz-dϕrdz, (11) where F is the total melt fraction and erupted compositions are an accumulation of the extracted melt: Ceruptk= ∫zbzLABFeCfkdz/∫zbzLABFedz, (12) where zb is the initial model depth, zLAB is the depth of final extraction (or LAB depth) and [Cf]k is the average melt concentration of component k at the time of extraction. Note that the above formulation is applicable to upward 1 D transport of melt only. It is believed that the production of MORB beneath spreading ridges has a more complicated 3 D character, in which melts generated off-axis migrate horizontally in the direction of younger lithosphere to produce axial crust (e.g. Canales et al., 2012). Such melt focusing should enrich the erupted composition with lower-degree melts from off-axis. However, we ignore this effect because it is unlikely to contribute significantly to ratios of incompatible over compatible elements. For example, we saw in Fig. 4 that La/Sm and Sm/Yb decrease only slightly above a depth of 25–40 km and this will not change in the case of disequilibrium melting. The case will be different for the major rock-forming elements and perhaps absolute concentrations of incompatible elements, but for ratios of incompatible trace elements, off-axis enrichment can be assumed to be unimportant. In our models we will treat the primary aggregate melt as the lava composition, which is then compared to the data. In reality, primary melts are sure to undergo further fractionation due to shallow-level crystallization during ascent and collection in magma chambers (Supplementary DataFig. S49, available for downloading at http://www.petrology.oxfordjournals.org). Because incompatible elements are enriched in the fractionation process, ignoring fractionation will result in an erroneously enriched composition in the source. However, if all OIB (and MORB, where relevant comparisons are made) experience the same amount of crystallization then all aspects of our models are unaffected except for source composition calculations. Thus we will only consider crystallization as a posterior correction in the analysis of our inverse models. Partition coefficients We have considered the experimental datasets of Salters et al. (2002) and Adam & Green (2006, 2010). The coefficients are listed in Table 3. Salters et al. (2002) reported experimental partitioning data for anhydrous spinel and garnet lherzolites (doped synthetic samples) from 1 to 3·4 GPa and 1350–1660°C. Adam & Green (2006, 2010) reported trace element partition coefficients between hydrous (5–10 wt % H2O) basanite melts and mica- and amphibole-bearing garnet lherzolite (doped natural samples). The experimental conditions are low temperature (1000–1200°C) with high volatile contents (5–10 wt % H2O), for pressures of 1–3·5 GPa. The main differences are thus the volatile contents and melting temperature. In addition, Salters et al. (2002) did not report values for La. In lieu of La values, we estimate them based on the regular dependence of partition coefficients on ionic radius in the Adam & Green’s (2010) database. Although these authors attempted to constrain dependences on pressure, temperature and composition, we assume constants for simplicity. We have preferred the values of Salters et al. (2002) in our models because they are more relevant to higher temperature, less hydrated melting conditions. Table 3: Partition coefficients used in this study Olivine/Liq Opx/Liq Cpx/Liq Garnet/Liq Reference La 0·00015 0·0033 0·0328 0·007625 Extrapolation (this work) Sm 0·0013 0·0188 0·1352 0·22125 Salters et al. (2002) Yb 0·0305 0·0988 0·3432 3·95425 Salters et al. (2002) La 0·000001 0·0006 0·03 0·005 Adam and Green (2010) Sm 0·0001 0·011 0·25 0·164 Adam and Green (2010) Yb 0·017 0·077 0·43 4·54 Adam and Green (2010) Olivine/Liq Opx/Liq Cpx/Liq Garnet/Liq Reference La 0·00015 0·0033 0·0328 0·007625 Extrapolation (this work) Sm 0·0013 0·0188 0·1352 0·22125 Salters et al. (2002) Yb 0·0305 0·0988 0·3432 3·95425 Salters et al. (2002) La 0·000001 0·0006 0·03 0·005 Adam and Green (2010) Sm 0·0001 0·011 0·25 0·164 Adam and Green (2010) Yb 0·017 0·077 0·43 4·54 Adam and Green (2010) View Large Table 3: Partition coefficients used in this study Olivine/Liq Opx/Liq Cpx/Liq Garnet/Liq Reference La 0·00015 0·0033 0·0328 0·007625 Extrapolation (this work) Sm 0·0013 0·0188 0·1352 0·22125 Salters et al. (2002) Yb 0·0305 0·0988 0·3432 3·95425 Salters et al. (2002) La 0·000001 0·0006 0·03 0·005 Adam and Green (2010) Sm 0·0001 0·011 0·25 0·164 Adam and Green (2010) Yb 0·017 0·077 0·43 4·54 Adam and Green (2010) Olivine/Liq Opx/Liq Cpx/Liq Garnet/Liq Reference La 0·00015 0·0033 0·0328 0·007625 Extrapolation (this work) Sm 0·0013 0·0188 0·1352 0·22125 Salters et al. (2002) Yb 0·0305 0·0988 0·3432 3·95425 Salters et al. (2002) La 0·000001 0·0006 0·03 0·005 Adam and Green (2010) Sm 0·0001 0·011 0·25 0·164 Adam and Green (2010) Yb 0·017 0·077 0·43 4·54 Adam and Green (2010) View Large It should be realized that although experimental partition coefficients represent the composition of minerals and liquids at equilibrium, they may be used to calculate the evolution of systems far from equilibrium as well, just as we have described above. This is strictly problematic because partition coefficients should depend on their concentration (as chemical potential depends on composition), but this is unlikely to be a significant problem for trace elements. Moreover, it may seem odd that we use partition coefficients taken from experiments which assume equilibrium is achieved over timescales of hours to days and yet claim that melting in real systems may occur far from equilibrium over timescales on the order of 103–105 years. However, the major difference between experimental and real systems is the grain size, which is on the order of 10 s of µm in experiments and 102–105 µm in real systems. Since the time required to achieve equilibrium increases as a function of the grain size squared r2, an increase in grain size corresponds with a squared increase in the time to chemical equilibrium. In our main analyses we will assume that partition coefficients are accurate (i.e. no disequilibrium in experiments) and apply to the melting process, although we will also investigate model results with partition coefficients as free variables. Diffusion coefficients A database of diffusion coefficients for La, Sm and Yb, has been compiled from the literature. Among the data for major upper mantle phases, only REE diffusion in diopside is reported to be significantly different for each element (Van Orman et al., 2001). As these authors do not report an estimate for Sm diffusion in cpx, we estimate Dcpx by interpolation between values for La and Yb using ionic radius. For the other phases, REE diffusion in garnet is based on Van Orman (2002b), diffusion in opx is based on Sano et al. (2011) and olivine on values from Cherniak (2010). We use an activation volume of 10 cm3 mol-1 (Van Orman et al., 2001) for all phases except garnet, which Carlson (2012) estimated to be 20 cm3 mol-1. The dependence of chemical diffusivity on temperature and pressure is given by Dqk= Aqexp⁡Eq+VqpRTk, (13) where k is the diffusing element, q is the phase, A is a pre-exponential coefficient, E is the activation energy, V is the activation volume, p is pressure, R is the gas constant and T is temperature. A tabulation of coefficients for each element and phase is provided in Table 4 and Supplementary Data Fig. S15, available for downloading at http://www.petrology.oxfordjournals.org illustrates the predictions as a function of temperature. Table 4: Chemical diffusivity coefficients used in this study Olivine Opx Cpx Garnet A(La) −9·097 −11·6 −4·22 −8·47 A(Sm) −9·097 −11·6 −4·453 −8·47 A(Yb) −9·097 −11·6 −4·64 −8·47 E(La) 289 233 466 322 E(Sm) 289 233 435 322 E(Yb) 289 233 411 322 V(La) 10 10 10 20 V(La) 10 10 10 20 V(La) 10 10 10 20 Olivine Opx Cpx Garnet A(La) −9·097 −11·6 −4·22 −8·47 A(Sm) −9·097 −11·6 −4·453 −8·47 A(Yb) −9·097 −11·6 −4·64 −8·47 E(La) 289 233 466 322 E(Sm) 289 233 435 322 E(Yb) 289 233 411 322 V(La) 10 10 10 20 V(La) 10 10 10 20 V(La) 10 10 10 20 View Large Table 4: Chemical diffusivity coefficients used in this study Olivine Opx Cpx Garnet A(La) −9·097 −11·6 −4·22 −8·47 A(Sm) −9·097 −11·6 −4·453 −8·47 A(Yb) −9·097 −11·6 −4·64 −8·47 E(La) 289 233 466 322 E(Sm) 289 233 435 322 E(Yb) 289 233 411 322 V(La) 10 10 10 20 V(La) 10 10 10 20 V(La) 10 10 10 20 Olivine Opx Cpx Garnet A(La) −9·097 −11·6 −4·22 −8·47 A(Sm) −9·097 −11·6 −4·453 −8·47 A(Yb) −9·097 −11·6 −4·64 −8·47 E(La) 289 233 466 322 E(Sm) 289 233 435 322 E(Yb) 289 233 411 322 V(La) 10 10 10 20 V(La) 10 10 10 20 V(La) 10 10 10 20 View Large It is known that diffusion coefficients in silicate melts are much higher than those in silicate minerals (Hofmann & Hart, 1978; Koepke & Behrens, 2001). However, high diffusion coefficients impart difficult constraints on numerical stability and computational expense in our models. In the current work, we assume that chemical diffusion in melt is infinite by equilibrating the compositions of all melt sites at each time step. Since diffusion over distance scales of 1 cm occurs over time-scales of years in melt and our simulation times are on the order of 105 years, we believe that this assumption is reasonable. PROBABILISTIC INVERSION To investigate the thermochemical properties in the sources of OIB, we employ a Markov-Chain Monte-Carlo (MCMC) method to fit our disequilibrium models to the global OIB dataset. We fit model results to five observables: La(t), Sm(t), Yb(t), La/Sm(t) and Sm/Yb(t), where t is a bin of seafloor age on which lavas erupt. Models are fit to bins of the data instead of the unfiltered data for computational simplicity and to avoid overrepresentation of age brackets with more samples in our dataset. The binned data, standard deviations and lithospheric thickness for three different values of TLAB are provided in Table 5 (unbinned data are provided in the Supplementary Data, available for downloading at http://www.petrology.oxfordjournals.org). We invert for a total of six unknowns: thermomechanical LAB temperature TLAB, OIB source temperature Tp, oib, mean grain radius r (for a given ascent velocity) and bulk source concentrations for three elements (La, Sm and Yb) Ck. All models assume that the potential temperature of ambient upper mantle is 1300 °C and that the critical porosity is 0·5%. Later we will also briefly discuss inversion results in the case where partition coefficients for all three elements in opx, cpx and garnet are also free parameters. Table 5: Parameters and data used for thermochemical inversion. Compositions are given in 1 Ma1/2 bins based on data from Humphreys and Niu (2009) and our revised seafloor ages t1/2 LT(t)a 1200°C LT(t)a 1225°C LT(t)a 1300°C La/Sm σ1 Sm/Yb σ1 La σ1 Sm σ1 Yb σ1 Ma1/2 km km km - - - - ppm ppm ppm ppm Ppm Ppm 1 16·5 19·4 28·2 1·903 0·432 1·412 0·389 11·341 8·234 4·234 2 2·801 1·177 2 27·9 31·1 40·9 2·502 0·546 1·678 0·249 17·927 9·127 5·529 2·175 3·246 1·369 3 39·2 43·2 55·2 3·188 0·623 2·148 0·26 26·193 11·725 6·829 2·17 3·276 1·171 4 51·1 55·7 69·3 4·158 0·947 2·791 0·374 35·305 16·36 7·724 2·159 2·92 0·893 5 62·8 67·9 83·4 4·61 1·192 3·544 0·487 45·661 22·392 8·64 2·269 2·589 0·839 6 74·5 80·1 97 4·823 1·319 3·957 0·644 49·478 25·071 9·457 2·927 2·638 1·268 7 85·6 91·1 107·5 4·655 1·428 3·987 0·814 45·381 23·597 9·343 3·587 2·677 1·633 8 94·8 99·5 113·8 4·836 1·612 4·09 0·903 45·447 24·613 8·85 3·417 2·363 1·192 9 101·2 105·2 117·2 4·835 1·709 4·27 0·94 47·499 24·07 8·717 2·974 2·133 0·748 10 105·1 108·6 118·9 4·59 1·517 4·234 1·068 42·153 17·819 8·566 2·911 2·138 0·72 11 107·5 110·6 119·8 4·425 1·232 4·092 1·047 40·977 16·207 8·982 3·172 2·375 0·925 12 109 111·9 120·5 5·133 1·622 4·19 0·933 52·246 25·398 9·991 3·327 2·639 1·185 13 110·2 112·9 121 5·821 1·933 4·349 0·941 60·278 32·332 10·483 3·316 2·704 1·247 t1/2 LT(t)a 1200°C LT(t)a 1225°C LT(t)a 1300°C La/Sm σ1 Sm/Yb σ1 La σ1 Sm σ1 Yb σ1 Ma1/2 km km km - - - - ppm ppm ppm ppm Ppm Ppm 1 16·5 19·4 28·2 1·903 0·432 1·412 0·389 11·341 8·234 4·234 2 2·801 1·177 2 27·9 31·1 40·9 2·502 0·546 1·678 0·249 17·927 9·127 5·529 2·175 3·246 1·369 3 39·2 43·2 55·2 3·188 0·623 2·148 0·26 26·193 11·725 6·829 2·17 3·276 1·171 4 51·1 55·7 69·3 4·158 0·947 2·791 0·374 35·305 16·36 7·724 2·159 2·92 0·893 5 62·8 67·9 83·4 4·61 1·192 3·544 0·487 45·661 22·392 8·64 2·269 2·589 0·839 6 74·5 80·1 97 4·823 1·319 3·957 0·644 49·478 25·071 9·457 2·927 2·638 1·268 7 85·6 91·1 107·5 4·655 1·428 3·987 0·814 45·381 23·597 9·343 3·587 2·677 1·633 8 94·8 99·5 113·8 4·836 1·612 4·09 0·903 45·447 24·613 8·85 3·417 2·363 1·192 9 101·2 105·2 117·2 4·835 1·709 4·27 0·94 47·499 24·07 8·717 2·974 2·133 0·748 10 105·1 108·6 118·9 4·59 1·517 4·234 1·068 42·153 17·819 8·566 2·911 2·138 0·72 11 107·5 110·6 119·8 4·425 1·232 4·092 1·047 40·977 16·207 8·982 3·172 2·375 0·925 12 109 111·9 120·5 5·133 1·622 4·19 0·933 52·246 25·398 9·991 3·327 2·639 1·185 13 110·2 112·9 121 5·821 1·933 4·349 0·941 60·278 32·332 10·483 3·316 2·704 1·247 a Lithospheric thickness for different ages (t) and isothermal temperatures. The complete dataset can be found in Supplementary data, available for downloading at http://www.petrology.oxfordjournals.org. View Large Table 5: Parameters and data used for thermochemical inversion. Compositions are given in 1 Ma1/2 bins based on data from Humphreys and Niu (2009) and our revised seafloor ages t1/2 LT(t)a 1200°C LT(t)a 1225°C LT(t)a 1300°C La/Sm σ1 Sm/Yb σ1 La σ1 Sm σ1 Yb σ1 Ma1/2 km km km - - - - ppm ppm ppm ppm Ppm Ppm 1 16·5 19·4 28·2 1·903 0·432 1·412 0·389 11·341 8·234 4·234 2 2·801 1·177 2 27·9 31·1 40·9 2·502 0·546 1·678 0·249 17·927 9·127 5·529 2·175 3·246 1·369 3 39·2 43·2 55·2 3·188 0·623 2·148 0·26 26·193 11·725 6·829 2·17 3·276 1·171 4 51·1 55·7 69·3 4·158 0·947 2·791 0·374 35·305 16·36 7·724 2·159 2·92 0·893 5 62·8 67·9 83·4 4·61 1·192 3·544 0·487 45·661 22·392 8·64 2·269 2·589 0·839 6 74·5 80·1 97 4·823 1·319 3·957 0·644 49·478 25·071 9·457 2·927 2·638 1·268 7 85·6 91·1 107·5 4·655 1·428 3·987 0·814 45·381 23·597 9·343 3·587 2·677 1·633 8 94·8 99·5 113·8 4·836 1·612 4·09 0·903 45·447 24·613 8·85 3·417 2·363 1·192 9 101·2 105·2 117·2 4·835 1·709 4·27 0·94 47·499 24·07 8·717 2·974 2·133 0·748 10 105·1 108·6 118·9 4·59 1·517 4·234 1·068 42·153 17·819 8·566 2·911 2·138 0·72 11 107·5 110·6 119·8 4·425 1·232 4·092 1·047 40·977 16·207 8·982 3·172 2·375 0·925 12 109 111·9 120·5 5·133 1·622 4·19 0·933 52·246 25·398 9·991 3·327 2·639 1·185 13 110·2 112·9 121 5·821 1·933 4·349 0·941 60·278 32·332 10·483 3·316 2·704 1·247 t1/2 LT(t)a 1200°C LT(t)a 1225°C LT(t)a 1300°C La/Sm σ1 Sm/Yb σ1 La σ1 Sm σ1 Yb σ1 Ma1/2 km km km - - - - ppm ppm ppm ppm Ppm Ppm 1 16·5 19·4 28·2 1·903 0·432 1·412 0·389 11·341 8·234 4·234 2 2·801 1·177 2 27·9 31·1 40·9 2·502 0·546 1·678 0·249 17·927 9·127 5·529 2·175 3·246 1·369 3 39·2 43·2 55·2 3·188 0·623 2·148 0·26 26·193 11·725 6·829 2·17 3·276 1·171 4 51·1 55·7 69·3 4·158 0·947 2·791 0·374 35·305 16·36 7·724 2·159 2·92 0·893 5 62·8 67·9 83·4 4·61 1·192 3·544 0·487 45·661 22·392 8·64 2·269 2·589 0·839 6 74·5 80·1 97 4·823 1·319 3·957 0·644 49·478 25·071 9·457 2·927 2·638 1·268 7 85·6 91·1 107·5 4·655 1·428 3·987 0·814 45·381 23·597 9·343 3·587 2·677 1·633 8 94·8 99·5 113·8 4·836 1·612 4·09 0·903 45·447 24·613 8·85 3·417 2·363 1·192 9 101·2 105·2 117·2 4·835 1·709 4·27 0·94 47·499 24·07 8·717 2·974 2·133 0·748 10 105·1 108·6 118·9 4·59 1·517 4·234 1·068 42·153 17·819 8·566 2·911 2·138 0·72 11 107·5 110·6 119·8 4·425 1·232 4·092 1·047 40·977 16·207 8·982 3·172 2·375 0·925 12 109 111·9 120·5 5·133 1·622 4·19 0·933 52·246 25·398 9·991 3·327 2·639 1·185 13 110·2 112·9 121 5·821 1·933 4·349 0·941 60·278 32·332 10·483 3·316 2·704 1·247 a Lithospheric thickness for different ages (t) and isothermal temperatures. The complete dataset can be found in Supplementary data, available for downloading at http://www.petrology.oxfordjournals.org. View Large Among the free parameters, each combination of grain size and source temperature requires an independently computed forward model of disequilibrium melting. Because the bulk source concentration and TLAB do not change the model behavior, they can be computed a posteriori from any case forward model (by an adjustment of concentrations and finding the cumulative melt concentration at the depth of a given isotherm, respectively). Treating the grain size and source temperature as free parameters requires some computational power because the sensitivity to these coefficients is highly nonlinear. We do not invert for diffusion coefficients because they are also coupled to the model behaviour, which would result in a significantly increased computational expense. Instead, diffusion coefficients are assumed known from experimental estimates. We obtain results for disequilibrium melting models with 15 irregularly spaced models by setting the model edge length Lg (Fig. 7) to the following values: 0·05, 0·1, 0·2, 0·3, 0·4, 0·6, 0·8, 1, 1·25, 1·5, 1·75, 2, 2·33, 2·66 and 3 cm. As r ∼Lg/13 in our models, this effectively defines the initial grain size. For each grain size we also vary the source temperature between 1300 °C and 1450 °C at regular intervals of 5 °C. Therefore we have 31 × 15 = 465 melting model results to serve as the basis for the inversion. Using the above computed models, a melting model for a given grain size in the range 0·05–3·0 cm and source temperature in the range 1300–1450°C can be found by interpolation. Since REE concentrations in all models are normalized to the bulk source concentration and since trace element transport is uncoupled (i.e. there are no cross effects), absolute concentrations are obtained by multiplication by bulk source concentrations for each element. Lastly, the REE composition of erupted lavas as a function of seafloor age is calculated by assuming that an isotherm defines the base of the lithosphere and upper mantle temperature uses the thermal plate model of Grose & Afonso (2013) for the case of ambient mantle potential temperature Tp = 1300 °C. We employ a straightforward MCMC algorithm using a random-walk Metropolis method to sample the model predictions with the above 6 free parameters and accumulate posterior distribution functions. Joint misfits to all five observations (La, Sm, Yb, La/Sm and Sm/Yb) control the acceptance rate. The inversion algorithm and modifications employed to resolve a more complicated inversion with free partition coefficients are described in Supplementary Data Appendix D, available for downloading at http://www.petrology.oxfordjournals.org. MAIN RESULTS Since we introduce a novel type of disequilibrium model, we first present results of the general microstructural and microchemical evolution observed in them. Next, we show the general predictions of our models for melt chemistry, which are the only model predictions used in the inversion. We then show results from our inversion study. Petrology and microstructural evolution To illustrate the microstructural and microchemical behavior of our model, we show snapshots of a typical model in Fig. 8 at different stages of evolution. The top row of Fig. 8 shows snapshots of microstructural evolution during decompression of fertile peridotite for a mantle potential temperature Tp = 1300 °C, 1% melting at the dry solidus depth (due to volatiles), ascent velocity vz = 0·65 m yr-1 and model edge length Lg = 0·65 cm. We also provide Supplementary Information (Figs S21 and S25, available for downloading at http://www.petrology.oxfordjournals.org) which depicts microchemistry for models with different kinematic parameters (vz = 0·05–1·0 m yr-1 and Lg = 0·05–1·0 cm). In addition, Supplementary DataFigs S26 and S27, available for downloading at http://www.petrology.oxfordjournals.org show illustrations of the extent of microchemical disequilibrium [C/Ceq]k to complement Figs 8 and 9. Since the evolution is partly driven by changes in modal abundance with decompression, these panels may be compared to the thermodynamic calculations of equilibrium phase abundance in Fig. 5. Also note that because of geometric self-similarity, the model dimensions (or grain size) and ascent velocity may be related by a coefficient for the extent of equilibrium as δe=[De]k/(vzr2), (14) where De is an effective rate of diffusive partitioning for the system. High values of δecorrespond to models closer to equilibrium. Thus, the model in Fig. 8 could also be taken to represent a domain where vz = 0·1 m yr-1 and model edge length Lg = 1·66 cm. Fig. 8. View largeDownload slide Results of a typical melting model with the following properties: vz = 0·65 m yr-1, Lg = 0·65 cm, Tp = 1300 °C and ϕr = 0·5 %. The top panels (row 1) show the microstructure with colours corresponding to phases as labelled in the first panel. The second row panels show the spin of each crystal in the model domain. Crystals existing at the initial condition of the model are represented in blue shades and subsequently grown interstitial crystals have yellow-red shades. The panels in the remaining rows show microchemical predictions with concentrations given relative to bulk source concentration. Fig. 8. View largeDownload slide Results of a typical melting model with the following properties: vz = 0·65 m yr-1, Lg = 0·65 cm, Tp = 1300 °C and ϕr = 0·5 %. The top panels (row 1) show the microstructure with colours corresponding to phases as labelled in the first panel. The second row panels show the spin of each crystal in the model domain. Crystals existing at the initial condition of the model are represented in blue shades and subsequently grown interstitial crystals have yellow-red shades. The panels in the remaining rows show microchemical predictions with concentrations given relative to bulk source concentration. Fig. 9. View largeDownload slide Same type of visualization as in Fig. 8, focussed on predictions of microstructure and microchemistry at 50 km depth for models with different grain sizes and ascent velocities (as given at the top). All properties of each of the five models are the same, except the parameters for model size Lg and ascent velocity vz. Illustrations showing the full range of depths for each of these models are provided in Supplementary Data Figs S14–S21, available for downloading at http://www.petrology.oxfordjournals.org. Fig. 9. View largeDownload slide Same type of visualization as in Fig. 8, focussed on predictions of microstructure and microchemistry at 50 km depth for models with different grain sizes and ascent velocities (as given at the top). All properties of each of the five models are the same, except the parameters for model size Lg and ascent velocity vz. Illustrations showing the full range of depths for each of these models are provided in Supplementary Data Figs S14–S21, available for downloading at http://www.petrology.oxfordjournals.org. The first column of panels (Fig. 8) shows the system at 300 km, which is the initial condition of the melting model. Crystal sizes are roughly proportional to the modes of the phases such that olivine grains (green) are largest, followed by garnet (red), opx (grey) and cpx (blue) with progressively smaller grains. The observation that grains of different phases adopt different relative sizes is important as it shows that high abundance phases will tend to equilibrate slower than others. Also, the fact that grain edges show no tendency to adopt particular orientations indicates that unwanted lattice pinning effects that often occur in Potts models (Janssens et al., 2007) are mitigated by our techniques and grain boundaries show smooth curvatures near triple junctions, as seen in nature. The panels in the second row of Fig. 8 show the spin (or crystal orientation) distribution of the microstructure shown in the first row. Crystals existing at the start of the melting model run are coloured in shades of blue and all subsequently formed grains are coloured in red shades. The remaining rows of panels show the microchemistry discussed later. The second column of panels (Fig. 8) shows that the system after 150 km of decompression has produced abundant interstitial crystals from solid-state phase transformations. The new crystals aggregate primarily at triple junctions and along the edges of decomposing garnets. Although the increasing abundance of orthopyroxene energetically favors growth of previously existing grains (versus the formation of new interstitial grains), phase transformation in our models occurs too rapidly to allow Ostwald ripening of opx. In fact, for the conditions of our models we can resolve three basic stages of microstructural evolution for Ostwald ripening during phase transformation. Initial heterogeneous nucleation from decomposing garnet results in an intergrowth all along garnet boundaries. In nature this intergrowth is known as kelyphite or garnet symplectite. However, in our simulations the timescales of coarsening are long, resulting in coarsening of the reaction rim. If the system evolves slowly enough for the microstructure to proceed further toward equilibrium, the intergrowth will be concentrated around triple junctions. Finally, if the system further proceeds toward equilibrium, small interstitial grains will be absorbed by larger neighbours. The microstructural evolution observed in Fig. 8 is thus far from equilibrium, except at the initial condition. Between Fig. 8 columns 2–3, the system continues to decompress from 150 km to 75 km. Garnet is breaking down more rapidly while small melt fractions are produced. Although the generated melt prefers to gather at triple junctions between solid phases, the melt is continually migrating along grain boundaries. The low energy of solid-liquid boundaries results in a capillary effect in which melt is drawn into solid-solid boundaries. Simultaneously, we can observe the effects of continued decomposition of garnet and growth of orthopyroxene. Because most possible sites of interstitial opx grains are now occupied and decomposition of garnet is accelerating, most grains, large and small, prefer to grow, although some are absorbed by larger orthopyroxenes in close proximity. Because garnet predominantly transforms into orthopyroxene, pre-existing opx grains tend to elongate in the direction of neighbouring garnets. The fourth column of panels (Fig. 8) shows the state of the system at 50 km after the garnet–spinel transition and the production of more melt. Although some of this melt has pooled at triple junctions, much of the melt continues to exist along the edges of two solid phases. All phases which had disequilibrium features (such as concavity) have started to round, although orthopyroxenes remain elongated in the direction of previously existing garnets. It is also interesting to see that some young and small clinopyroxene grains have grown at the expense of large ones. This probably reflects a number of factors, especially the presence of triple junctions where cpx does not exist. These junctions readily accept cpx inherited from neighbouring unstable garnets. The formation of many new triple junctions from the precipitation of opx and olivine further increases this effect. Note that small garnet crystals still persist at a depth of 50 km even though garnet has not been stable over the previous 10 km of decompression (Fig. 5). In fact, while these garnets continue to gradually dissolve, they still remain even at 25 km (Fig. 8, column 5). This behavior is due to two factors. Firstly, the ascent rate in the model shown is high (0·65 myr-1). Models with ascent rates much lower than this are closer to phase equilibrium. The second factor is that our algorithm for phase evolution is not entirely thermodynamically consistent (Supplementary Data Appendix B, available for downloading at http://www.petrology.oxfordjournals.org). In our current model, phase transformations are driven by the difference between the existing and equilibrium phase fractions. In reality, phase transformation is driven by free energy changes associated with nucleation and inter-phase grain boundary migration (Porter & Easterling, 1981). Nevertheless, our tests indicate that the effect on trace element partitioning is small since the amount of residual garnet is much less than 1% (Supplementary DataFigs S11 and S12, available for downloading at http://www.petrology.oxfordjournals.org). In the final snapshot (Fig. 8, column 5), the system has decompressed to a depth of 25 km. Nearly half of the cpx and some opx has dissolved into melt, decreasing their grain sizes. Some opx has been incorporated into olivine, resulting in larger grains. Total melt fraction has increased significantly and melt has formed a thick film on the edges of all grains. The pooling of melt around grains and at triple junctions has resulted in some grain rounding. Equilibrium microstructures produced by the Potts model, with and without melt, have clear similarities to microstructures for melt-bearing silicates and metals observed in the laboratory (e.g. Cmiral et al., 1998; Zhu et al., 2011; Garapic et al., 2013). However, syn-transformational microstructures for silicates (e.g. Cabane et al., 2005; Solferino et al., 2015; Faul, 2000; Kubo et al., 2008; Nishi et al., 2013b) show some apparent differences. The main differences are that exsolution lamellae are not produced and reaction rims formed in solid-state transformations readily coarsen in our models. The reason for both is that our models model coarsening and the temperature-pressure change over thousands of years, whereas in laboratory settings model processes occur over hours to days. Phase transformations occurring over laboratory timescales produce microstructures which reflect nucleation and early coarsening processes, not the steady-state coarsening of a microstructure close to chemical equilibrium. Rapid changes in pressure and temperature also facilitate homogeneous nucleation, which in our models is prevented because homogeneous nucleation is energetically costly. Microchemical disequilibria Snapshots of trace element disequilibria for La, Sm, Yb, La/Sm and Sm/Yb in the microstructure discussed above are also shown in Fig. 8. At the start of the simulation (300 km), we assume that trace elements are distributed at equilibrium according to their partition coefficients. Concentrations are normalized to the initial bulk composition (e.g. [La/Sm]/[La/Sm]bsrc and [La/La]bsrc, where the subscript bsrc denotes bulk source composition). Thus, if the bulk source element concentrations are 1 ppm, then results are in units of ppm and no scaling is required. With the assumption of chemical equilibrium, each phase has a single concentration everywhere. Subsequently, solid-state phase transformations take place and the new phase inherits the local trace element composition of the old phase. The site is no longer in equilibrium, resulting in a flux of mass throughout the system until a new equilibrium is achieved. However, due to the kinetics of diffusion the system remains far from equilibrium. Figure 9 illustrates results for models at 50 km depth with a wide range of grain sizes and ascent velocity. As expected, for small grains and low ascent velocity, the system microchemistry is close to equilibrium. Phase abundances have also been affected by the extent of equilibrium because large grains and higher ascent rates allow less time for phase transformation. In particular, the abundance of garnet at 50 km increases with the extent of disequilibrium, even though garnet is not stable at 50 km at equilibrium. Illustrations of the chemical history over all depths for the models shown in Fig. 9 are provided in the Supplementary Data (Figs S21–S25, available for downloading at http://www.petrology.oxfordjournals.org). Most microchemical disequilibrium features that arise in the models can be linked to the inheritance of chemical signatures, diffusive interaction with a gradually increasing volume of melt, interaction with neighbouring grains in states of disequilibrium, variation in rates of diffusion among phases and with pressure and complex convolutions of these effects over time. In Supplementary Data Appendix E, available for downloading at http://www.petrology.oxfordjournals.org we discuss some of these behaviors for each element. General features of near-fractional disequilibrium primary magmas Figure 10 shows predictions of La, Sm, Yb, La/Sm and Sm/Yb in total accumulated primary magmas as a function of final melting depth for models with Tp, oib = 1300–1450°C, ϕr = 0·5%, r = 0·11–6·71 mm and vz = 0·1 m yr-1. The positions of the wet solidus, dry solidus and garnet-out depth are also shown, as in Fig. 5. A more detailed survey of these basic model results is given in Supplementary DataFig. S29, available for downloading at http://www.petrology.oxfordjournals.org. Fig. 10. View largeDownload slide Disequilibrium melting predictions of La, Sm, Yb, La/Sm and Sm/Yb in cumulative primary melt as a function of final melting depth, source potential temperature and grain size (radius) r. All results are for a critical melt porosity ϕr = 0·5 % and solid ascent velocity vz = 0·1 m yr-1. Results in the left panel are for the smallest grain size and are nearest to equilibrium. Fig. 10. View largeDownload slide Disequilibrium melting predictions of La, Sm, Yb, La/Sm and Sm/Yb in cumulative primary melt as a function of final melting depth, source potential temperature and grain size (radius) r. All results are for a critical melt porosity ϕr = 0·5 % and solid ascent velocity vz = 0·1 m yr-1. Results in the left panel are for the smallest grain size and are nearest to equilibrium. For low extents of melting (below the dry solidus), REE concentrations and their ratios depend strongly on the extent of equilibrium. Models near equilibrium (r = 0·11 mm) show enrichment of progressively incompatible elements (La>Sm>Yb). Models farther from equilibrium (e.g. the case of r = 6·71 mm) are significantly less enriched in incompatible elements. In contrast, compatible elements (Yb) exhibit enriched signatures in low-degree melts for the case far from equilibrium. The sensitivity to disequilibrium in model predictions is related to the competition between the inheritance of compositions from dissolving solids (generally cpx) and the time required for melt to absorb incompatible elements by diffusion. In contrast, concentrations in high-degree primary melts (at shallow final depths of melting) depend little on disequilibrium in the source. In high-degree melts REE ratios quickly approach a constant value equal to the bulk source ratio, while absolute concentrations of incompatible elements are continuously diluted with melting (Supplementary Fig. S29, available for downloading at http://www.petrology.oxfordjournals.org). Final absolute concentrations in primary melts (∼0·3 melt fraction) are about 3·5 times the bulk source value. This relates to the fact that most REE mass exists in garnet and cpx, the mass fraction of garnet+cpx in the source is about 1/3·5 and cpx and garnet are almost completely removed at shallow melting depths in our models (Fig. 5). Inversion results Lithospheric thickness and mantle potential temperature in the OIB source Figure 11 shows Probability Density Function (PDF) contours and a scatter plot of parameters (TLAB, Tp, oib and r) from all model realizations. Overall, the best fitting models occur in a narrow band between where TLAB = 1220 °C and Tp=1300°C and where TLAB = 1300 °C and Tp=1375°C. This band of preferred models is labelled as a red dashed line. Grain size tends to decrease with Tp. We find no realizations in which equilibrium models fit the data (<1·5 mm). The PDF density is generally constant throughout this band, except where TLAB approaches 1300 °C, where density is high. However, we found that the high density of realizations at high TLAB is a misleading artefact of using isotherms in a cooling plate model (which includes thermal effects of melting) to define the base of the lithosphere. Because the latent heat of melting and adiabatic cooling reduces mantle temperature, high temperature isotherms occur deeper in the mantle. For example, the 1300 °C isotherm occurs at about 25 km below the ridge axis in our cooling model. Thus, using 1300 C as the base of the lithosphere is equivalent to the assumption that an OIB source is above 25 km depth below the ridge axis. Nevertheless, because these models result in marginally improved fits, more realizations are found with high TLAB. In our assessment, all models within the band of dense realizations (red dashed line) are equally viable. Figure 11 also shows a secondary band of realizations shifted to TLAB about 50 °C lower. However, this band originates from models which have minimized misfit by achieving excellent fits to old-age data (majority of the data) at the cost of a poor fit to young-age data (minority of the data). Therefore we consider them unreliable. Fig. 11. View largeDownload slide Results of Markov-Chain Monte-Carlo inversion with experimental partition coefficients, showing contours of the 2D PDF (at intervals of 1/6, normalized to the maximum density) overlying a scatter plot of 250 000 realizations. Mantle potential temperature of ambient mantle (which controls the thermal structure of the lithosphere and upper mantle) is 1300 °C for all models. Fig. 11. View largeDownload slide Results of Markov-Chain Monte-Carlo inversion with experimental partition coefficients, showing contours of the 2D PDF (at intervals of 1/6, normalized to the maximum density) overlying a scatter plot of 250 000 realizations. Mantle potential temperature of ambient mantle (which controls the thermal structure of the lithosphere and upper mantle) is 1300 °C for all models. The correlation of Tp, oib with TLAB is explained in Fig. 12. In our models, the demarcation age for thin- and thick-plate OIB (∼12 Ma) occurs above lithosphere where the final melting depth approximately intersects the dry solidus. More specifically, inverse models predict that the dry solidus intersects the LAB below lithosphere of about 16 Ma (black vertical bar), or slightly older than our demarcation age (blue vertical bar). This is expected because although Light/Heavy REE values are rapidly diluted as the rate of melt production with decompression greatly increases around the dry solidus depth, dilution is not instantaneous. For Tp, oib = 1300 °C the dry solidus depth below the demarcation age intersects a mantle temperature of about 1225 °C. Hence, models with TLAB = 1225 °C provide the best fits (A’ in Fig. 12). For higher Tp, oib the dry solidus depth increases and thus TLAB must also increase (B’ in Fig. 12). Eventually the solidus is so deep that the required isotherm no longer tracks the base of the lithosphere (C’ in Fig. 12). Fig. 12. View largeDownload slide Illustration showing the relationship between the choice of LAB isotherm, melt source potential temperature (which constrains the dry solidus depth), and explanation of the chemical data. (a) Sm/Yb in OIB as a function of lithospheric age as in Fig. 2 with 1 My1/2 bins as black circles. Also shown is the approximate age at which the dry solidus intersects the LAB and the estimated demarcation age for thin- and thick-plate OIB. (b) Upper mantle thermal structure and LAB depths (for the possible temperatures 1200°C, 1300°C, and 1350°C) as a function of lithospheric age. The large white circles show the depths at which the dry solidus occurs for three possible OIB source temperatures (A'=1300°C, B'=1400°C, and C'=1475°C). Dry solidus depths use our melting parameterization. Isotherms use the thermal plate model of Grose and Afonso (2013). See text for discussion. Fig. 12. View largeDownload slide Illustration showing the relationship between the choice of LAB isotherm, melt source potential temperature (which constrains the dry solidus depth), and explanation of the chemical data. (a) Sm/Yb in OIB as a function of lithospheric age as in Fig. 2 with 1 My1/2 bins as black circles. Also shown is the approximate age at which the dry solidus intersects the LAB and the estimated demarcation age for thin- and thick-plate OIB. (b) Upper mantle thermal structure and LAB depths (for the possible temperatures 1200°C, 1300°C, and 1350°C) as a function of lithospheric age. The large white circles show the depths at which the dry solidus occurs for three possible OIB source temperatures (A'=1300°C, B'=1400°C, and C'=1475°C). Dry solidus depths use our melting parameterization. Isotherms use the thermal plate model of Grose and Afonso (2013). See text for discussion. Predicted compositions of OIB lavas The compositions of OIB predicted by our inverse models are shown for La, Sm, Yb, La/Sm and Sm/Yb in Fig. 13a–e. The white–blue–red colour scale is shown in Fig. 13e and indicates the probability density of all model realizations. These model results are compared to both the OIB data (red and blue markers with error bars), 1 Ma1/2 bins of the OIB data (black lines with filled circles) and the MORB distribution (pink histogram at the left of each panel). Although there is general agreement with the data, the fits are imperfect. Compared to the binned data, La is underpredicted over most seafloor. Although a better fit to the La data might have been provided by adjusting the extent of equilibrium (via the grain size), these models were rejected in the inversion because such an adjustment results in poorer fits to more compatible elements. The mismatch may thus be attributed to another factor, such as the choice of partition coefficients. Sm, Yb and their ratio also appear to modestly underpredict the data over the thinnest lithosphere (<5 Ma), but their fits are overall excellent. The transition between thick-plate and thin-plate signatures appears considerably more gradual in all of the data compared to our models. This is perhaps most clear in the case of La/Sm, which features a nearly instantaneous transition in our models around 15 Ma. The sharp transitions in our models are a result of idealistic assumptions regarding the chemical and thermomechanical properties in the melt source, ambient upper mantle and oceanic lithosphere. Regional variations in such properties and interactions are likely to affect the final depth of melting. We have not accounted for such an uncertainty in our models. More concerning is that our models markedly overestimate La/Sm in OIB erupted on the thinnest lithosphere (<5 Ma), discussed below. Fig. 13. View largeDownload slide Comparison of empirical data (blue diamonds and red squares with error bars), bins of these data (connected black circles) and the predictions of our inverse models. From the top to bottom rows panels show data and predictions for La, Sm, Yb, La/Sm and Sm/Yb in lavas. The left column panels show the actual a posteriori model distributions, the middle column panels show the model distributions convolved with the empirical distributions observed in MORB and the right column panels show the models convolved with the empirical distributions for thick-plate OIB (t > 12 Ma). The colour scale for the probability density (normalized to the maximum density) is shown in panel (e). The MORB distributions are shown at the left of each panel. Fig. 13. View largeDownload slide Comparison of empirical data (blue diamonds and red squares with error bars), bins of these data (connected black circles) and the predictions of our inverse models. From the top to bottom rows panels show data and predictions for La, Sm, Yb, La/Sm and Sm/Yb in lavas. The left column panels show the actual a posteriori model distributions, the middle column panels show the model distributions convolved with the empirical distributions observed in MORB and the right column panels show the models convolved with the empirical distributions for thick-plate OIB (t > 12 Ma). The colour scale for the probability density (normalized to the maximum density) is shown in panel (e). The MORB distributions are shown at the left of each panel. Convolutions with empirical distributions Since the binned data represent averages, the models from our inversion (Fig. 13a–e) represent predictions of mean OIB concentration. Scatter in the data within the thin- and thick-plate OIB datasets may be attributed to melting and crystallization processes and compositional heterogeneity in the source, effects not accounted for in our models. As discussed earlier, compositional distributions of MORB and thick-plate OIB are remarkably similar, differing significantly only in the most incompatible elements (e.g. La). Because processes unaccounted for in our models may be captured by the distributions in these datasets, illustrating model predictions with the distributions mapped onto them may provide insight. We accomplish this with a type of logarithmic convolution (*l) of the empirical distribution function gk(x) with the inverse model distribution fk(x) for each element k (Fig. 13a–e) using the equation gk,n*lfkx= ∫0∞fkτgkxnkτdτ, (15) where the subscript n indicates that g(x) is mean-normalized, x and τ are concentration and nk is the mean concentration of the function g(x): nk= ∫0∞xgk(x)dx/∫0∞gk(x)dx. (16) A straightforward example of convolution is provided in Supplementary DataFig. S48, available for downloading at http://www.petrology.oxfordjournals.org. The result of convolving the inverse model distributions (Fig. 13a–e) with the MORB (Fig. 3c) and thick-plate OIB distributions (Fig. 3d) are shown in Fig. 13f–j and k–o, respectively. There are at least two main insights provided by this exercise. Firstly, the overlap of the distributions with the OIB data after convolution with either dataset is substantial, showing that employing our techniques can potentially explain the bulk of the OIB data. Secondly, because the MORB distributions are nearly the same as the thick-plate OIB distributions, the MORB convolution has nearly the same explanatory power. Convolutions also show how the means do not always coincide with the modes. For example, because the La and La/Sm distributions are skewed to low concentrations, sample density is expected to be higher below the mean. For the La dataset in particular this reveals a problem, as our predictions seem to noticeably underpredict the concentrations in OIB, particularly for the convolution with the MORB distribution. This is probably due to a combination of the previously remarked underprediction of our models for mean La in thick-plate lavas, as well as the greater skewing of La in MORB. Predicted MORIB and observed MORB values Model compositions in PDFs at zero age (Fig. 13) are the predicted compositions for OIB erupted on ridge axes, or MORIB (Mid-Ocean Ridge Island Basalt). These functions are also shown in Fig. 14a compared to mean MORB from Gale et al. (2013) and enriched MORB (E-MORB) from Sun & McDonough (1989). The usefulness of the MORIB concept is that although OIB sources may be thermally or chemically different from MORB sources, the screening effect of lithospheric thickness and the associated variations in extent of melting complicate direct comparison of chemical (or isotopic) datasets between these settings. The removal of the screen in nature is only accomplished for OIB erupted on ridge axes. In this way the concept of a MORIB geochemical signature has some analogies to the geophysical concept of a mantle potential temperature. Overall, our model MORIB values show mixed agreements with observed MORB. Observed mean La in MORB is low, but still within our posterior distribution. However, observed mean Sm and Yb in MORB are significantly higher than our predictions. Posterior distributions for Sm/Yb are in good agreement with the data, while predicted La/Sm is significantly higher than the observed mean. We tested whether or not the misfit could be due to disequilibrium by assuming that OIB ascent is fast (0·5 m yr-1) compared to a hypothetical ascent rate of the MORB source (0·05 m yr-1). This is done by recalculating all models which fit the OIB data to the reduced ascent rate (Eq. 14). Although such an adjustment markedly affects low-degree melt (OIB) compositions (Supplementary DataFig. S33, available for downloading at http://www.petrology.oxfordjournals.org), high-degree melts (MORB, MORIB) are negligibly affected (Fig. 14a;Supplementary DataFig. S33, available for downloading at http://www.petrology.oxfordjournals.org). There are several other possible contributors to the misfit, such as inaccuracies in our model of melting processes or the choice of partition coefficients. On the other hand, the discrepancy between observed MORB and modelled MORIB values may be real, which could imply that the OIB source is more enriched than the source of MORB. The markedly better fit of our MORIB calculations to E-MORB values from Sun & McDonough (1989) agrees with this view. E-MORB La fits our MORIB prediction exactly, E-MORB Sm is within our MORIB distribution and E-MORB Yb is only marginally above our models. Fig. 14. View largeDownload slide Statistics for models from our Markov-Chain Monte-Carlo inversion. (a) Predictions of the concentrations in quickly ascending OIB sources erupting on ridge axes (MORIB, solid lines) and for a slowly ascending MORB (dashed lines). Because MORIB may be equivalent to MORB if source compositions are the same, predictions are compared to empirical MORB and E-MORB concentrations from G13 (Gale et al., 2013) and SM89 (Sun & McDonough, 1989), respectively. (b) Predictions of bulk source concentration for La, Sm, Yb, La/Sm and Sm/Yb, uncorrected for crystallization effects. Distributions are compared to predictions of previous studies for DM from reference WH05 (Workman & Hart, 2005), PM from SM89 (Sun & McDonough, 1989) and PM from LK07 (Lyubetskaya & Korenaga, 2008). See text for discussion. Fig. 14. View largeDownload slide Statistics for models from our Markov-Chain Monte-Carlo inversion. (a) Predictions of the concentrations in quickly ascending OIB sources erupting on ridge axes (MORIB, solid lines) and for a slowly ascending MORB (dashed lines). Because MORIB may be equivalent to MORB if source compositions are the same, predictions are compared to empirical MORB and E-MORB concentrations from G13 (Gale et al., 2013) and SM89 (Sun & McDonough, 1989), respectively. (b) Predictions of bulk source concentration for La, Sm, Yb, La/Sm and Sm/Yb, uncorrected for crystallization effects. Distributions are compared to predictions of previous studies for DM from reference WH05 (Workman & Hart, 2005), PM from SM89 (Sun & McDonough, 1989) and PM from LK07 (Lyubetskaya & Korenaga, 2008). See text for discussion. It is worth mentioning that for all geochemical signatures, poor fits from our models to MORB seems to coincide with poor fits to the OIB data. In other words, the OIB data seem to trend toward the MORB distribution, suggesting that a better model fit to the OIB data would result in models with better fits to MORB. For example, the first (1 Ma1/2) bin for Yb coincides perfectly with the MORB mean of about 3 ppm (Fig. 13c), but our models are much lower (∼2 ppm) at zero-age than MORB. Similar misfits are seen in Sm and La/Sm. This indicates that if our melting models better fit the OIB data, they would also better fit MORB. It will be worth revisiting this problem in the future with improved melting models and data analysis techniques. Bulk source composition and fractional crystallization PDFs for bulk source concentration predicted by our inversion are shown in Fig. 14b and are compared to compositions for Primitive Mantle (PM) from Sun & McDonough (1989) and Lyubetskaya & Korenaga (2008) and Depleted Mantle (DM) from Workman & Hart (2005). Predictions for Yb are in good agreement or are slightly higher than DM and PM estimates. Sm is significantly higher than DM estimates and slightly higher than PM. Predicted La is far higher than DM and still significantly higher than the PM value. These progressively higher concentrations may originate from model assumptions (e.g. source homogeneity), elevated source concentrations, or crystallization effects. The effects of crystallization on our model calculations of bulk source concentration can be tested using an inversion of the Rayleigh equation for fractional crystallization, CC0k=1-XKk,bulk-1, (17) where [C]k is the uncorrected bulk source concentration of element k as predicted by our inverse models, [C0]k is the corrected prediction of bulk source concentration, X is the mass fraction of crystallization and Kk, bulk is the solid/melt bulk partition coefficient. Notice that incompatible element concentration decreases with increasing crystallization as we are computing the effect on source composition, not lava composition. We assume bulk mineral/melt partition coefficients of KLa = 0·048, KSm = 0·179 and KYb = 0·242, which were reasoned to apply to crystallization of a 0·42/0·41/0·18 cpx/plag/ol solid matrix by O’Neill & Jenner (2012). Calculations are shown in Fig. 15. The three clouds of markers show 250 000 inverse model realizations of bulk source composition for 0%, 60% and 90% fractional crystallization. Markers are coloured according to their mean grain size. Because crystallization increases the concentration of incompatible elements in erupted lavas, increasing extents of crystallization require that the bulk source concentration decreases. Thus, corrected models predict lower values which are generally in better agreement with the candidate mantle compositions. In detail, however, calculations show disagreements. No high probability (>1/8 maximum normalized probability density, Fig. 15) models agree with either DM, PM, or C1 chondrite values and no models agree with DM or PM compositions. Only C1 chondrite can fit model compositions with no adjustments to bulk partition coefficients and those models require about 75% fractional crystallization. PM compositions from Sun & McDonough (1989) are barely on the fringes of our models if crystallization is low (<20%) and grain size is small (∼1 mm). PM compositions from Lyubetskaya & Korenaga (2008) are also on the fringes of our model predictions, but require larger crystallization amounts (35–50%). Fig. 15. View largeDownload slide Plots showing how inverse models of bulk source concentration are corrected for fractional crystallization, using bulk partition coefficients from O’Neill & Jenner (2012). The cloud of markers around the open circle labelled 0% crystallization represents uncorrected model realizations from the inversion (compare to PDFs of these models in Fig. 14b). Each marker is coloured according to the mean grain size (radius) for the model, contours (normalized to the maximum density and at intervals of 0.125) indicate the density of the uncorrected models and the maximum density is indicated by the small open circle. All models corrected for 60% and 90% crystallization are also shown, although contour plots are not shown to reduce clutter. Additional predictions for the maximum probability model are shown at 5% intervals (small filled black circles). Candidate ‘mantle compositions’ are also shown for Depleted Mantle (DM) from reference WH05 (Workman & Hart, 2005), values for C1 chondrite and Primitive Mantle (PM) from SM89 (Sun & McDonough, 1989) and PM from LK07 (Lyubetskaya & Korenaga, 2008). See text for discussion. Fig. 15. View largeDownload slide Plots showing how inverse models of bulk source concentration are corrected for fractional crystallization, using bulk partition coefficients from O’Neill & Jenner (2012). The cloud of markers around the open circle labelled 0% crystallization represents uncorrected model realizations from the inversion (compare to PDFs of these models in Fig. 14b). Each marker is coloured according to the mean grain size (radius) for the model, contours (normalized to the maximum density and at intervals of 0.125) indicate the density of the uncorrected models and the maximum density is indicated by the small open circle. All models corrected for 60% and 90% crystallization are also shown, although contour plots are not shown to reduce clutter. Additional predictions for the maximum probability model are shown at 5% intervals (small filled black circles). Candidate ‘mantle compositions’ are also shown for Depleted Mantle (DM) from reference WH05 (Workman & Hart, 2005), values for C1 chondrite and Primitive Mantle (PM) from SM89 (Sun & McDonough, 1989) and PM from LK07 (Lyubetskaya & Korenaga, 2008). See text for discussion. An alternative approach to evaluating our results is to calculate the partition coefficients required for models to fit the candidate mantle compositions which, from Eq. 17, can be calculated as Kk,bulk= log⁡CC0klog⁡(1-X)+1. (18) The results of this experiment are shown in Fig. 16. Bulk partition coefficients which fit the source composition [C]k are calculated as a function of crystallization amount for four candidate source compositions and each element. We show calculations for the maximum probability model (coloured lines) and calculations for all model realizations are shown as probability density functions for each crystallization amount. A self-consistent set of partition coefficients is found on any given vertical line. For example, for a crystallization amount of 80%, the maximum probability bulk partition coefficients for an assumed source concentration of PM from Lyubetskaya & Korenaga (2008) are KLa = 0·25, KSm = 0·6 and KYb = 0·8. However, low probability realizations allow for a significant amount of uncertainty about these values. The main result of this experiment is to show that La requires large crystallization amounts for any candidate source composition and that partition coefficients for Sm and Yb must be significantly higher than assumed by O’Neill & Jenner (2012) to fit our models. Fig. 16. View largeDownload slide Plots showing the bulk partition coefficient required to satisfy four candidate source compositions (as indicated for each panel) as a function of crystallization amount. Coloured lines show functions for the maximum probability model (Figs 14 and 15) and the greyscale gradients shows the probability densities from all 250 000 model realizations. Horizontal lines for La, Sm and Yb are the bulk partition coefficients from O’Neill & Jenner (2012). Note that the density functions for each element are stacked to save space. Fig. 16. View largeDownload slide Plots showing the bulk partition coefficient required to satisfy four candidate source compositions (as indicated for each panel) as a function of crystallization amount. Coloured lines show functions for the maximum probability model (Figs 14 and 15) and the greyscale gradients shows the probability densities from all 250 000 model realizations. Horizontal lines for La, Sm and Yb are the bulk partition coefficients from O’Neill & Jenner (2012). Note that the density functions for each element are stacked to save space. Inversion with variable partition coefficients In our previous inversion we decided to treat partition coefficients as experimentally known constants. However, the values of partition coefficients are critical for chemical evolution whether the system is far from equilibrium or not. Moreover, there are at least three reasons that a direct application of experimental partition coefficients may result in a model which does not represent decompression melting in the Earth. Firstly, the uncertainties in the calculation of partition coefficients from experimental measurements can be large, even several factors away from the means (Adam et al., 2014). Second, although we have used constants, partition coefficients vary with pressure, temperature and compositional changes in the melt and solid phases. In the future it may be useful to employ lattice strain models which account for such effects (e.g. Sun & Liang, 2014). Thirdly, experimental determination of equilibrium partition coefficients requires that the system is at equilibrium, which may not be correct for slowly diffusing trace elements in a system where crystals may grow rapidly (e.g. Schneider & Eggler, 1986; Adam et al., 1993; 1997). As part of our investigation we also performed an inversion treating partition coefficients as free parameters (Supplementary Data Appendix D, available for downloading at http://www.petrology.oxfordjournals.org). All aspects of the inversion are the same as before, except that OIB source temperature is fixed at 1300 °C and partition coefficients for La, Sm and Yb each in opx, cpx and garnet are free parameters (between 0·05 and 5 times the experimental values). Although this inversion is far more expensive, the lack of cross-diffusion effects among REEs allows some simplifications for computational tractability. As all partition coefficients in olivine are low, we only consider coefficients in opx, cpx and garnet, thus adding three free variables to the inversion. Removing Tp, oib as a free parameter this is reduced to only two new parameters. The main results for partition coefficients are shown in Fig. 17 and a fuller range of results are explored in Supplementary Data Figs S43–S47, available for downloading at http://www.petrology.oxfordjournals.org. In summary, we obtained models which provide much improved fits to some of the available data. For absolute concentrations and ratios in OIB and MORB, misfits are reduced, although there remains some underprediction of Yb and overprediction of La in MORIB. The largest misfit is in the MORIB prediction of La/Sm which, although improved, is still about 40% higher than the observed value. The posterior distributions for partition coefficients (Fig. 17) show significant disagreements with the experimental values by a factor of up to 5 (although this was the maximum deviation from experimental values considered). The most robust predictions are that partition coefficients KLacpx/melt and KSmcpx/melt are significantly higher and KSmgrt/melt and KYbgrt/melt are significantly lower than experimental values (Fig. 17). The resulting models show relatively good fits to REE concentrations in abyssal peridotite cpx grains (Supplementary Data Appendix L, Fig. S47, available for downloading at http://www.petrology.oxfordjournals.org). Also, the extent of equilibrium in the models is markedly greater than predicted by our inversion with experimental coefficients (Supplementary Fig. S45, available for downloading at http://www.petrology.oxfordjournals.org). In fact, the inversion with free partition coefficients includes near-equilibrium models. The thick lines in Fig. 17 represent the samples closest to equilibrium in the distribution, showing that KYb coefficients are more compatible and KLa coefficients are more incompatible than the unconstrained distribution. Consequently, the inability of our original equilibrium models (Fig. 6) to fit the data may simply be due to the use of partition coefficients which do not describe the effective values of the melting system. Fig. 17. View largeDownload slide Results of inversion with variable partition coefficients. Each panel shows a posteriori probability density functions (PDF) for opx/melt (a–c), cpx/melt (d–f) and grt/melt (g–i) partition coefficients for La (a, d, g), Sm (b, e, h) and Yb (c, f, i). The vertical scale is the probability density and the horizontal scale is the partition coefficient (top of each panel) and partition coefficient adjustment from default values (bottom of each panel). In each panel the thin line is the probability density of all inverse model realizations, whereas the thick line labelled ‘eq’ represents the near-equilibrium samples. The default values in our model (partition coefficient adjustment = 1) are the values from Salters et al. (2002) and our extrapolation for La in Table 3. Partition coefficients from Adam & Green (2006) and Suhr et al. (1998) are also shown as coloured diamonds and circles, respectively. Note that 5x the default value is the maximum value allowed by our inversion. Fig. 17. View largeDownload slide Results of inversion with variable partition coefficients. Each panel shows a posteriori probability density functions (PDF) for opx/melt (a–c), cpx/melt (d–f) and grt/melt (g–i) partition coefficients for La (a, d, g), Sm (b, e, h) and Yb (c, f, i). The vertical scale is the probability density and the horizontal scale is the partition coefficient (top of each panel) and partition coefficient adjustment from default values (bottom of each panel). In each panel the thin line is the probability density of all inverse model realizations, whereas the thick line labelled ‘eq’ represents the near-equilibrium samples. The default values in our model (partition coefficient adjustment = 1) are the values from Salters et al. (2002) and our extrapolation for La in Table 3. Partition coefficients from Adam & Green (2006) and Suhr et al. (1998) are also shown as coloured diamonds and circles, respectively. Note that 5x the default value is the maximum value allowed by our inversion. Although fits to the data have improved by treating partition coefficients as free parameters, it is not clear how meaningful the results are because now there are at least two types of coefficients which are unknown (partition coefficients and the extent of equilibrium) and can affect the composition of low-degree melts in similar ways. Basically, more compatible bulk partition coefficients for La and Sm and more incompatible coefficients for Yb can produce the effect of disequilibrium on melt compositions. Nevertheless, the results demonstrate that it is not possible to reasonably constrain the extent of equilibrium in the source of OIB without constraints on partition coefficients. Moreover, our models show that if melting actually occurs near equilibrium, partition coefficients may be significantly different from typically reported experimental values. In turn, because our models demonstrate that melting is a disequilibrium process if grain sizes are large (r > 0·1 mm), we may provisionally conclude that if equilibrium melting applies, then mantle grain size is small, the diffusion coefficients in our database are erroneously low, or important aspects of melting processes have not been captured by our models. DISCUSSION Is disequilibrium a controlling factor in compositional distributions in lavas? Our models show that the effect of disequilibrium on melt compositions is a significant adjustment in the concentration of incompatible elements in thick-plate OIB. Equilibrium models predict high concentrations of incompatible elements, whereas disequilibrium models do not allow the time necessary for this partitioning to occur, thus lowering incompatible element concentrations in the low-degree melts which form thick-plate OIB. From these facts it may be hypothesized that variations in OIB signatures might be explained by variations of grain size in the source. However, disequilibrium imparts no control on Light/Heavy REE in high-degree melts (Fig 10 and Supplementary DataFig. S29, available for downloading at http://www.petrology.oxfordjournals.org). Thus, if lavas represent well-mixed aggregate melts, then variations due to disequilibrium will be present in thick-plate OIB, minor or negligible in thin-plate OIB and completely absent in MORB. Our observation that compositional distributions are nearly identical in MORB and thick-plate OIB (Fig. 3), would seem to confirm that disequilibrium is not responsible for variance in either dataset. This would also imply that variations in mantle grain size are not significant over length-scales larger than the melting region (∼10–100 km). On the other hand, if compositional distributions originate from incomplete mixing of primary magmas, part of the distribution must originate from the sampling of low-degree melts, whose signatures are sensitive to disequilibrium in the source. It has become widely accepted that much of the compositional diversity in basalts originates from incomplete mixing (e.g. Maclennan, 2008; O’Neill & Jenner, 2012; Rudge et al., 2013; Shorttle, 2015), although the degree to which low-degree melts are sampled in basalts is debated. It will be worthwhile to extend such studies in the future to include the effects of disequilibrium. Are major elements at thermodynamic equilibrium during melting? Our models roughly assume that major element partitioning and phase transformations occur rapidly, so that phase abundances are close to equilibrium. The effective assumption is thus that because trace elements do not thermodynamically participate in interface migration and volumetric free energy changes during phase transformation, trace elements can maintain a disequilibrium state according to their diffusion kinetics, while major elements are near equilibrium. However, it is possible that major element partitioning is a disequilibrium process as well, in which case our models may significantly underestimate disequilibrium effects. Equilibrium compositions of alkali basalts and low pressure pyroxenes require that Al2O3 partitions into basaltic liquid during melting and into pyroxenes during garnet dissolution (Supplementary DataFig. S13, available for downloading at http://www.petrology.oxfordjournals.org). Thus, major element disequilibrium should be mainly controlled by Si-Al inter-diffusion in these phases. Relevant data are limited, but Nishi et al. (2013a) and Van Mierlo et al. (2013) measured Si–Al inter-diffusion coefficients in the reaction of majorite and pyropic garnet and found rates substantially lower than for REE. Although their measurements were obtained at high pressure (12–20 GPa), a low activation volume measured by Van Mierlo et al. (2013), applied to both of their estimates, results in diffusion rates at zero pressure consistent with REE values (Supplementary DataFig. S15, available for downloading at http://www.petrology.oxfordjournals.org). These estimates suggest that major element diffusion and phase transformation during peridotite decompression may be as sluggish as REE diffusion. Because a physical model of grain-scale, non-equilibrium thermodynamics is beyond the scope of the present work, we assumed that major elements are not far from equilibrium. Although this is a reasonable first step in this research direction, it will be important to elucidate the details of microstructural evolution with a self-consistent, non-equilibrium thermodynamics modelling framework. Temperature in the source of OIB A significant prediction of our inversion is that the temperature in the source of OIB is unlikely to be >1400°C (Fig. 11). However, it is important to recognize that this is a statistical inference based on fitting melting model predictions to the global OIB dataset. Overall, the prediction arises due to the position of the chemical step in lava compositions around 12 Ma and the depth of the dry solidus in our melting parameterization (Fig. 12). We cannot necessarily infer the source temperature for individual islands because bulk source compositions and crystallization amounts are unknown. Our approach of fitting a single melting model to the global dataset assumes that the source temperature, bulk source composition and its variance are similar for all island groups. The absence of any significant overlap of MORB and thin-plate OIB with thick-plate OIB (Table 1) suggests that this assumption is correct. In contrast, Hawaii and Reunion are notable exceptions showing La/Sm and Sm/Yb distributed toward lower values than typical thick-plate OIB (Supplementary DataFigs S7 and S8, available for downloading at http://www.petrology.oxfordjournals.org), consistent with higher degree melting. A possible scenario is that the majority of islands originate from low (<1400°C) potential temperature sources, whereas the classical hotspots (e.g. Iceland, Galapagos, Hawaii, Reunion) originate from thermal, chemical, or thermochemical plumes. In addition, we note that the thin-plate OIB data are dominated by the Galapagos Island group. Because the consistently depleted signatures in these lavas contribute to the presence of the step-like feature in the data, this is a nontrivial observation. If the depleted signatures in the Galapagos Island lavas originate from a high temperature and compositionally (in terms of trace elements) MORB-like source, the position of the chemical step may be a less meaningful feature, or even an artefact of limited data. In the future it will be important to expand the OIB dataset, particularly for islands on young seafloor. The best case scenario is a single hotspot track which traverses a wide span of seafloor ages, including the youngest (<12 Ma) and much older seafloor. A discontinuity in signatures as a function of plate thickness might then constrain the temperature of the source. Constraints on the thermomechanical lithosphere–asthenosphere boundary Deformation observed to occur in geodynamic models with realistic viscosities is consistent with TLAB of about 1200 °C (e.g. Afonso et al., 2008a; Ballmer et al., 2011; Afonso et al., 2016), which agrees well with our predictions (Fig. 11). The band of dense realizations which we consider to be most reliable is approximately fit by the polynomial (Fig. 11, red dashed line) TLAB(Tp=1300°C)=0·000159075Tp,oib3−0·654916Tp,oib2+899Tp,oib−410159, (19) with model realizations distributed over a range of about ±10°C, skewed to high values. However, these predictions are only for an ambient mantle potential temperature Tp of 1300 °C. For a different Tp, the adjustment to our predictions might be reasonably approximated by a proportional shift to higher TLAB as TLAB(Tp,a)= TLAB(Tp=1300°C)×Tp,a/1300, (20) where Tp,a is the actual mantle potential temperature of ambient mantle and Tp is the value used in our inversions (1300 °C). Equation 20 should apply to any reasonable Tp (perhaps 1250–1450°C), including cases where the OIB source is colder than ambient mantle (Tp, oib<Tp). For example, for Tp=1400°C and Tp,oib = 1300 C, TLAB = 1315 °C. It is likely that our calculations are sensitive to the detailed nature of the mechanical interaction between the melting source material (whether it is a thermal plume or otherwise) and the base of the lithosphere, as well as the details of melt production and segregation. For example, a thermal plume may heat and mechanically erode the base of the lithosphere, decreasing lithospheric thickness and allowing it to rise to a shallower depth. Also, not all parts of the source will melt to the same degree and aggregate melts may sample a large range of depths. These details of mechanical interaction and melt migration are ultimately important for linking the parameter TLAB predicted by our geochemical inversion to the thermo-mechanical LAB in geodynamic models. Radiogenic isotope systematics and lithological heterogeneity of the mantle We have shown that REE systematics of OIB and MORB may be explained by melting of lithologically homogeneous fertile peridotite. However, based mainly on isotopic data it has become widely recognized that the mantle is at least compositionally and may also be lithologically, heterogeneous (Stracke et al., 2003; Willbold & Stracke, 2006). Predicting isotopic signatures of lavas consistent with heterogeneous sources, disequilibrium melting and the lithospheric thickness effect is reserved for later study. However, as we found for REE, isotopic signatures in OIB exhibit striking correlations with lithospheric thickness which must be addressed by future studies with capable models. Using the dataset of Niu et al. (2011) with our seafloor ages and thermal plate model we plot 143Nd/144Nd, 176Hf/177Hf and 87Sr/86Sr as a function of lithospheric thickness in Fig. 18. MORB and thick- and thin-plate OIB distributions are also shown. As in the REE dataset, lavas tend to adopt MORB to moderately enriched isotopic values for eruptions on thin lithosphere, while lavas on thick plates exhibit more enriched signatures (high 87Sr/86Sr and low 143Nd/144Nd and 176Hf/177Hf). As previously recognized (Ito & Mahoney, 2005a, 2005b), the range of isotopic scatter is larger in islands formed on thick lithosphere compared to islands on thin lithosphere. Variance of the distributions in the MORB and thin- and thick-plate OIB subsets can be explained as a lithospheric thickness effect on melts generated from sources which are compositionally similar and contain an enriched component. The enriched component is minor in terms of mass and is more strongly sampled in thick-plate OIB. Thin-plate OIB signatures are diluted by a more predominant depleted component and MORB is a maximally diluted end-member. Plots of isotopic ratios versus La/Sm and Sm/Yb are also shown in Fig. 18. Fig. 18. View largeDownload slide Plots of isotopic ratios for all islands in the dataset of Niu (2011). Error bars are 1σ standard deviations of the intra-island dataset. The pink density function projecting to the left of each panel represents values observed in MORB from Gale et al. (2013). As in previous figures, red points represent lavas erupted on thin lithosphere and blue points represent lavas erupted on thick lithosphere. To the right of plots (a–c) are shown distribution functions for lavas on thin and thick lithosphere. Panels (d–f) show La/Sm vs isotopic ratios and panels (g-i) show Sm/Yb vs isotopic ratios. MORB values are shown in panels (d–i) as black dots. Lithospheric thickness is calculated assuming TLAB = 1225 °C. Fig. 18. View largeDownload slide Plots of isotopic ratios for all islands in the dataset of Niu (2011). Error bars are 1σ standard deviations of the intra-island dataset. The pink density function projecting to the left of each panel represents values observed in MORB from Gale et al. (2013). As in previous figures, red points represent lavas erupted on thin lithosphere and blue points represent lavas erupted on thick lithosphere. To the right of plots (a–c) are shown distribution functions for lavas on thin and thick lithosphere. Panels (d–f) show La/Sm vs isotopic ratios and panels (g-i) show Sm/Yb vs isotopic ratios. MORB values are shown in panels (d–i) as black dots. Lithospheric thickness is calculated assuming TLAB = 1225 °C. Detailed studies of isotopic composition spaces (e.g. Stracke, 2012) have shown that almost all of the data can be explained by mixing along an array between a Depleted Mantle (DM) component and a PREMA (Prevalent Mantle) or FOZO (Focus Zone) component, which in Sr–Nd–Hf space can be extended toward an Enriched Mantle (EM) component(s) (Supplementary DataFig. S5, available for downloading at http://www.petrology.oxfordjournals.org). Sun and Liang (2014) considered that the high frequency of OIB with 143Nd/144Nd∼0·5130 and 87Sr/86Sr∼0·7033 indicated that a mantle component with this signature (PREMA) was more likely than an overrepresentation of an intermediate value between DM and an EM component. However, our plots show that the high frequency mainly reflects the restricted composition space of the thin-plate OIB subset. On the other hand, there is some ambiguity as the data show that isotopic signatures trend back toward thin-plate values on the thickest lithosphere (>90 km), although more data should be accumulated to demonstrate statistical significance. Contrary to previous studies (Park, 1990; Niu et al., 2011), our plots of Pb isotopes show no lithospheric thickness effect (Supplementary DataFig. S6, available for downloading at http://www.petrology.oxfordjournals.org). Thin- and thick-plate OIB generally overlap, forming distributions clearly shifted to more enriched values compared to MORB. Exceptional island groups (e.g. Hawaii and Iceland) are noticeably MORB-like while other groups (e.g. Austral–Cook and St. Helena) adopt the opposite Pb space (called HIMU, for high 238U/204Pb). Harrison et al. (2017) studied Northwest Hawaiian Ridge lavas and also found no lithospheric thickness effect on Pb isotopes. Disequilibrium effects on isotopic signatures Because isotopic signatures in melts aggregated from heterogeneous sources depend on the relative contributions of trace element mass from each mantle component, disequilibrium may significantly affect signatures. The primary means of this effect is a difference in partitioning kinetics within different mantle components, which can be facilitated by differences in grain size, mineral modes, or composition. For example, a mantle component far from equilibrium (large grain size) will contribute a melt with a comparatively low concentration of incompatible elements. This will dampen signatures of disequilibrium components in low-degree melts without affecting their signatures in high-degree melts. Thus, the aggregate concentration will not only be a function of the volume and composition of each mantle component, but also the extent of equilibrium in each. Little is known about the details of thermodynamically-consistent coarsening in compositionally realistic mantle materials (e.g. Solomatov and Reese, 2008), but it may be reasonable to hypothesize that rates of grain coarsening over billions of years will differ among mantle components with significantly different major element compositions (e.g. fertile peridotite, depleted peridotite, recycled oceanic crust and metasomatic domains). Such lithologies are also likely to differ in mechanical properties, resulting in differences in deformation behavior, including grain size reduction and non-equilibrium (e.g. damaged) textures. Also, chemical and interfacial mobilities, as well as viscosities, are thought to depend on water content (e.g. Hirth & Kohlstedt, 1996; Hier-Majumder et al., 2005; Hosoya et al., 2005; Kubo et al., 2008). Thus, regional variations in water content could result in regional covariations in grain size. If water also covaries with other incompatible trace elements, these regions will act as mantle components with different extents of equilibrium during melting. To address such settings it will be important to develop new computational tools which address the details of textural evolution in mantle rocks. We have attempted a step in this direction and expect that future efforts linking more comprehensive analyses of the geochemical diversity of oceanic basalts with more physically explicative models of melting processes (e.g. Oliveira et al., 2018) and kinetics are likely to provide critical insights on the composition, dynamics and history of the mantle. CONCLUSION The lithospheric thickness effect on the chemistry of ocean island basalts is a major feature of the data. Because final melting depths can be inferred from seafloor ages using geophysical models, a careful examination of the data can provide critical constraints on the progression of melting processes in the mantle and the thermochemical properties of the melting region. We have presented equilibrium and disequilibrium models for decompression melting in fertile peridotite. Equilibrium models can fit the data if partition coefficients are allowed to differ substantially from experimental values, whereas disequilibrium models can fit the data simply with an appropriate grain size. Although the role of diffusion kinetics in melting processes is underdetermined, our results indicate that unless diffusion coefficients are significantly higher than values used here, disequilibrium may be a likely feature of mantle melting, even if it does not manifest precisely as our forward and inverse modelling results suggest. Inverse model predictions of LAB temperature and OIB source temperature are significant geophysical conclusions and do not depend on the role of disequilibrium. A LAB temperature of 1200–1300°C is in excellent agreement with geodynamic and rheological studies and suggests that geochemical data might be useful as a constraint on geophysical models of the oceanic upper mantle. Temperature in the source of OIB was predicted to be <1400°C, which is in the range of ambient mantle. Such a provocative conclusion will need to be tested by further inspections of more complete OIB databases. Because lava compositions represent evolved melts, whereas our models calculate aggregates of primary melts, it was necessary to correct bulk source compositions for crystallization effects to arrive at realistic estimates of source chemistry. Corrected models can fit a C1 chondrite mantle if the crystallization amount is large (75%), but only barely fit Primitive Mantle compositions for 20–50% crystallization amounts and cannot fit depleted mantle compositions at all. We posited the concept of a Mid-Ocean Ridge Island Basalt (MORIB) in order to examine how melting models which fit the OIB data can be used to predict compositions in OIB erupted on mid-ocean ridges, which can be compared to MORB. Although misfits between MORIB and MORB may be explained by differences in source composition or melting and crystallization processes, we found that misfits were in fact small. Such similarities are strong indications that sources and processes between the settings of OIB and MORB may be more similar than conventionally thought. Our analysis has shown the promise of using chemical signatures as a function of lithospheric thickness as an essential constraint on melting processes in the Earth. Although we examined the lithospheric thickness effect in isotopes, we did not include them in models. However, because REE data can be explained without reference to source heterogeneity and isotopic evidence seems to demand it, future models must incorporate isotopic signatures, although this task may not be trivial as new factors enter and existing factors increase in complexity. We have also shown that disequilibrium melting processes can be addressed in new levels of detail using lattice-based microdynamics methods. Nevertheless, our models remain incomplete views into the physics of decompression melting. The main shortcoming of our melting model is that microstructural evolution is not governed by a self-consistent, non-equilibrium thermodynamics framework. In the future this shortcoming will be remedied by a generalization of our approach for major elements. Acknowledgments C. J. Grose thanks Jeffrey G. Ryan for encouragement to investigate links between geophysical observations, thermodynamics and basalt chemistry. Thanks to Yaoling Niu for discussions on the lithospheric thickness effect in the early stages of this research and for providing his database of major, trace and radiogenic isotopes in OIB. Thanks to Paul Asimow and William L. Griffin for providing comments on early drafts. Yaoling Niu, Hugh O’Neill and two anonymous reviewers are thanked for critical reviews of the submitted manuscript. The editorial handling and commentary of Wendy Bohrson is greatly appreciated. This is contribution 1332 from the Australian Research Council Centre of Excellence for Core to Crust Fluid Systems (http://www.ccfs.mq.edu.au) and 1306 in the GEMOC Key Centre (http://www.gemoc.mq.edu.au). FUNDING C J. Grose acknowledges support from National Science Foundation (grant number 1826310). The work of J. C. Afonso has been supported by an Australian Research Council Discovery Grant (DP160103502). JCA also acknowledges support from the Research Council of Norway through its Centers of Excellence funding scheme, Project Number 223272. REFERENCES Adam J. , Green T. ( 2006 ). Trace element partitioning between mica- and amphibole- bearing garnet lherzolite and hydrous basanitic melt: 1. Experimental results and the investigation of controls on partitioning behavior . Contributions to Mineralogy and Petrology 152 , 1 – 17 . Google Scholar Crossref Search ADS Adam J. , Green T. ( 2010 ). Trace element partitioning between mica- and amphibole- bearing garnet lherzolite and hydrous basanitic melt: 2. Tasmanian Cainozoic basalts and the origins of intraplate basaltic magmas . Contributions to Mineralogy and Petrology 161 , 883 – 899 . Google Scholar Crossref Search ADS Adam J. , Green T. H. , Sie S. H. ( 1993 ). Proton microprobe determined partitioning of Rb, Sr, Ba, Y, Zr, Nb and Ta between experimentally produced amphiboles and silicate melts with variable F content . Chemical Geology 109 , 29 – 49 . Google Scholar Crossref Search ADS Adam J. , Green T. H. , Sie S. H. , Ryan C. G. ( 1997 ). Trace element partitioning between aqueous fluids, silicate melts and minerals . European Journal of Mineralogy 9 , 569 – 584 . Google Scholar Crossref Search ADS Adam J. , Locmelis M. , Afonso J. C. , Rushmer T. , Fiorentini M. L. ( 2014 ). The capacity of hydrous fluids to transport and fractionate incompatible elements and metals within the Earth’s mantle . Geochemistry Geophysics Geosystems 15 , 2241 – 2253 . Google Scholar Crossref Search ADS Afonso J. C. , Fernandez M. , Ranalli G. , Griffin W. L. , Connolly J. A. D. ( 2008a ). Integrated geophysical- petrological modeling of the lithosphere and sublithospheric upper mantle: Methodology and applications . Geochemistry Geophysics Geosystems 9 , 5 . Google Scholar Crossref Search ADS Afonso J. C. , Moorkamp M. , Fullea J. ( 2016 ). Imaging the lithosphere and upper mantle: Where we are at and where we are going. In: Moorkamp M. , Lelievre P. , Linde N. , Khan A. (eds) Wiley : AGU Monograph , in press. Afonso J. C. , Zlotnik S. , Fernandez M. ( 2008b ). Effects of compositional and rheological stratifications on small-scale convection under the oceans: Implications for the thickness of oceanic lithoshpere and seafloor flattening . Geophysical Research Letters 35 , 1 – 15 . Google Scholar Crossref Search ADS Anderson M. P. , Grest G. S. , Strolovitz D. J. ( 1989 ). Computer simulation of normal grain growth in three dimensions . Philosophical Magazine Part B , 59 , 293 – 329 . Google Scholar Crossref Search ADS Armienti P. , Gasperini D. ( 2007 ). Do we really need mantle components to define mantle composition? Journal of Petrology 48 , 693 – 709 . Google Scholar Crossref Search ADS Ballmer M. D. , Ito G. , van Hunen J. , Tackley P. J. ( 2011 ). Spatial and temporal variability in Hawaiian hotspot volcanism induced by small-scale convection . Nature Geoscience 4 , 457 – 460 . Google Scholar Crossref Search ADS Cabane H. , Laporte D. , Provost A. ( 2015 ). An experimental study of ostwald ripening of olivine and plagioclase in silicate melts: implications for the growth and size of crystals in magmas . Contributions to Mineralogy and Petrology 150 , 37 – 53 . Google Scholar Crossref Search ADS Canales J. P. , Carton H. , Carbotte S. M. , Mutter J. C. , Nedimovic M. R. , Xu M. , Aghaei O. , Marjanovic M. , Newman K. ( 2012 ). Network of off-axis melt bodies at the East Pacific Rise . Nature Geoscience 5 , 279 – 283 . Google Scholar Crossref Search ADS Carlson R. L. , Johnson H. P. ( 1994 ). On modeling the thermal evolution of the oceanic upper mantle: An assessment of the cooling plate model . Journal of Geophysical Research: Solid Earth 99 , 3201 – 3214 . Google Scholar Crossref Search ADS Carlson W. D. ( 2012 ). Rates and mechanism of Y, REE and Cr diffusion in garnet . American Mineralogist 97 , 1598 – 1618 . Google Scholar Crossref Search ADS Cherniak D. J. ( 2010 ). REE diffusion in olivine . American Mineralogist 95 , 362 – 368 . Google Scholar Crossref Search ADS Cmiral M. , Gerald J. D. F. , Faul U. H. , Green D. H. ( 1998 ). A close look at dihedral angles and melt geometry in olivine-basalt aggregates: A TEM study . Contributions to Mineralogy and Petrology 130 , 336 – 345 . Google Scholar Crossref Search ADS Coogan L. A. , O’Hara M. J. ( 2015 ). MORB differentiation: In situ crystallization in replenished-tapped magma chambers . Geochimica et Cosmochimica Acta 158 , 147 – 161 . Google Scholar Crossref Search ADS Crosby A. G. , McKenzie D. ( 2009 ). An analysis of young ocean depth, gravity and global residual topography . Geophysical Journal International 178 , 1198 – 1219 . Google Scholar Crossref Search ADS Crosby A. , McKenzie D. , Sclater J. G. ( 2006 ). The relationship between depth, age and gravity in the oceans . Geophysical Journal International 166 , 553 – 573 . Google Scholar Crossref Search ADS Dasgupta R. , Jackson M. G. , Lee C. A. ( 2010 ). Major element chemistry of ocean island basalts conditions of mantle melting and heterogeneity of mantle source . Earth and Planetary Science Letters 289 , 377 – 392 . Google Scholar Crossref Search ADS Davis E. E. , Lister C. R. B. ( 1974 ). Fundamentals of ridge crest topography . Earth and Planetary Science Letters 21 , 405 – 413 . Google Scholar Crossref Search ADS Doin M. P. , Fleitout L. ( 1996 ). Thermal evolution of the oceanic lithosphere: An alternative view . Earth and Planetary Science Letters 142 , 121 – 136 . Google Scholar Crossref Search ADS Ellam R. M. ( 1992 ). Lithosphere thickness as a control on basalt chemistry . Geology 20 , 153 – 156 . Google Scholar Crossref Search ADS Faul U. H. ( 2000 ). Constraints on the Melt Distribution in Anisotropic Polycrystalline Aggregates Undergoing Grain Growth. In: Bagdassarov N. , Laporte D. & Thompson A. B. (eds) Physics and Chemistry of Partially Molten Rocks. Petrology and Structural Geology , vol 11 . Dordrecht : Springer . Gale A. , Dalton C. A. , Langmuir C. H. , Su Y. , Schilling J.-G. ( 2013 ). The mean composition of ridge basalts . Geochemistry Geophysics Geosystems 14 , 30 . Google Scholar Crossref Search ADS Garapic G. , Faul U. H. , Brisson E. ( 2013 ). High-resolution imaging of the melt distribution in partially molten upper mantle rocks: evidence for wetted two-grain boundaries . Geochemistry Geophysics Geosystems , 14 , 556 – 566 . Google Scholar Crossref Search ADS Gast P. ( 1968 ). Trace element fractionation and the origin of tholeiitic and alkaline magma types . Geochimica et Cosmochimica Acta 32 , 1057 – 1086 . Google Scholar Crossref Search ADS Goutorbe B. , Hillier J. K. ( 2013 ). An integration to optimally constrain the thermal structure of oceanic lithosphere . Journal of Geophysical Research: Solid Earth 118 , 432 – 446 . Google Scholar Crossref Search ADS Grose C. J. ( 2012 ). Properties of oceanic lithosphere: revised plate model predictions . Earth and Planetary Science Letters 333–334 , 250 – 264 . Google Scholar Crossref Search ADS Grose C. J. , Afonso J. C. ( 2013 ). Comprehensive plate models for the thermal evolution of oceanic lithosphere . Geochemistry, Geophysics, Geosystems 14 , 3751 – 3778 . Google Scholar Crossref Search ADS Grose C. J. , Afonso J. C. ( 2015 ). The hydrothermal power of oceanic lithosphere . Solid Earth 6 , 1131 – 1155 . Google Scholar Crossref Search ADS Haase K. M. ( 1996 ). The relationship between the age of the lithosphere and the composition of the oceanic magmas: Constraints on partial melting, mantle sources and the thermal structure of plates . Earth Planet. Sci. Lett 144 , 75 – 92 . Google Scholar Crossref Search ADS Harrison L. N. , Weis D. , Garcia M. O. ( 2017 ). The link between Hawaiian mantle plume composition, magmatic flux and deep mantle geodynamics . Earth and Planetary Science Letters 463 , 298 – 309 . Google Scholar Crossref Search ADS Hasterok. ( 2013 ). A heat flow based cooling model for tectonic plates . Earth and Planetary Science Letters 311 , 386 – 395 . Herzberg C. ( 2004 ). Geodynamic information in peridotite petrology . Journal of Petrology 45 , 2507 – 2530 . Google Scholar Crossref Search ADS Hier-Mujumder S. , Anderson I. M. , Kohlstedt D. L. ( 2005 ). Influence of protons on Fe-Mg interdiffusion in olivine . Journal of Geophysical Research 110 , doi:10.1029/2004JB003292. Hillier J. K. , Watts A. B. ( 2005 ). Relationship between depth and age in the North Pacific Ocean . Journal of Geophysical Research 110 , B02405 . Google Scholar Crossref Search ADS Hirschmann M. M. , Kogiso T. , Baker M. B. , Stolper E. M. ( 2003 ). Alkalic magmas generated by partial melting of garnet pyroxenite . Geology 31 , 481 – 484 . Google Scholar Crossref Search ADS Hirth G. , Kohlstedt D. L. ( 1996 ). Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere . Earth and Planetary Science Letters 144 , 93 – 108 . Google Scholar Crossref Search ADS Hofmann A. W. ( 2006 ). Sampling mantle heterogeneity through oceanic basalts: Isotopes and trace elements. In: Holland, H. D. & Turekian, K. K. (eds) Treatise on Geochemistry , Vol. 2, Pergamon. pp. 1 – 44 . Hofmann A. W. , Hart S. R. ( 1978 ). An assessment of local and regional isotopic equilibrium in the mantle . Earth and Planetary Science Letters 38 , 44 – 62 . Google Scholar Crossref Search ADS Holm E. A. , Glazier J. A. , Srolovitz D. J. , Grest G. S. ( 1991 ). Effects of lattice anisotropy and temperature on domain growth in the two-dimensional Potts model . Physical Review A Atomic, Molecular and Optical Physics 43 , 2662 – 2668 . Google Scholar Crossref Search ADS Hosoya T. , Kubo T. , Ohtani E. , Sano A. , Funakoshi K. ( 2005 ). Water controls the field of metastable olivine in cold subducting slabs . Geophysical Research Letters 32 , doi:10.1029/2005GL023398. Hui L. , Guanghou W. , Feng D. , Xiufang B. , Pederiva F. ( 2003 ). Monte Carlo simulation of three-dimensional polycrystalline material . Materials Science and Engineering A357 , 153 – 158 . Google Scholar Crossref Search ADS Humphreys E. R. , Niu Y. ( 2009 ). On the composition of ocean island basalts (OIB): The effects of lithospheric thickness variation and mantle metasomatism . Lithos 112 , 118 – 136 . Google Scholar Crossref Search ADS Ito G. , Mahoney J. J. ( 2005a ). Flow and melting of a heterogeneous mantle: 1. Method and importance to the geochemistry of ocean island and mid-ocean ridge basalts . Earth and Planetary Science Letters 230 , 29 – 46 . Google Scholar Crossref Search ADS Ito G. , Mahoney J. J. ( 2005b ). Flow and melting of a heterogeneous mantle: 2. Implications for a chemically nonlayered mantle . Earth and Planetary Science Letters 230 , 47 – 63 . Google Scholar Crossref Search ADS Iwamori H. ( 1992 ). Melt-solid flow with diffusion-controlled chemical reaction . Geophysical Research Letters 19 , 309 – 312 . Google Scholar Crossref Search ADS Iwamori H. ( 1993 ). A model for disequilibrium mantle melting incorporating melt transport by porous and channel flows . Nature 366 , 734 – 737 . Google Scholar Crossref Search ADS Janssens K. G. F. , Raabe D. , Kozeschnik E. , Miodownik M. A. , Nestler B. ( 2007 ). Computational Materials Engineering: An Introduction to Microstructure Evolution . Cambridge : Elsevier Academic Press . Koepke J. , Behrens H. ( 2001 ). Trace element diffusion in andesitic melts: An application of synchrotron X-ray fluorescence analysis . Geochimica et Cosmochimica Acta 65 , 1481 – 1498 . Google Scholar Crossref Search ADS Kogiso T. , Hirschmann M. M. , Reiners P. W. ( 2004 ). Length scales of mantle heterogeneities and their relationship to ocean island basalt geochemistry . Geochimica et Cosmochimica Acta 68 , 345 – 360 . Google Scholar Crossref Search ADS Korenaga T. , Korenaga J. ( 2008 ). Subsidence of normal oceanic lithosphere, apparent thermal expansivity and seafloor flattening . Earth and Planetary Science Letters 268 , 41 – 51 . Google Scholar Crossref Search ADS Kubo T. , Ohtani E. , Kato T. , Kondo T. , Hosoya T. , Sano A. , Kikegawa T. ( 2008 ). Kinetics of the post-garnet transformation: Implications for density and rheology of subducting slabs . Physics of the Earth and Planetary Interiors 170 , 181 – 192 . Google Scholar Crossref Search ADS Liang Y. , Liu B. ( 2016 ). Simple models for disequilibrium fractional melting and batch melting with application to REE fractionation in abyssal peridotites . Geochimica et Cosmochimica Acta 173 , 181 – 197 . Google Scholar Crossref Search ADS Lyubetskaya T. , Korenaga J. ( 2008 ). Chemical composition of Earth’s primitive mantle and its variance: Method and results . Journal of Geophysical Reserch 112 , B03211 . Maclennan J. ( 2008 ). Concurrent mixing and cooling of melts under Iceland . Journal of Petrology 49 , 1931 – 1953 . Google Scholar Crossref Search ADS Mason J. K. , Lind J. , Li S. F. , Reed B. W. , Kumar M. ( 2015 ). Kinetics and anisotropyp of the Monte Carlo model of grain growth . Acta 82 , 155 – 166 . Mineralia, Marty J. C. , Cazenave A. ( 1989 ). Regional variations in subsidence rate of oceanic plates: a global analysis . Earth and Planetary Science Letters 94 , 301 – 315 . Google Scholar Crossref Search ADS McKenzie D. , Bickle M. J. ( 1988 ). The volume and composition of melt generated by extension of the lithosphere . Journal of Petrology 29 , 625 – 679 . Google Scholar Crossref Search ADS McKenzie D. , Jackson J. , Priestley K. ( 2005 ). Thermal structure of oceanic and continental lithosphere . Earth and Planetary Science Letters 233 , 337 – 349 . Google Scholar Crossref Search ADS McKenzie D. , O'Nions R. K. ( 1995 ). The source regions of Ocean Island Basalts . Journal of Petrology 36 , 133 – 159 . Google Scholar Crossref Search ADS Morgan W. J. ( 1972 ). Deep mantle convection plumes and plate motions . American Association of Petrology Geologists Bulletin 56 , 203213 . Müller R. D. , Sdrolias M. , Gaina C. , Roest W. R. ( 2008 ). Age, spreading rates and spreading symmetry of the world’s ocean crust . Geochemistry Geophysics Geosystems 9 , Q04006 . Google Scholar Crossref Search ADS Nishi M. , Kubo T. , Ohfuji H. , Kato T. , Nishihara Y. , Irifune T. ( 2013a ). Slow Si-Al interdiffusion in garnet and stagnation of subducting slabs . Earth and Planetary Science Letters 361 , 44 – 49 . Google Scholar Crossref Search ADS Nishi M. , Kubo T. , Ohfuji H. , Kato T. , Nishihara Y. , Irifune T. ( 2013b ). Slow Si-Al interdiffusion in garnet and stagnation of subducting slabs . Earth and Planetary Science Letters 361 , 44 – 49 . Google Scholar Crossref Search ADS Niu Y. , Langmuir C. H. , Kinzler R. J. ( 1997 ). The origin of abyssal peridotites: a new perspective . Earth and Planetary Science Letters 152 , 251 – 265 . Google Scholar Crossref Search ADS Niu Y. , Wilson M. , Humphreys E. R. , O'Hara M. J. ( 2011 ). The origin of intra-plate Ocean Island Basalts (OIB): the Lid Effect and its geodynamic implications . Journal of Petrology 52 , 1443 – 1468 . Google Scholar Crossref Search ADS Niu Y. , Wilson M. , Humphreys E. R. , OHara M. J. ( 2012 ). A trace element perspective on the source of ocean island basalts (OIB) and fate of subducted ocean crust (SOC) and mantle lithosphere (SML ). Episodes 35 , 310 – 327 . O’Neill H. S. C. ( 2016 ). The smoothness and shapes of chondrite-normalized rare earth element patterns in basalts . Journal of Petrology 57 , 1463 – 1508 . Google Scholar Crossref Search ADS O’Neill H. S. C. , Jenner F. E. ( 2012 ). The global pattern of trace-element distributions in ocean floor basalts . Nature 491 , 698 – 704 . Google Scholar Crossref Search ADS PubMed Oliveira B. , Afonso J. C. , Zlotnik S. , Diez P. ( 2018 ). Numerical modelling of multi-phase multi-component reactive transport in the earth’s interior . Geophysical Journal International 212 , 345 – 388 . Google Scholar Crossref Search ADS Park K.-H. ( 1990 ). Sr, Nd and Pb isotope studies of ocean island basalts: Constraints on their origin and evolution. PhD thesis, Columbia University, New York. Parsons B. , Sclater J. G. ( 1977 ). An analysis of the variation of ocean floor bathymetry and heat flow with age . Journal of Geophysical Research 82 , 803 – 827 . Google Scholar Crossref Search ADS Peate D. W , Hawkesworth C. J. ( 2005 ). U Series disequilibria: Insights into mantle melting and the timescales of magma differentiation . Reviews of Geophysics 43 , RG1003 . Google Scholar Crossref Search ADS Piazolo S. , Jessell M. W. ( 2008 ). Subgrain growth Potts Model. In: Bons P. D. , Koehn D. , Jessell M. W (eds) Microdynamics Simulation. Lecture Notes in Earth Sciences , vol. 106. Springer-Verlag Berlin Heidelberg, pp. 97 – 103 . Porter D. A. , Easterling K. E. ( 1981 ). Phase Transformations in Metals and Alloys . UK : Van Nostrand Reinhold Co. Ltd ., p. 445 . Qin Z. ( 1992 ). Disequilibrium partial melting model and its implications for trace element fractionations during mantle melting . Earth and Planetary Science Letters 112 , 75 – 90 . Google Scholar Crossref Search ADS Rudge J. F. , Maclennan J. , Stracke A. ( 2013 ). The geochemical consequences of mixing melts from a heterogeneous mantle . Geochimica et Cosmochimica Acta 114 , 112 – 143 . Google Scholar Crossref Search ADS Salters V. J. M. , Longhi J. E. , Bizimis M. ( 2002 ). Near mantle solidus trace element partitioning at pressures up to 3.4 GPa . Geochemistry Geophysics Geosystems 3 , 7 . Google Scholar Crossref Search ADS Sano J. , Ganguly J. , Hervig R. , Dohmen R. , Zhang X. ( 2011 ). Neodymium diffusion in orthopyroxene: Experimental studies and applications to geological and planetary problems . Geochimica et Cosmochemica Acta 75 , 4684 – 4698 . Google Scholar Crossref Search ADS Schneider M. E. , Egg.er D. H. ( 1986 ). Fluids in equilibrium with peridotite minerals: Implications for mantle metasomatism . Geochimica et Cosmochemica Acta 50 , 711 – 724 . Google Scholar Crossref Search ADS Sclater J. G. , Francheteau J. ( 1970 ). The implications of terrestrial heat flow observations on current tectonic and geochemical models of the crust and upper mantle of the earth . Geophysical Journal of the Royal Astronomical Society 20 , 509 – 542 . Google Scholar Crossref Search ADS Shaw D. M. ( 2006 ). Trace Elements in Magmas: A Theoretical Treatment . Cambridge : Cambridge University Press , pp. 243 . Shorttle O. ( 2015 ). Geochemical variability in MORB controlled by concurrent mixing and crystallisation . Earth and Planetary Science Letters 424 , 1 – 14 . Google Scholar Crossref Search ADS Sleep N. S. ( 2011 ). Small-scale convection beneath oceans and continents . Chinese Science Bulletin 56 , 1292 – 1317 . Google Scholar Crossref Search ADS Solferino G. F. D. , Golabek G. J. , Nimmo F. , Schmidt M. W. ( 2015 ). Fast grain growth of olivine in liquid Fe-S and the formation of pallasites with rounded olivine grains . Geochimica et Cosmochimica Acta 162 , 259 – 275 . Google Scholar Crossref Search ADS Solomatov V. S. , El-Khozondar R. , Tikare V. ( 2002 ). Grain size in the lower mantle: constraints from numerical modeling of grain growth in two-phase systems . Physics of the Earth and Planetary Interiors 129 , 265 – 282 . Google Scholar Crossref Search ADS Solomatov V. S. , Reese C. C. ( 2008 ). Grain size variations in the Earth’s mantle and the evolution of primordial chemical heterogeneities . J. Geophys. Res 113 , B07408 . Google Scholar Crossref Search ADS Spiegelman M. , Elliot T. ( 1993 ). Consequences of melt transport for uranium series disequilibrium in young lavas . Earth and Planetary Science Letters 118 , 1 – 20 . Google Scholar Crossref Search ADS Spiegelman M. , Kenyon P. ( 1992 ). The requirements for chemical disequilibrium during magma migration . Earth and Planetary Science Letters 109 , 611 – 620 . Google Scholar Crossref Search ADS Stein C. A. , Stein S. ( 1992 ). A model for the global variation in oceanic depth and heat flow with lithospheric age . Nature 359 , 123 – 129 . Google Scholar Crossref Search ADS Stracke A. ( 2012 ). Earth's heterogeneous mantle: A product of convection-driven interaction between crust and mantle . Chemical Geology 330–331 , 274 – 299 . Google Scholar Crossref Search ADS Stracke A. , Bizimis M. , Salters V. J. M. ( 2003 ). Recycling oceanic crust: Quantitative constraints . Geochemistry Geophysics Geosystems 8003 , 3 . Suhr G. , Seck H. A. , Shimizu N. , Gunther D. , Jenner G. ( 1998 ). Infiltration of refractory melts into the lowermost oceanic crust: evidence from dunite- and gabbro-hosted clinopyroxenes in the Bay of Islands Ophiolite . Contributions to Mineralogy and Petrology 131 , 136 – 154 . Google Scholar Crossref Search ADS Sun C. , Liang Y. ( 2014 ). An assessment of subsolidus re-equilibration on REE distribution among mantle minerals olivine, orthopyroxene, clinopyroxene, and garnet in peridotites . Chemical Geology 372 , 80 – 91 . Google Scholar Crossref Search ADS Sun S.-S. , McDonough W. F. ( 1989 ). Chemical and isotopic systematics of oceanic basalts: Implications for mantle composition and processes . Geological Society, London, Special Publications 42 , 313 – 345 . Google Scholar Crossref Search ADS Tikare V. ( 1995 ). Numerical simulation of grain growth in liquid phase sintered materials. Ph.D. thesis, Case Western Reserve University. Tikare V. , Holm E. A. , Fan D. , Chen L. Q. ( 1998 ). Comparison of phase-field and Potts models for coarsening processes . Acta Materialia 47 , 363 – 371 . Google Scholar Crossref Search ADS Van Mierlo W. L. , Langenhorst F. , Frost D. J. , Rubie D. C. ( 2013 ). Stagnation of subducting slabs in the transition zone due to slow diffusion in majoritic garnet . Nature Geoscience , DOI: 10.1038/NGEO1772. Van Orman J. A. , Grove T. L. , Shimizu N. ( 1998 ). Uranium and thorium diffusion in diopside . Earth and Planetary Science Letters 160 , 505 – 519 . Google Scholar Crossref Search ADS Van Orman J. A. , Grove T. L. , Shimizu N. ( 2001 ). Rare earth element diffusion in diopside: Influence of temperature, pressure and ionic radius and an elastic diffusion model for diffusion in silicates . Contributions to Mineralogy and Petrology 141 , 687 – 703 . Google Scholar Crossref Search ADS Van Orman J. A. , Grove T. L. , Shimizu N. ( 2002a ). Diffusive fractionation of trace elements during production and transport of melting in the earth’s upper mantle . Earth and Planetary Science Letters 198 , 93 – 112 . Google Scholar Crossref Search ADS Van Orman J. A. , Grove T. L. , Shimizu N. , Layne G. D. ( 2002b ). Rare earth element diffusion in a natural pyrope single crystal at 2.8 GPa . Contributions to Mineralogy and Petrology 142 , 416 – 424 . Google Scholar Crossref Search ADS Willbold M. , Stracke A. ( 2006 ). Trace element composition of mantle end-members: Implications for recycling of oceanic and upper and lower continental crust . Geochemistry Geophysics Geosystems , 7 , Q04004 . Google Scholar Crossref Search ADS Workman R. K. , Hart S. R. ( 2005 ). Major and trace element composition of the depleted MORB mantle (DMM) . Earth and Planetary Science Letters 231 , 53 – 72 . Google Scholar Crossref Search ADS Yang Z. , Sista S. , Elmer J. W. , Debroy T. ( 2000 ). Three dimensional Monte Carlo simulation of grain growth during GTA welding of titanium . Acta Materialia 48 , 4813 – 4825 . Google Scholar Crossref Search ADS Zhong S. , Ritzwoller M. , Shapiro N. , Landuyt W. , Huang J. , Wessel P. ( 2007 ). Bathymetry of the pacific plate and its implications for thermal evolution of the lithosphere and mantle dynamics . Journal of Geophysical Research 112 , B06412 . Zhu W. , Gaetani G. A. , Fusseis F. , Montesi L. G. J. , Carlo F. D. ( 2011 ). Microtomography of partially molten rocks: Three-dimensional melt distribution in mantle peridotite . Science 332 , 88 – 91 . Google Scholar Crossref Search ADS PubMed © The Author(s) 2019. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Chemical Disequilibria, Lithospheric Thickness, and the Source of Ocean Island Basalts
JF - Journal of Petrology
DO - 10.1093/petrology/egz012
DA - 2019-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/chemical-disequilibria-lithospheric-thickness-and-the-source-of-ocean-ddsGOGD6F3
SP - 755
VL - 60
IS - 4
DP - DeepDyve
ER -