## Abstract

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.

## INTRODUCTION

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 (*2*–*4*)—more consistent with a phase change, that is, melting, elastic softening, changes in grain size, or presence of fluids (for example, H_{2}O, CO_{2}, etc.) (*5*–*8*). 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 (*10*–*12*). 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 (*16*–*18*). 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).

## RESULTS

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).

The reduction in compressional (*P* wave) and shear (*S* wave) velocities (Δ*V*_{P} and Δ*V*_{S}, respectively) as a function of melt fraction (φ) is(1)(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 (*K*_{S}) of the solid to the liquid, γ is ratio of the shear modulus (*G*) to *K*_{S} 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 by(3)(4)where *K*_{F} and *N* are the bulk and shear moduli, respectively, of the crystalline framework (*11*). *K*_{F} 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 (*10*–*12*). Λ_{K} and Λ_{G} are the slopes of *K*_{F}(φ,α)/*K*_{S} and *N*(φ,α)/*G*, respectively, as a function of φ and at relatively low values (<10%) can be approximated as a linear functions of φ. *K*_{F}(φ,α)/*K*_{S} 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 Δ*V*_{P} is dependent on β, whereas Δ*V*_{S} 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 Δ*V*_{P} but does not influence Δ*V*_{S}. Moreover, both Δ*V*_{P} and Δ*V*_{S} are proportional to φ/2, so the relative changes in *P* and *S* wave velocities (*R*_{SP}) can be used to isolate the effect of melt distribution. *R*_{SP} is found by dividing Eq. 2 to Eq. 1(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*/*K*_{S} 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 (*K*_{T}) = 18 to 19 GPa, the pressure derivative (*K*′) = 4.5 to 5 GPa (*24*, *25*), and *K*_{S} = 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 ρ_{L}/ρ_{S} is on average ~0.9 for upper mantle conditions.

Figure 3 compares the *K*_{S} for the solid mantle according to the ak135-f model and the BM EoS and NFD models for silicate melt. Although *K*_{S} 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 *K*_{S} 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 *K*_{S} 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.

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 Δ*V*_{P}, Δ*V*_{S}, and *R*_{SP}, the only difference in the calculations is how *K*_{S} for the melt phase is modeled (Fig. 2). We first examine model results for Δ*V*_{P} and Δ*V*_{S} 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 Δ*V*_{P} and Δ*V*_{S} 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 Δ*V*_{S} than Δ*V*_{P}. Likewise, for a given melt fraction, Δ*V*_{S} is relatively constant over a wide range in pressure.

The most striking differences in the models can be seen for Δ*V*_{P}. Above 0.5 GPa, Δ*V*_{P} 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 Δ*V*_{P} lead further to very different relationships between *R*_{SP} and pressure. Figure 4C shows *R*_{SP} 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 *R*_{SP}, whereas those aligned in sheets/films (geometry 3) have the highest values. For all geometries, *R*_{SP} increases with pressure using the BM EoS model, whereas it remains relatively constant for the NFD model. Hence, *R*_{SP} 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 *R*_{SP} is close to unity at the onset of melting of peridotite (*27*, *28*).

## DISCUSSION

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 *R*_{SP} 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 *R*_{SP} 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 *R*_{SP} 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 (*2*–*4*), 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 *R*_{SP} 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 Δ*V*_{P}, Δ*V*_{S}, and *R*_{SP} 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 Δ*V*_{P} 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.

## MATERIALS AND METHODS

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

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/3/12/e1701312/DC1

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 *R*_{SP} 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.

## REFERENCES AND NOTES

**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.

- Copyright © 2017 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).