Research ArticleGEOPHYSICS

Elastic properties of silicate melts: Implications for low velocity zones at the lithosphere-asthenosphere boundary

See allHide authors and affiliations

Science Advances  13 Dec 2017:
Vol. 3, no. 12, e1701312
DOI: 10.1126/sciadv.1701312


Low seismic velocity regions in the mantle and crust are commonly attributed to the presence of silicate melts. Determining melt volume and geometric distribution is fundamental to understanding planetary dynamics. We present a new model for seismic velocity reductions that accounts for the anomalous compressibility of silicate melt, rendering compressional wave velocities more sensitive to melt fraction and distribution than previous estimates. Forward modeling predicts comparable velocity reductions for compressional and shear waves for partially molten mantle, and for low velocity regions associated with the lithosphere-asthenosphere boundary (LAB), melt present at <5% distributed in near-textural equilibrium. These findings reconcile seismic observations for the LAB regionally and locally and favor models of strong coupling across the LAB rather than melt channeling due to shear deformation.


Seismic velocities constrain the physical, thermal, and chemical state of the Earth’s interior. Typically, velocities increase with depth where discontinuities are associated with marked compositional and/or phase changes. Of particular interest are regions where seismic velocities are slower than the global average, as has been documented crossing the mantle lithosphere-asthenosphere boundary (LAB). In the theory of plate tectonics, rigid plates move over convecting mantle with apparent minimal friction (1). The LAB by definition is the transition in thermal regime dominated by conduction above and advection below. However, recent high-resolution seismic studies show that the LAB can be sharp (24)—more consistent with a phase change, that is, melting, elastic softening, changes in grain size, or presence of fluids (for example, H2O, CO2, etc.) (58). Although a variety of mechanisms can lower seismic velocities, we consider the effects of silicate melt on the seismic properties of the mantle in the vicinity of the LAB drawing on our recent experimental work (9) that indicates that the elastic properties of amorphous silicates depend weakly on density, as well as by extension pressure, at upper mantle conditions. Our work implies a correspondingly weak pressure dependence for compressional wave velocities at least at pressures corresponding to the upper mantle for silicate melts (up to ~10 GPa). Here, we consider the consequences of these seismic properties for partially molten rock, in particular for the abundance and distribution of the silicate melt. We show that the anomalous compressional and shear wave velocities in the vicinity of the LAB are consistent with the presence of texturally equilibrated melt at melt fractions <0.05.

Commonly, the compressibility of silicate melts and their seismic properties at elevated pressures are predicted from an equation of state (EoS) such as the Birch-Murnaghan (BM) and Vinet formulations and extensions thereof (1012). These EoS implicitly adhere to Birch’s Law, which establishes a near-linear relation of the compressional wave velocity and density for crystalline materials of similar mean atomic weight (13). The obvious utility of Birch’s Law for seismology derives from the irrefutable inverse relationship between density and pressure for minerals and rocks that greatly facilitates predictions of compressional wave velocities with depth in the Earth. However, as discussed below, Birch’s law is not necessary appropriate for noncrystalline materials such as silicate glasses and melts.

By contrast, compressional wave velocities of many amorphous materials, specifically the silicate glasses and melts of interest here, exhibit a very weak and sometimes negative change in compressional wave velocity with increasing density (and pressure) (see Fig. 1) (9, 14, 15). This behavior is most pronounced at crustal and upper mantle pressures where changes in density are accommodated by distortion of the aluminosilicate network structure without changes in nearest-neighbor bond distances and/or coordination number (1618). That is, a change in volume (densification) is accommodated by the rotation of rigid polyhedra with little or no changes in the cation-oxygen bond length. This rotational freedom is far more restricted in a mineral phase because distortion of the crystal lattice has to occur mainly through changes in bond length. This unique mode of densification of silicate melts and glasses permits volume reductions that do not arise from the compression of interatomic bonds (that is, Si–O bonds). Densification accommodated by flexibility of the silicate network is not without limit because eventually atoms are so closely packed that further compression must involve bond shortening and coordination changes. However, it can account for anomalous compressibility and elastic properties of silicate liquids reported up to at least 10 GPa (19, 20), where increasing the density of amorphous silicates does not result in the corresponding increase in the elastic properties predicted on the basis of relationships developed for crystalline materials, as shown in Fig. 1. We refer to the low-energy process of densification via rotation of the rigid aluminosilicate polyhedral as network flexibility densification (NFD).

Fig. 1 Relationship between density and velocity for silicate crystals, glasses, and melts.

Dashed lines are lines of constant mean atomic weight (m); the molecular weight is divided by the number of atoms in the chemical formula. Birch’s Law states that for a constant m, there is a linear relationship between density and velocity. One-atm (1-atm) and high-pressure data for crystalline silicates follow Birch’s Law (dashed lines). One-atm data for glasses and melts are shown as labeled fields, with 1-atm data for melts shown as red diamonds [see the studies of Clark et al. (9) and Rivers and Carmichael (26) for further information]. High-pressure data (<6 to 10 GPa) for silicate glasses are shown as symbols with the compositions labeled (9, 15, 35, 36).


To explore the effects of NFD of silicate melts on the seismic properties of partially molten mantle, we draw on the work of Takei (11) that presents analytical solutions for seismic properties of aggregates having a range of melt fractions and geometric distributions. Takei provides a continuum solution for melt distribution, and here, we focus on three end-member melt geometries: (i) isolated spherules (unlikely, unless melt-solid equilibrium dihedral angles approach 90°), (ii) interconnected (cuspate) pores reflecting textural equilibrium for intermediate melt-solid dihedral angles, and (iii) elongate channels to subparallel films/sheets of melt expected as melt-solid dihedral angles approach zero or the aggregate is actively sheared (21, 22). The aspect ratio of the melt pockets (α = b/a) lie along a continuum where a value of unity corresponds to geometry 1, 0.1 corresponding to textural equilibrium (geometry 2), and 0.01 depicting the case for melt films or bands (geometry 3) (Fig. 2).

Fig. 2 Conceptual models of the melt geometry used in these calculations.

Geometry 1 is based on the studies of Berryman (37, 38). Geometry 2 is modified from Takei (11), Mavko (12), and Schmeling (39). Geometry 3 is modified from Holtzman and coworkers (6, 22). The effect of any intermediate geometry on seismic velocities can be modeled by modifying the aspect ratio (α) in Eqs. 3 and 4 [see the study of Takei (11) for further discussion of the derivation of α].

The reduction in compressional (P wave) and shear (S wave) velocities (ΔVP and ΔVS, respectively) as a function of melt fraction (φ) isEmbedded Image(1)Embedded Image(2)where V(0) is the velocity of the crystalline mantle, dV is the reduction the velocity, β is the ratio of the adiabatic bulk moduli (KS) of the solid to the liquid, γ is ratio of the shear modulus (G) to KS for the solid phases, and ρL and ρS are density of the silicate melt and crystalline mantle, respectively (11). ΛK and ΛG are geometric functions of pore shape and can be approximated byEmbedded Image(3)Embedded Image(4)where KF and N are the bulk and shear moduli, respectively, of the crystalline framework (11). KF and N are calculated for a crystalline matrix at a given φ, assuming the pore space is empty (the dry rock moduli) and the shape is described by α, as discussed above (1012). ΛK and ΛG are the slopes of KF(φ,α)/KS and N(φ,α)/G, respectively, as a function of φ and at relatively low values (<10%) can be approximated as a linear functions of φ. KF(φ,α)/KS and N(φ,α)/G are unity for melt fractions of zero and decrease with increasing porosity. With increasing connectivity (smaller α), both ΛK and ΛG increase (11). The main difference between Eqs. 1 and 2 is that ΔVP is dependent on β, whereas ΔVS is not. This is simply a reflection of the fact that compressional waves can travel through liquids, whereas shear waves cannot. Thus, in the implementation of Eqs. 1 and 2, it will be found that the choice of EoS for the liquids can have a potential consequence for ΔVP but does not influence ΔVS. Moreover, both ΔVP and ΔVS are proportional to φ/2, so the relative changes in P and S wave velocities (RSP) can be used to isolate the effect of melt distribution. RSP is found by dividing Eq. 2 to Eq. 1Embedded Image(5)

We use the ak135-f one-dimensional (1D) global average (23) to model the properties of the solid mantle because it does not impose a low velocity zone in the upper mantle (similar to ISAP91 reference model) and have self-consistent velocity and density models [similar to Preliminary Reference Earth Model (PREM)]. Note that the parameters in this model that depend on the 1D global average model (for example, γ and β) are the ratio of significantly different values, and so, the results presented here are not significantly influenced by variations in the 1D model values. For example, the parameter γ (G/KS of solid) is allowed to vary as a function of pressure as prescribed by the ak135-f model and falls in a narrow range from 0.48 to 0.54 up to 10 GPa, consistent with the γ values for calculated for PREM (see the Supplementary Materials for an expanded discussion of the comparison of 1D global average model values). For melt density, we consider two cases. The first being the more traditional approach using third-order BM EoS to model melt density, where the isothermal bulk moduli (KT) = 18 to 19 GPa, the pressure derivative (K′) = 4.5 to 5 GPa (24, 25), and KS = 18 GPa (26). The second is an empirical fit to the available high-temperature experimental density measurements for basalt (fig. S1). Regardless of EoS, the ratio ρLS is on average ~0.9 for upper mantle conditions.

Figure 3 compares the KS for the solid mantle according to the ak135-f model and the BM EoS and NFD models for silicate melt. Although KS for melts is markedly smaller than the solid mantle, the differences between the BM EoS and NFD models are striking. The similarities in the slope of BM EoS for the melt and for ak135-f arise from the implicit assumption of collinearity between density and elastic properties from Birch’s law, whereas the very weak pressure dependence of KS shown for the NFD model connects directly to the failure of silicate melts to obey Birch’s law, as discussed earlier and shown in Fig. 1. The impact of these differences in the pressure dependence of KS on seismic velocities of partially molten mantle can be further appreciated by comparing β for the BM EoS and NFD models in the insert of Fig. 3. β from the NFD model remains relatively constant between 0.5 and 10 GPa, whereas the BM EoS model predicts a rapid decrease in β between 0.5 and 5 GPa and a more modest decline at higher pressures.

Fig. 3 Bulk modulus (KS) as a function of depth for silicate melt and the crystalline mantle.

Silicate melts modeled using the NFD model (solid curve) and using the third-order BM EoS (dashed curve). KS for the solid is calculated from the ak135-f 1D global model (dotted curve) (23). Envelopes for the BM EoS model and the NFD model are in green and blue, respectively, and account for variations in density model and elastic properties for the melt phase as reported in the literature. (Inset) Ratio of the solid to liquid bulk moduli (β) as a function of depth for silicate melts modeled using the NFD model and using the third-order BM EoS. The thin curves correspond to variations due to density model (fig. S1) and elastic properties for the melt phase reported herein.

The determination of β provides sufficient constraint to solve Eqs. 1 and 5, permitting an assessment of the effects of melt fraction and geometry on reducing seismic velocities. For the following comparison of ΔVP, ΔVS, and RSP, the only difference in the calculations is how KS for the melt phase is modeled (Fig. 2). We first examine model results for ΔVP and ΔVS for melt fractions of 0.01 to 0.05 (Fig. 4 and figs. S4 to S6). As stated above, the controls of melt geometry (ΛK and ΛG) increase with increasing melt connectivity and are proportional to the velocity reduction, as shown in Eqs. 1 and 2. So, it follows that the magnitude of the velocity reduction for a given melt fraction at a given pressure is smallest for melts in isolated melt pockets (geometry 1) and greatest for melts aligned in films/sheets (geometry 3) (fig. S4). Although the absolute velocity reductions vary between the melt geometries, the general trends for ΔVP and ΔVS with pressure and melt fraction are similar for BM EoS and the NFD models. In both cases, ΔV for both P and S wave velocities increase with melt fraction. It can also be seen that regardless of the model, the effects of melt are greater on ΔVS than ΔVP. Likewise, for a given melt fraction, ΔVS is relatively constant over a wide range in pressure.

Fig. 4 Effect of varying melt fraction and geometry on seismic velocity reductions calculated as a function of depth.

(A) Reduction in seismic velocity (ΔV) for P and S waves (solid and dashed curves, respectively) as a function for pressure in equilibrated melt texture (geometry 2) for 0.01, 0.03, and 0.05 melt fraction shown as blue, green, and red curves, respectively, using the NFD model for the compressibility of the silicate liquid. (B) ΔV using BM EoS for the silicate liquid compressibility. (C) RSP calculated as a function of depth for the NFD model and BM EoS shown as black and gray curves, respectively. Melt geometry 1, 2, and 3 are shown as dotted, solid, and dashed curves, respectively. For clarity, the curves shown are calculated using the average value in Fig. 3. Curves with analysis for possible variations in density and KS for the melt phase are given in the Supplementary Materials. Boxes indicate the values of ΔV and RSP from recent seismic studies (24). The bold box for ΔV is the reported value, whereas the shaded box is for error. For seismic studies, the assumed transition width (thickness) of the low velocity layer affects the magnitude of ΔV (see the Supplementary Materials for geometries 1 and 3).

The most striking differences in the models can be seen for ΔVP. Above 0.5 GPa, ΔVP increases for the NFD model and decreases for the BM EoS model. This is due directly to the differences in bulk moduli shown in Fig. 3. Thus, melt fractions required to explain P wave velocity reductions at modest pressures by the BM EoS model are significantly greater than melt fractions required by the NFD model (for example, 20% at 3 GPa), and the deviation in estimated melt fractions increases at higher pressures (for example, 40% at 10 GPa).

Differences in ΔVP lead further to very different relationships between RSP and pressure. Figure 4C shows RSP calculated for melt geometries 1, 2, and 3 (Fig. 2) for both the NFD and BM EoS models. For both cases, melts isolated in pockets (geometry 1) yield the lowest RSP, whereas those aligned in sheets/films (geometry 3) have the highest values. For all geometries, RSP increases with pressure using the BM EoS model, whereas it remains relatively constant for the NFD model. Hence, RSP for the NFD model is most sensitive to assumptions about how the melt is distributed. The most marked changes occur because the melt distribution evolves from textural equilibrium (α = 0.1) to planar shear bands (α = 0.01). Note that NFD model results are consistent with experiments showing that RSP is close to unity at the onset of melting of peridotite (27, 28).


Figure 4 illustrates the marked effects and assumptions regarding the properties and distributions of melt can have seismic properties of partially molten rocks at upper mantle conditions. It is instructive to compare these predictions to observed seismic anomalies associated with the LAB (6, 22). Globally, S wave velocity anomalies associated with the LAB are typically 3 to 8% slower than those predicted by a 1D global average model. Determining melt fraction and geometry require both P and S wave velocity reductions, and so, we focus on the LAB below the western Pacific plate where both P and S wave reductions have been reported in the literature. Stern et al. (3) observed an 8 ± 2% P wave velocity reduction from their active source seismic study near New Zealand at 70 to 80 km (~2.5 GPa)—depths corresponding to the LAB. Similarly, receiver function studies below the western Pacific plate find comparable velocity reductions for S waves in the vicinity of the LAB (2, 4). These observations suggest that values of RSP are close to unity, which at LAB pressures (~2-3 GPa) is better accounted for by the NFD model than the BM EoS model. Moreover, at the depth of the LAB, this low RSP points to either isolated melt pockets (geometry 1) or texturally equilibrated melt (geometry 2) for the NFD model and cannot be achieved for any geometry using the BM EoS model (Fig. 4C). The failure of the BM EoS model to predict such a low RSP is simply a reflection of the fact that, under similar conditions, the assumption that the melt obeys Birch’s Law leads to P wave velocities for a partially molten aggregate that are markedly greater than the S wave velocities. In practical terms, what this means is that estimates for melt fraction made independently from P and S wave anomalies could differ by as much as 70% if the BM EoS model is used. What we regard as most significant is that, by considering the anomalous properties of the melt (NFD model), we can reconcile the similarities in the magnitude of the velocity reductions in the vicinity of the LAB below the western Pacific plate.

The volume and distribution of melt at the LAB has important implications for the nature of the LAB. Both geophysical and experimental studies have argued that melt is localized into subhorizontal channels along the LAB due to shear from the relative motions of the lithosphere and asthenosphere (1, 2, 6, 22). The large velocity reduction and sharp boundary (strong signal) at the LAB determined in receiver function studies have been used to argue for the presence of horizontally aligned melt channels (2). Likewise, Holtzman and coworkers (6, 21, 22) show that melt aligns in shear bands during LAB deformation of partial molten aggregates, approaching aspect ratios and melt connectivity depicted by geometry 3. For this melt geometry, 8% reductions in P and S wave velocity would imply very low melt fractions (fig. S4) and require P wave velocity reductions a factor of 2 higher than S wave velocity reductions (Fig. 4C). This is not what is observed seismically (Fig. 4), suggesting that if the velocity anomalies at the LAB are due to silicate melt, then the fraction of melt is higher than 1 to 2%. In these regions, seismic anisotropy arising melt shear bands can be excluded and rather indicates that seismic anisotropy arises from the deformation of the crystalline mantle [for example, in the study of Higgie and Tommasi (29)].

We can estimate melt fraction from the magnitude of the observed P and S wave velocity anomalies predicted for the NFD model in Fig. 4A. For example, the seismic velocities in the vicinity of the LAB below the western Pacific plate are reported to be 8 ± 2% slower than the global average (24), which according to the NFD model would correspond to melt fractions of 0.03 to 0.05. We regard that these are maximum melt fractions given that mineralogical variations (30), presences of volatile species (8), and anelasticity enhanced by high temperatures (7), among other factors, may contribute to a lowering seismic. These values exceed melt fractions expected to be trapped interstitially and may point to ponding of melt at the LAB, perhaps due to the competing effects of a melt-solid density contrast and low melt viscosity as proposed by Sakamaki et al. (31). A second line of evidence for appreciable melt in the vicinity of the LAB comes from studies of young (<10 Ma) petit-spot volcanism on old (~140 Ma) thick lithospheric plates on the western Pacific plate (32, 33) that seem to require a relatively extensive reservoir of partial melt immediately beneath the plate (34). If melts do pond at the LAB and are not localized there by the effects of shear deformation, then that may explain why the values of RSP are close to unity (reflecting textural equilibrium) rather than substantially larger as predicted if melts were segregated into subhorizontal channels. These conjectures are consistent with the recent work of Machida et al. (34) showing that the spatial progression of the petrogenically isolated magmatic systems feeding petit-spot volcanism is similar to the plate velocity implying that the lithosphere and asthenosphere are well coupled—at least in this tectonic setting.

In conclusion, we show that accounting for the anomalous density and elastic properties of silicate melts at upper mantle conditions can reconcile ΔVP, ΔVS, and RSP for seismic observations for the LAB, and further constrain melt fractions of 0.04 ± 0.01 for regions of the LAB exhibit P and S wave velocity reductions of 8 ± 2%. We demonstrate that ΔVP can provide a critical constraint on how melt is distributed at the LAB and that, together, the observed P and S wave velocity reductions suggest that melt assumes an equilibrium distribution. If this is true, then the lithospheric plate must be relatively well coupled to the underlying convecting asthenosphere despite the presence of melt.


Details on the calculations are given in the main text and the Supplementary Materials.


Supplementary material for this article is available at

fig. S1. Density of a mafic (basalt) melt as a function of pressure.

fig. S2. Comparison of velocity profiles for 1D global average models.

fig. S3. Comparison of model parameters for the ak135-f and PREM 1D global average reference models.

fig. S4. P and S wave velocity reductions for all melt geometries.

fig. S5. Analysis of variation in melt properties on the calculated P wave velocity reduction.

fig. S6. Analysis of variation in melt properties on the calculated S wave velocity reduction.

fig. S7. Analysis of variation in melt properties on the calculated RSP values.

Reference (40)

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 are grateful for numerous discussions during the development of this work with T.-R. A. Song and Y. Takei. We would also like to thank S. Petitgirard, the anonymous reviewers, and the editor for comments that helped to improve this paper. Funding: This research was supported in part by NSF grants EAR-1215714 to C.E.L. and a grant from the University of California Laboratory Fees Research Program (12-LR-237546) to C.E.L. C.E.L. also acknowledges support from the Danish National Research Foundation for the Niels Bohr Professorship at Aarhus University, Denmark. A.N.C. also acknowledges support from le Domaine d’Intérêt Majeur OxyMORE from la région d’Île-de-France, France and the NSF Division of Earth Sciences Postdoctoral Fellowship (award 1625205). Author contributions: A.N.C. and C.E.L. contributed equally to the manuscript. 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