Research ArticleGEOCHEMISTRY

The prevalence of kilometer-scale heterogeneity in the source region of MORB upper mantle

See allHide authors and affiliations

Science Advances  22 Nov 2017:
Vol. 3, no. 11, e1701872
DOI: 10.1126/sciadv.1701872


The source regions of mid-ocean ridge basalts (MORB) are heterogeneous, consisting of chemically and lithologically distinct domains of variable size. Partial melting of such heterogeneous mantle sources gives rise to diverse isotopic compositions of MORB and abyssal peridotites. Variations in radiogenic isotope ratios in MORB are attributed to mixing of melts derived from enriched and depleted mantle components. However, melt mixing alone cannot fully account for the difference between the average 143Nd/144Nd in abyssal peridotites and their spatially associated MORB. We show that the more depleted Nd isotope composition in abyssal peridotites is a natural consequence of melt migration–induced mixing or smearing in the melting column. Sub-kilometer scale enriched mantle components or heterogeneities are significantly damped or homogenized in both the residue and erupted melt during their transit through the melting region. Heterogeneities with larger size and higher incompatible trace element abundance are more resistive to the mixing processes. The size-sensitive mixing depends on a parameter called the enrichment strength, which is the product of the heterogeneity size and the ratio between incompatible trace element abundance in the enriched and depleted mantle sources. Observed Nd-Hf isotope variations in MORB and abyssal peridotites can be reproduced if the enrichment strength is 20 to 60 km. These heterogeneities could be on the kilometer scale and have similar isotope ratios to but less incompatible trace element abundances than recycled oceanic crust.


Isotopic compositions of basalts and residual peridotites have been widely used to fingerprinting mantle source compositions (1). Because melting cannot fractionate isotopes of heavy elements (for example, 143Nd/144Nd or 176Hf/177Hf), the extracted instantaneous melt would have the same isotopic composition as its residue (2, 3). An important observation is that the average 143Nd/144Nd in abyssal peridotites is higher than their spatially associated mid-ocean ridge basalts (MORB) (referred to as the “MORB-AP offset” hereafter; Fig. 1A). In individual locations where abyssal peridotites are sampled, such offset also exists (fig. S1). The MORB-AP offset becomes more noticeable if refertilized samples are excluded (Fig. 1A). The offset of the mean 143Nd/144Nd in all abyssal peridotites relative to MORB has been attributed to the preferential melting of isotopically enriched lithology (4, 5) or higher Nd abundance in isotopically enriched melt (6). In the former case, radiogenic isotope ratio would correlate with the abundance of incompatible trace element in the residue. However, existing data do not demonstrate such a correlation at global or local scale (fig. S2). Instead, the lack of correlation between incompatible trace elements, such as Yb and 143Nd/144Nd in residual clinopyroxene, suggests that an isotopically distinct source could have experienced variable degrees of melting (fig. S2). In the latter case, residues inherit both the isotope ratio of enriched and depleted sources, whereas extracted melts are mixtures of the two sources. Because melts derived from the enriched lithology generally have higher incompatible trace element abundance, the isotope signature of the melt would be dominated by the enriched component. Although mixing has been widely used to explain variations in radiogenic isotope ratios observed in MORB and abyssal peridotites (1, 79), it is difficult to distinguish whether the mixing happens in the mantle source (for example, a mixed mantle source) before partial melting or during melt extraction and aggregation.

Fig. 1 Observations and the model for melting of an upwelling and chemically heterogeneous mantle beneath mid-ocean spreading centers.

(A) Histogram of 143Nd/144Nd in abyssal peridotites (AP) from Mid-Atlantic Ridge (37), Mid-Cayman Rise (37), and Southwest Indian Ridge (4, 6, 37); selected residual abyssal peridotites; and spatially associated MORB (see data compilation). (B) Cartoon illustrating the double-porosity ridge model for melting and melt migration in a heterogeneous mantle beneath a mid-ocean spreading center. The ridge model consists of bundles of double-porosity columns. Enriched heterogeneities are sketched as orange blobs in one column.

To distinguish where and when mixing takes place during melting and melt migration in the mantle, one should consider the physical size of mantle source heterogeneity and the rate of chemical exchange between the partial melt and its coexisting residue. During decompression melting of a spatially heterogeneous mantle, melts generated in the lower part of the melting column percolate upward and may encounter a residue of different isotopic composition. Chemical exchange then takes place, which tends to restore chemical and isotopic equilibrium between the melt and the residual solid (Fig. 2). Recent studies have suggested that small extent of chemical disequilibrium should be present during decompression melting beneath mid-ocean ridge spreading centers (1014) because the rates of Nd and Hf diffusion in pyroxene are relatively slow compared to the rate of melting (15, 16). The effective percolation velocity of a trace element in the partially molten mantle depends on its partition coefficient (17, 18). Chromatographic fractionation during melt percolation in the mantle can result in decoupling of Sr, Nd, and Hf isotope ratios in residual peridotites when there are spatial variations in these trace elements and isotopes in the mantle source (18). The extent of fractionation depends on the degree of chemical disequilibrium between the percolating melt and residual solid and the size of mantle source heterogeneity.

Fig. 2 Stretching and smearing of an enriched heterogeneity in an upwelling melting column.

Spatial variations of (A and B) 143Nd/144Nd and (C and D) 176Hf/177Hf at three selected times are shown by blue (residue) and red lines (matrix melt). The initial size of the heterogeneity is 4.4 km (A and C) or 1.1 km (B and D). The initial Nd and Hf in the heterogeneity are five times more enriched than in the depleted mantle. Pairs of facing triangles are material markers that track the locations of major element particles in the residue at the selected times. The small reduction in the marker distance is due to partial melting (fig. S3B). The degree of melting is linearly proportional to the vertical distance above the solidus (z = 0). The dashed lines mark the isotope ratio of the enriched endmember (see movie S1).

Here, we use a simplified double-porosity ridge model to investigate how the size and composition of mantle source heterogeneity affect the variation and correlation of Nd and Hf isotopes in basalts and residual peridotites during decompression melting of a two-component mantle beneath mid-ocean ridge spreading centers (Fig. 1B). The enriched component (19) has distinct isotope signature as well as higher concentrations of Nd and Hf than the depleted component (20, 21). We focus on the roles of the size and trace element abundances of heterogeneities under constraints of 143Nd/144Nd and 176Hf/177Hf systematics in MORB and abyssal peridotites. We are not considering Sr isotope in this study because Sr is less resistive to hydrothermal alteration than Nd and Hf in peridotites (22).


Stretching and smearing isotopic heterogeneity by melt migration in the mantle

A key result of the present study is that an enriched isotopic signal in the mantle source can be significantly deformed by melt migration–induced stretching and smearing during decompression melting in the mantle. This can be illustrated using a simple one-dimensional (1D) double-porosity model that treats the melting column as two overlapping continua consisting of low-porosity residual peridotites and high-porosity interconnected channels (2325). It has been suggested that replacive dunites observed in the mantle sections of ophiolites and massif peridotites may serve as high-porosity pathways for efficient melt extraction in the mantle (25). To model the spatial and temporal evolution of the enriched isotopic signal in the melting column, we use a time-dependent advection-reaction model that includes the effect of slow diffusion in pyroxene to calculate Nd and Hf concentrations and their isotopic ratios in the melts and the solids. To distinguish melts and solids in the two continua, we refer the low-porosity continuum as the matrix and the high-porosity continuum as the channel in the description below. During decompression melting, the porosity and the velocity of the matrix melt increase, whereas the velocity of the matrix solid decreases along the melting column (fig. S3), giving rise to increasing velocity difference between the matrix melt and the matrix solid. Without loss of generality, we consider spatial and temporal variations of two Gaussian-shaped isotopic heterogeneities (4.4 and 1.1 km long) embedded in the background mantle (Fig. 2). Nd and Hf abundances in the enriched source are five times more than in the background mantle.

After entering the melting column, the isotope signature in a layer of enriched heterogeneity splits into three distinct signals: one in the slow-moving matrix solid, one in the slightly faster matrix melt (Fig. 2), and one in the much faster channel melt (movie S1). As the enriched melt percolates upward, its size as marked by Nd and Hf isotope ratios increases due to the upward acceleration of the matrix melt, whereas its amplitude decreases due to the diffusive exchange between the melt and the residue in the matrix. Because the enriched melt percolates upward, it metasomatizes the residue immediately ahead of the original heterogeneity, increasing the size of the enriched signal in the matrix solid. However, the enriched signal in the residue has larger amplitude and lags behind that in the matrix melt [compare blue and red curves at 0.65 million years (My) in Fig. 2]. Hence, the residue could record more extreme enriched isotopic signals and larger isotope ranges than the matrix melt. Larger disequilibrium in Hf isotope between the matrix melt and matrix solid is expected (Fig. 2C) because the diffusion of Hf in pyroxene is slower than that of Nd (15, 16). The Hf isotopic signal in the residue is more resistive to amplitude reduction than Nd. The stretching- and smearing-induced chemical deformation (that is, shape change) of isotopic heterogeneity in the matrix melt and the residue is cumulative and accelerates upward in melting column (Fig. 2 and movie S1). This chemical deformation is more effective to disperse smaller heterogeneities than larger ones (Fig. 2). Depending on Nd and Hf abundances in their sources, sub-kilometer heterogeneities may not survive their transit through the melting column.

During near fractional melting, most of the melt generated by decompression melting is extracted through high-porosity channels. In contrast to that in the matrix melt, the enriched signal in the channel melt is quickly dispersed throughout the melting column (in less than 0.06 My) by the much faster flow of the channel melt (movie S1). Continuous extraction of matrix melts to the channel results in mixing between the depleted matrix melt and the enriched channel melt. The enriched isotope ratio in the channel is therefore gradually diluted upward. It is very difficult to preserve and sample extreme isotope compositions of the channel melt unless the melting column is very short. Therefore, isotopic signals of the channel melt and residual solid are decoupled. This has important implications for the observed MORB-AP offset.

The importance of the size and composition of mantle heterogeneity to the interpretation of the MORB-AP offset

To understand the MORB-AP offset, we first examine the effects of the size and composition of mantle heterogeneity on the isotope ratios in the residues and melts using the double-porosity ridge model (Fig. 1B). The triangular melting region consists of bundles of vertically upwelling double-porosity melting columns made of a low-porosity matrix and high-porosity channel. Melt generated in the matrix (blue regions in Fig. 1B) is extracted into nearby channels (green) and percolates upward until diverted and focused to the ridge axis along the base of the lithosphere where conductive cooling produces a permeability barrier (26, 27). The channel melt and the remaining matrix melt are collected below the permeability barrier where lateral melt transport and extensive mixing of melt take place (26). Melts derived from different melting columns are assumed to be well mixed to form pooled melt. The pooled melt is further mixed in a crustal magma chamber before eruption. The volume of the magma chamber is 1 km2 per unit length of the ridge axis, which is on the upper end of geophysical observations (28).

Figure 3 compares four cases of melting a two-component mantle that consists of variable proportion of depleted MORB mantle (DMM) and enriched mantle (EM). The isotopic ranges (19, 20) of EM and DMM are the same in the four cases. In each case, enriched heterogeneities of a specific size (2.2 or 8.9 km) are embedded randomly in the background DMM. The trace element abundance of EM is not well constraint. Here, we consider two cases: one is 3 times and the other is 12 times that of DMM (20). This enrichment factor is designated as “EMX” hereafter. The proportion of EM is adjusted so that the concentration-averaged isotope ratio of the source in each case fits the average of global MORB data. The concentration-averaged 143Nd/144Nd in the source, for example, is the total mass of 143Nd in the source divided by the total mass of 144Nd in the source. The calculated proportion of EM decreases from 11% to 3% as EMX increases from 3 to 12.

Fig. 3 The effect of the size and trace element abundance of mantle heterogeneity on the variations of isotope ratios in the pooled melt, erupted melt, and residues.

(A to D) The enrichment factor (EMX), proportion of EM in the source, and the initial size of enriched heterogeneities are listed above each panel. Gray outlines at the lower left and upper right corner are the range of isotope ratios of EM and DMM, respectively. Gray dots are global MORB data (39). The pentagram is the average of global MORB. Contours show the probability density function of the composition of the erupted melt predicted by the model. Histograms show the marginal distribution of 143Nd/143Nd and 176Hf/177Hf in the source, residue, pooled melt, and erupted melt. Dashed lines on histograms mark the concentration average of the source.

As demonstrated in the preceding section, melt migration–induced chemical deformation of a heterogeneous mantle source will produce a range of isotope ratios in the residue. If the size of heterogeneity and the enrichment factor are large (for example, 8.9 km and EMX = 12; Fig. 3D), the residue presents a wide range of isotope ratios in between DMM and EM in the source. The residue has little modification of the isotopic composition compared to the source. Because DMM has more than 96% volume fraction in the source, the most common isotope signature of the residue is that of DMM. If the size of heterogeneity and the enrichment factor are small (for example, 2.2 km and EMX = 3; Fig. 3A), the 143Nd/144Nd of the residue no longer resembles the source. Melt migration–induced smearing or mixing has produced a significant amount of the residue with intermediate 143Nd/144Nd, which converges to the concentration-averaged value of the source. Increasing Nd abundance in EM makes the enriched source more resistive to mixing and thus favors the preservation of the original mantle source 143Nd/144Nd in the residue (compare Fig. 3, A and C, or Fig. 3, B and D). Mixing of 176Hf/177Hf in the residue also occurs for the small and less enriched heterogeneities (Fig. 3A) but at a smaller extent compared to 143Nd/144Nd. This is due to the slower diffusion of Hf in pyroxene (Fig. 2, B and D). Thus, the most common 176Hf/177Hf in the residue does not converge to the concentration average of the source. The residue has the most common 176Hf/177Hf in between DMM and the concentration average of the source.

Mixing of partial melts derived from a wide range of depth takes place in the high-porosity channel network during melt migration in the mantle and accumulation of melts in the magma chamber process (Fig. 1B). Increasing the abundance of trace element in the enriched heterogeneity could partly offset the dilution effect. Because the chemical signal in the matrix melt is already damped, the composition of channel melt will inherit the size-sensitive mixing effect in the matrix. Because the channel melt is a mixture of all melts extracted from the matrix, all enriched chemical signals in a melting column are integrated in the channel. As the size of heterogeneity increases, the isotopic enrichment in the channel would be less diluted (movie S1). Consequently, the channel melt will be mixed less and samples a wider range of isotopic variability if the size of heterogeneity is larger.

The pooled melt mixes channel melt from all melting columns in the triangular melting region. Thus, the composition of pooled melt will also inherit the size-sensitive mixing effect and be further mixed under the permeability barrier. As the size and the enrichment decrease, mixing and dilution in the matrix and channel melts strengthens, which makes the isotopic composition of pooled melt to converge to the concentration average of the source. The magma chamber process leads to further homogenization, but this effect is very small due to the small magma chamber residence time, which is defined as the ratio between the magma chamber volume and the melt volume flux (29). In Nd-Hf isotopic space, the erupted melt has a wide range of composition around the concentration average of the source or the average MORB in all cases unless the enrichment is high and the size is large (Fig. 3D). The range of melt composition shrinks as the size of heterogeneity decreases. Increasing the trace element abundance of EM or EMX could enlarge the range of compositions in the erupted melt. In all cases, the erupted melt does not represent either the EM or the DMM source. The entire field of EM and part of the depleted isotope ratios of DMM are underrepresented in the erupted melt.

The relation between the isotope ratios in the erupted melt and the residue changes systemically with the size and the enrichment factor. In a case of large size and high enrichment (8.9 km and EMX = 11.9; Fig. 3D), the erupted malt peaks at roughly the same 143Nd/144Nd as the residue. Hence, there is no MORB-AP offset in this case. Because the size or the enrichment decreases, the composition of the erupted melt converges to the concentration average, whereas the composition of the residue peaks at DMM-like 143Nd/144Nd. The MORB-AP offset is produced. When the size is small and the enrichment of EM is low (2.2 km and EMX = 3.04; Fig. 3A), both the melt and residue converge to the source-averaged 143Nd/144Nd. The MORB-AP offset disappears. The most common 176Hf/177Hf in the residue is still higher than the erupted melt (Fig. 3A). Therefore, the 143Nd/144Nd MORB-AP offset as observed in natural samples should correspond to a set of moderate size and enrichment. Because the Hf offset is observed more often than the Nd offset in our simulations, the residue could be plotted above the melt trend in Nd-Hf isotopic space (fig. S7).

There is a trade-off between the size (L) and the enrichment factor (EMX) of EM in producing the Nd and Hf isotopic variations presented in the preceding section: The probability density function (PDF) of the erupted melt and the residue will not change so long as the product L*EMX remains the same. This is demonstrated in Fig. 3 (B and C) (L*EMX = 27 km and 26 km, respectively). We refer the product, L*EMX, as the enrichment strength. Figure 4 (A and B) summarizes our simulation results for the most common 143Nd/144Nd in the residue and the erupted melt (color-coded dots) for a wide range of EM size and enrichment factor. The range or variability of 143Nd/144Nd (marked by vertical bars) is defined by half-height width of the PDF in the residue or erupted melt. If the enrichment strength is larger than 30 km (as in the case of Fig. 3D), the composition of the residue is almost the same as DMM. As the enrichment strength decreases to 10 km, by reducing either the size or the Nd abundance in EM, the enriched signals are smeared by melt migration–induced mixing, the marginal PDF becomes broader, and the peak location migrates to the concentration average of the source. As the enrichment strength decreases below 10 km (as in the case of Fig. 3A), strong mixing in the residue leads to narrowing of the marginal PDF. For the erupted melt, reducing the enrichment strength from 200 to 1.8 km shrinks the marginal PDF (Fig. 4B). The 143Nd/144Nd of EM (lower than 0.51295) is underrepresented in the marginal PDF of the erupted melt. The average DMM-like 143Nd/144Nd (0.5132) is within the half-height width of the marginal PDF of the erupted melt only if the enrichment strength is larger than 80 km. At such large enrichment strength, the marginal PDF of the erupted melt centers at a value larger than the concentration average of the source. If the enrichment strength is smaller than 80 km, the marginal PDF of melt is symmetric and centers around the source average.

Fig. 4 Size-sensitive mixing for the residue and the melt.

The 143Nd/144Nd of the residue (A) and the erupted melt (B) as a function of the enrichment strength of the heterogeneity (L*EMX). Colored dots are the 143Nd/144Nd at the peak of marginal distribution in the residue or the erupted melt. Vertical bars are the half-height width of the marginal distribution of 143Nd/144Nd. The thick gray lines are the concentration average of the source. Dashed lines in (B) are the range of MORB. If the enrichment strength is smaller than 20 km [arrow in (A)], the peak of marginal distribution in the residue would converge to the source average, and no offset would be observed. If the enrichment strength is larger than 60 km [arrow in (B)], the range of isotopic compositions in the erupted melt would be larger than that observed in the MORB.

The observed MORB-AP offset requires that the 143Nd/144Nd of melt is relatively well mixed to represent the concentration average of EM and DMM, whereas the 143Nd/144Nd of the residue is not efficiently mixed so that the residue still has a DMM-like signature. To produce an MORB-AP offset (0.51306 versus 0.51320), the enrichment strength should be greater than 20 km. The abyssal peridotites seem to sample a wider range of 143Nd/144Nd than the locally associated MORB. The range of observed MORB composition requires an enrichment strength of 60 km. The range of melt composition will expand if the size of the magma chamber is smaller than used here or the mixing during melt pooling is not complete (see the next section). Therefore, the enrichment strength of 60 km should be taken as an upper bound. The lower bound (20 km) of the enrichment strength for preserving the MORB-AP offset varies only two times, and the upper bound (60 km) for the restricted range of melt composition is invariable for a range of model parameters (fig. S5).


In our model, the pooled melt is a well-mixed aggregation of melts coming out of each melting column in the triangular melting region. This is likely an oversimplification of the mixing process during melt focusing because the isotopic variability of melt inclusions suggests incomplete mixing of mantle-derived melts (30, 31). We test a case in which melt focusing does not mix melts derived from vertical melt channels. The erupted melt is a collection of melts with different compositions before entering the magma chamber. We can sample the collection at a probability proportional to the melt flux to obtain a statistical representation of the melt collection (fig. S6). As expected, the composition of the melt collection is more variable than the well-mixed pooled melt. The lowest degree (2%) channel melt in the melt collection is barely diluted by melt extraction and contributes to the most depleted isotope ratios in the source that pooled melt cannot sample. Because of the small volume proportion of this low degree melt in all melts generated in the melting region, the chance of sampling extreme isotope ratios in the melt collection is very low. The marginal PDF of 143Nd/144Nd in the melt collection tends to peak at more depleted isotope ratios than the pooled melt (fig. S6). For the case shown in fig. S6, the peak of 143Nd/144Nd in the melt collection is indistinguishable from the residue, whereas the pooled melt has noticeably lower 143Nd/144Nd than the residue. To produce the MORB-AP offset, some mixing during melt focusing or in the magma chamber is required. The relatively low viscosity of basaltic melt favors rapid overturn and efficient homogenization in the magma chamber (32). The erupted melt is sufficiently, if not completely, mixed to approach the concentration-averaged isotope ratios of the source. Because inefficient mixing during melt pooling makes it more difficult to produce the MORB-AP offset, the constraints for the enrichment strength could only become tighter.

The 20- to 60-km enrichment strength and the form of enriched heterogeneities

The Nd-Hf isotopic composition of the heterogeneity or EM is related to the unique elemental ratios in oceanic crust, including a small portion of sediments (19, 33). The abundance of Nd in the recycled oceanic crust is 25 to 12 times higher than DMM, depending on the proportion of sediments, the age of recycling, the extent of subduction modification, and the choice of DMM composition (19, 34). The proportion of oceanic crust should be 3% to fit the average MORB (Fig. 3, C and D). The subducted oceanic crust is likely thinned or stretched by mantle convection, resulting in a distribution of size and shape of the heterogeneity (35). A significant amount of recycled crustal materials would melt before the onset of peridotite melting (36), producing melt with enriched incompatible trace element abundance and isotope ratios. Such enriched melt could metasomatize the overlying depleted mantle, producing heterogeneities with a range of EMX. The isotope ratio of this metasomatized mantle is expected to be dominated by the recycled materials. The size of these metasomatized heterogeneities is expected to be larger than unmelted crustal materials because of percolation of small-degree melts derived from them. Because the extent of melt migration–induced mixing after the onset of peridotite melting is scaled with the enrichment strength and the observed MORB-AP offset requires a moderate extent of mixing, most of the heterogeneities cannot have an enrichment strength outside the range of 20 to 60 km. The enrichment factor of EM could be in the range of 1 to 25. Therefore, the size of heterogeneities responsible for the observed isotope ratios in MORB and abyssal peridotites could be on the kilometer scale. Our model does not imply the absence of smaller heterogeneities in the mantle source as the sub-kilometer class does not contribute to the observed MORB-AP offset. The size and composition of mantle heterogeneity are important to the interpretation of isotope and incompatible trace element data in basalts and residual peridotites and must be taken into consideration in geochemical studies of mantle-derived rocks. Systematic studies of high-resolution spatially correlated mantle samples may help to constrain the size, distribution, and nature of chemical heterogeneities in the mantle.


Data compilation

We have compiled 96 published 143Nd/144Nd data of abyssal peridotites from Mid-Atlantic Ridge (37), Southwest Indian Ridge (4, 6, 37), and Mid-Cayman Rise (37). If the TiO2 content in the spinel is larger than 0.12 weight % or the La/Ce ratio normalized by CI chondrite (38) is greater than 2, the sample was considered refertilized or metasomatized (13, 37). The remaining 67 samples were referred to as residual abyssal peridotites. We used a geochemical database, PetDB (Petrological Database) (39), to choose any MORB data as the local MORB along the ridge within 550 km of the dredge site of each residual abyssal peridotite. The data were summarized in fig. S1. Figure 1A combines all abyssal peridotites and their local MORB. To account for the different size of residual abyssal peridotite data set and local MORB data set, the frequency of local MORB was weighted by the number of local residual abyssal peridotite samples.

Model summary

In Fig. 2, we modeled the spatial and temporal variations of 143Nd/144Nd and 176Hf/177Hf heterogeneity during concurrent melting and melt migration using a 1D disequilibrium double-porosity model, which is built on previous studies (14, 23, 24). We assumed that the physical state of the system is at steady state. Darcy’s law constrains the velocity difference between the melt and solid residue. The diffusion of Nd and Hf in clinopyroxene limits the rate of chemical exchange between the matrix melt and the residue (18). The time-dependent mass conservation equations for 143Nd, 144Nd, 176Hf, and 177Hf were solved using a Runge-Kutta discontinuous Galerkin method (40). The error for an isotope was 10−9 parts per million (ppm). The time-dependent boundary condition has a Gaussian-shaped enrichment with a background value of DMM and a peak value of a specific EM. In Figs. 3 and 4, the ridge model consisted of nine melting columns with the maximum degree of melting 2%, 4%, …, 18% stacked together as half of the triangular melting region. Gaussian peaks of EM signature were placed randomly in the ambient DMM. The composition of DMM varied within extreme values from Salters and Stracke (20) and Workman and Hart (21). The isotope ratios of EM were related to the age of recycled oceanic crust (1 to 3 billion years ago) and the proportion of sediments in the crust (0 to 10%). The abundance of trace element in EM was determined by mixing recycled oceanic crust and DMM. The range of trace element abundance in EM was from 3 to 12 times of DMM. The DMM composition was from a bivariate normal distribution within the range of extreme values (20, 21). In the ridge model, all residues with degrees of melting between 5% and 18%, equivalent to the degree of melting of abyssal peridotites (14), were included. Channel melt and matrix melt from all melting columns were assumed to mix instantaneously below the permeability barrier and form the pooled melt. The composition of the magma chamber was modeled using pooled melt as input (29). See the Supplementary Materials for the detailed model setup and the numerical method.


Supplementary material for this article is available at

table S1. Errors and numerical orders.

table S2. Physical parameters in working equations (Eqs. 33 to 35 and 41 to 43).

fig. S1. Histograms of Nd isotopes in abyssal peridotites, residual abyssal peridotites, and local MORB at individual locations.

fig. S2. Lack of correlation between the REE concentration and 143Nd/144Nd of clinopyroxenes in abyssal peridotites.

fig. S3. Spatial variations of porosity, residue, and matrix melt velocities in the 1D upwelling melting column.

fig. S4. The 176Hf/177Hf of the residue and the erupted melt as a function of the enrichment strength of the heterogeneity (L*EMX).

fig. S5. The effect of chemical exchange rate and melt extraction on 143Nd/144Nd of the residue and the erupted melt.

fig. S6. The composition of melt collection, residues, pooled melt, and sources in the Hf-Nd isotope correlation diagram.

fig. S7. A comparison of the predicted compositions of the residue and erupted melt in a case with Hf offset but without Nd offset.

movie S1. Spatial and temporal variations of 143Nd/144Nd and 176Hf/177Hf in the melting column.

References (4153)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank C. Huber, S. Mallick, and A. Saal for useful discussion. This paper benefited from thoughtful reviews by G. Ito, M. Jackson, and the associate editor C.-T. Lee. Funding: This work was supported by a grant from the U.S. NSF (EAR-1624516). Author contributions: Y.L. conceived the project and secured its funding. B.L. and Y.L. formulated the problem, developed the model, and wrote the manuscript. B.L. developed the numerical scheme and carried out the modeling. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article