Significance Although plate tectonics has seen broad acceptance for Earth, the manner in which lithospheric plates are coupled to Earth’s deeper interior is still heavily debated. In particular, recent seismological observations suggest a sharp, flat base of the lithosphere, whereas thermal models suggest a gradational boundary that deepens with age. Based on laboratory experiments, we suggest that thermal models are most appropriate and that seismic studies are detecting features frozen into the lithosphere after melting at midocean ridges. Experiments on olivine aggregates demonstrate that the seismic characteristics of deforming upper mantle are dramatically different between melt-free and low-melt-fraction aggregates. A model of upper-mantle flow incorporating these results predicts seismological features in excellent agreement with observations beneath the Pacific Ocean basin.

Abstract Tectonic plates are a key feature of Earth’s structure, and their behavior and dynamics are fundamental drivers in a wide range of large-scale processes. The operation of plate tectonics, in general, depends intimately on the manner in which lithospheric plates couple to the convecting interior. Current debate centers on whether the transition from rigid lithosphere to flowing asthenosphere relates to increases in temperature or to changes in composition such as the presence of a small amount of melt or an increase in water content below a specified depth. Thus, the manner in which the rigid lithosphere couples to the flowing asthenosphere is currently unclear. Here we present results from laboratory-based torsion experiments on olivine aggregates with and without melt, yielding an improved database describing the crystallographic alignment of olivine grains. We combine this database with a flow model for oceanic upper mantle to predict the structure of the seismic anisotropy beneath ocean basins. Agreement between our model and seismological observations supports the view that the base of the lithosphere is thermally controlled. This model additionally supports the idea that discontinuities in velocity and anisotropy, often assumed to be the base of the lithosphere, are, instead, intralithospheric features reflecting a compositional boundary established at midocean ridges, not a rheological boundary.

Although plate tectonics is the unifying paradigm in the Earth sciences, important questions remain regarding the physical nature of a tectonic plate. A cornerstone of plate tectonic theory states that plates translate across Earth’s surface in a relatively rigid and coherent fashion, with deformation largely concentrated at plate boundaries. Restated, the plates are taken to have a high viscosity with a relatively sharp transition to the less viscous convecting mantle beneath. Referring to plates as the lithosphere and to the underlying rock as the asthenosphere [terminology which predates plate tectonic theory (1)] is now common. Partial decoupling of the asthenosphere from the lithosphere appears essential to plate-tectonic-like behavior (2), but whether a change in material properties at the base of the lithosphere arises from increasing temperature (3) or from transitions in melt (4⇓⇓–7) or water content (8, 9) remains unclear, even in the tectonically simple ocean basins.

Recent seismological investigations of oceanic upper mantle indicate a change in composition at depths often associated with the lithosphere−asthenosphere boundary (LAB). Receiver-function studies and underside reflections of SS precursors detect a sharp velocity discontinuity (the Gutenberg discontinuity) at a depth of 40 to 150 km (4, 6, 10⇓–12). The sharpness of this discontinuity may indicate the presence of a small amount of melt (4), and the lack of a strong dependence of the discontinuity depth on plate age (6, 10, 13) is roughly consistent with a carbonated peridotite fluid (7, 14). Alternatively, a transition from water-poor to water-rich conditions could establish a sharp reduction in seismic velocity. Extrapolations from laboratory experiments (9) suggest that a dehydration front established by melting near the ridge axis can produce a water-poor layer whose thickness would not depend on plate age (8). In both interpretations, the sharp drop in velocity is often taken to indicate a change in viscosity structure and therefore the base of the lithosphere.

Instead of compositional variations, recent investigations of seismic anisotropy (13, 15, 16) provide evidence for a lithospheric base that is thermally controlled and measurably deeper than the sharp seismic discontinuities. Beghein et al. (13) highlighted two observations. First, the azimuth of the seismically fastest direction, ψ, appears to be well aligned with the plate motion direction at depths of >100 km, consistent with motion of the overlying lithosphere accommodated through viscous flow of peridotite. In shallower portions of the system, ψ is at a high angle to the plate motion direction, suggesting that plate motion direction has changed since it was “frozen in” during thickening of the lithosphere. Importantly, the depth to which ψ is well aligned with the plate motion direction appears to depend on plate age, roughly in accord with standard thermal models for oceanic upper mantle. Based on these motivating observations, we test the previously suggested hypothesis (13, 15, 17, 18) that the well-aligned region represents a relatively low-viscosity asthenosphere, whose upper boundary is thermally controlled.

Second, recent seismological studies (13, 15) suggest a transition in the magnitude of radial anisotropy (ξ, the ratio between the velocity of horizontally polarized S waves and vertically polarized S waves) as a function of depth. Although these surface-wave studies have not included the high-frequency data required to fully resolve anisotropy at the shallowest depths, a transition from low values of ξ at shallow depths to higher values of ξ at greater depths appears to occur near 80 km. The depth of this change in ξ exhibits very little—if any—dependence of depth on plate age and correlates well with estimates for the location of the Gutenberg discontinuity (13). This discontinuity in ξ has been attributed to ubiquitous melt structures in a partially molten asthenosphere (13, 19), but that hypothesis does not explain the presence of strong azimuthal anisotropy described above unless melt-rich layers are significantly inclined from horizontal (5). Alternatively, this discontinuity has been attributed to processes occurring near midocean ridges (13, 20). We support the latter hypothesis by specifically suggesting that this discontinuity results from the crystallographic alignment of olivine grains, which forms in partially molten rocks at shallow depths beneath the ridge axis. This process imparts a seismologically observable weakening in the alignment of the fast axes of olivine, reducing the magnitude of anisotropy, but does not dramatically modify the viscosity structure of the mantle.

Laboratory Constraints on Upper-Mantle Anisotropy To test these hypotheses, we draw on the results from a combination of new (Dataset S1) and recently published (21, 22) laboratory deformation experiments conducted on (i) nominally melt-free (<1% melt) and (ii) melt-present (>1% melt) olivine aggregates. These torsion experiments (Supporting Information), which reached shear strains as high as 18, provide insight into the evolution of olivine textures, likely the primary source of upper-mantle anisotropy (23). This data set is unique in that most previous estimates of the evolution of anisotropy have been based on experiments conducted in transpression (24, 25) rather than simple shear, as achieved here. At high shear strain (>10), all samples exhibit steady-state textures (Fig. 1). Nominally melt-free samples have [100] dominantly parallel to the shear direction and [010] dominantly normal to the shear plane, consistent with previous investigations of experimentally (21, 22, 24, 25) and naturally (26, 27) deformed olivine (Fig. 1A). Melt-present samples exhibit weaker textures than nominally melt-free samples, with [100] and [001] distributed in a horizontal girdle, similar to previous experimental results (28, 29) (Fig. 1B). Fig. 1. Characteristic olivine textures in experimental samples. Olivine textures are presented for (A) nominally melt-free aggregates and (B) aggregates of olivine + 4% basalt. Samples were deformed in torsion to approximately the same shear strain. Each data point corresponds to the orientation of one grain, and each pole figure depicts 1,000 individual grain orientations. Data points are colored by the corresponding value of the orientation distribution function. The addition of a basaltic melt weakens the overall texture, leading to decreased magnitudes of seismic anisotropy. The observed textural evolution in these samples can be used to predict systematic changes in seismic anisotropy as deformation progresses (Fig. 2). Nominally melt-free samples (red squares) exhibit a systematic increase in ξ and alignment of ψ with the shear direction with increasing strain, until a steady state is reached at ξ ≈ 1.06 and ψ ≈ 0° (parallel to flow). This value of ξ is close to maximum observed values in the upper mantle, which are on the order of ξ ≈ 1.08 (13, 15), where ξ =1 indicates isotropy. Melt-present samples (4% melt fraction, blue circles) all exhibit reduced values of ξ (ξ ≈ 1.03) relative to nominally melt-free samples, because textures developed in this deformation regime are weaker. Thus, even after a melt phase has been extracted or crystallized, textures related to earlier melt-present deformation could explain the relatively low values of ξ observed at shallow depths. A similar reduction in the magnitude of azimuthal anisotropy might be expected, a feature observed in tomographic models (13, 15). In addition, melt-present samples exhibit values of ψ > 60°, typically about 90° (normal to flow). Fig. 2. Calculated values of seismic anisotropy as a function of strain. Values are given for (A) the orientation of azimuthal anisotropy (fastest direction in the horizontal plane) relative to the shear direction and (B) the magnitude of radial anisotropy taking the shear plane to be horizontal. Vertical lines indicate 1 SD. Weaker anisotropy in melt-present samples than melt-free samples results in larger error bars. Thin red squares are from Hansen et al. (22) and thin blue circles are from Qi et al. (29). All other data are from this study. Fits to data (solid lines) from melt-free aggregates and numerical mixtures are provided with 95% confidence intervals (dashed lines). A combination of laboratory observations can explain the seismological observations. In a qualitative sense, high values of ξ observed at depths of >80 km in the mantle are consistent with values in nominally melt-free samples, and low values of ξ observed at depths of <80 km are consistent with values in melt-present samples. However, fast shear-wave propagation directions observed in the upper mantle are generally much less than 60° to the absolute plate motion direction (13), a proxy for flow direction, suggesting that values of ψ both above and below 80 km are only consistent with values observed in melt-free samples. Thus, at depths of <80 km, values of ξ are consistent with melt-present samples, but values of ψ are consistent with melt-free samples. To resolve this discrepancy, we propose that regions with crystallographic textures similar to those in our melt-present samples are heterogeneously distributed throughout the mantle at depths of <80 km. If seismological surveys sample at length scales larger than these heterogeneities, the measured anisotropy would result in reduced values of ξ, but relatively unchanged values of ψ relative to the melt-free case. (See Fig. S1 for an exploration of this phenomenon.) In support of this proposition, observations of natural samples reveal heterogeneity in the distribution of melt-free and melt-present texture types throughout major portions of the lithosphere on length scales of <1 km (30), whereas tomographic models of upper-mantle anisotropy constructed from surface waves tend to have depth resolutions of >10 km. Fig. S1. Sensitivity of seismic anisotropy to fraction of textures associated with melt-present deformation. Data and models are presented as (A) the orientation of azimuthal anisotropy and (B) the magnitude of radial anisotropy, both as a function of the total shear strain. Black lines indicate the results of numerically mixing textures developed during deformation of nominally melt-free (red) and melt-present (blue) aggregates of olivine. Symbol shapes are as in Fig. 2. Individual curves are presented for different proportions of melt-present relative to melt-free textures. The solid line indicates the mixture of 80% of the anisotropy from melt-present samples and 20% of the anisotropy from melt-free samples used to predict seismic anisotropy in Figs. 3 and 4.

Laboratory Experiments This study combines data from new and previously published experiments on aggregates of olivine (Dataset S1). Similar materials and methods were used in both the old and the new experiments. Nominally melt-free samples consisted of hot-pressed aggregates of Fo 50 olivine. Fo 50 (rather than Fo 90 ) was used for its low strength relative to San Carlos olivine (Fo 90 ), which reduces the likelihood of slip at piston−sample interfaces during torsion experiments. As demonstrated by Hansen et al. (22), the olivine iron content does not significantly affect the evolution of anisotropy. Melt-present samples were hot-pressed aggregates of San Carlos olivine (Fo 90 ) and midocean ridge basalt. Melt-present samples were weak enough that it was not necessary to use Fo 50 . Textures from olivine aggregates used in this study were measured in samples deformed in torsion in a gas-medium apparatus at the University of Minnesota. Experiments were conducted at temperatures of 1,200 °C and confining pressures of 300 MPa. All samples were deformed in torsion at a constant twist rate, which corresponds to a constant strain rate for any given annulus within the sample. Shear strain rates at the maximum radius were between 2 × 10−5 and 1 × 10−2 s−1. Further details of sample preparation and experimental parameters are described in Hansen et al. (22) for nominally melt-free samples and in Qi et al. (29) for melt-present samples. We additionally examine the role of preexisting textures in subsequent textural evolution in torsion. In a midocean ridge setting, upwelling directly beneath the ridge imparts a texture before horizontal shearing as material moves away from the ridge. To simulate this effect, we initially deformed two samples in tension to a true strain of ∼80% before deformation in torsion. The observed evolution of texture and seismic anisotropy in these samples is effectively indistinguishable from that in samples deformed solely in torsion. This similarity implies that implementation of more realistic flow paths near the ridge axis in our large-scale model (e.g., corner flow) would not significantly alter the results we obtain using only simple shear (see description of model setup in Basic Flow Calculation).

Characterization of Anisotropy All olivine textures were measured from polished sections using electron backscatter diffraction (EBSD) systems at Stanford University and the University of Minnesota. EBSD data were processed by first removing wild spikes and then interpolating empty pixels with at least six nearest neighbors. All textures presented are one data point per grain, with grain boundaries defined to be misorientations of >10°. Each data set contains at least 600 grains. Seismic properties of olivine aggregates were estimated using the orientation distribution measured by EBSD and the elastic properties of olivine single crystals. The elastic moduli of olivine single crystals were taken from Abramson et al. (32). Calculation of aggregate elasticity tensors was implemented with the MTEX toolbox for MATLAB (33). The Voigt−Reuss−Hill average was implemented for all calculations (26). To calculate metrics for radial and azimuthal anisotropy, we calculated the parameters L, N, G c , and G s from aggregate elasticity tensors as described by Montagner and Nataf (34). To simulate the effect of large-scale heterogeneity in melt distribution, we combined the anisotropy of the melt-free and melt-present samples by taking the Voigt−Reuss−Hill average of their elasticity tensors. Because the elastic properties of melt-present samples do not exhibit significant evolution with increasing strain, elasticity tensors of each melt-free sample were averaged with the elasticity tensor of a typical melt-present sample (PT0767). The error in each calculation was estimated with the bootstrapping technique used by Hansen et al. (22). Each orientation distribution was undersampled to 100 orientations, and the azimuthal and radial anisotropies were calculated from this subset. This procedure was repeated 100 times to build a distribution of estimates of the magnitude and orientation of anisotropy, the SDs of which were used as an indicator of error in the measurement. Our dataset examines azimuthal anisotropy at low strains in torsion and demonstrates that azimuthal anisotropy in nominally melt-free samples is at a high angle to the flow direction at low strains. This observation results from the dominant [100] direction being >45° to the shear plane at low strain. The orientation of azimuthal anisotropy at low strains is counterintuitive in the context of previous results, which suggest the fast Raleigh wave direction will always be parallel to the shear direction (24). However, previous estimates of azimuthal anisotropy at low strain have either come from experiments conducted in transpression (24, 25) or simulations (35, 36) calibrated using these experiments, despite models that demonstrate significant differences in textural evolution between transpression and simple shear (37). Our observation of anisotropy at low strains in simple shear demonstrates this difference, due to vertically polarized Rayleigh waves becoming the fastest waves when [100] is at a high angle to the shear plane. We weighted the properties of both elasticity tensors to simulate a material with 80% of the volume characterized by deformation of melt-present samples and 20% of the volume characterized by deformation of nominally melt-free samples. This ratio yielded an evolution of azimuthal anisotropy similar to the melt-free case and an evolution of radial anisotropy similar to the melt-present case. If the volume characterized by melt-present deformation is reduced to 50%, maximum value of ξ increases to 1.04%, and if the volume is increased to 90%, the strain to reach a particular value of ψ increases by a factor of ∼5. Neither case significantly changes our results. The sensitivity of the calculated seismic characteristics to different mixing proportions is presented in Fig. S1. We parameterized the evolution of azimuthal anisotropy in simple shear using an inverse exponential. A least-squares fit to the data from nominally melt-free samples yielded ψ = ( 80.8 ± 7.7 ) exp [ − ( 0.75 ± 0.12 ) γ ] , [S1]where γ is the shear strain in simple shear; confidence intervals are at the 95% level. A least-squares fit to the data from the numerical mixture of melt-free and melt-present samples yielded ψ = ( 88.6 ± 6.1 ) exp [ − ( 0.41 ± 0.05 ) γ ] . [S2]We parameterized the evolution of radial anisotropy in simple shear using a hyperbolic tangent function. A least-squares fit to the data from nominally melt-free samples yielded ξ = ( 0.079 ± 0.005 ) tanh [ ( 0.28 ± 0.05 ) γ ] + ( 0.982 ± 0.005 ) . [S3]A least squares fit to the data from melt-present samples yielded ξ = ( 0.017 ± 0.001 ) tanh [ ( 0.30 ± 0.04 ) γ ] + ( 1.015 ± 0.001 ) . [S4]

Basic Flow Calculation To simulate flow in the upper mantle associated with the horizontal translation of a plate, we set up a simple 1D model for simple shear (Couette flow) with constant velocity boundary conditions similar to those used by Podolefsky et al. (38) and Behn et al. (39) (Fig. S2). We took into account the upper mantle from 0 to 400 km depth, assuming a horizontal velocity of 100 mm/y maintained at the upper boundary and a velocity of 0 mm/y maintained at the lower boundary. The shear stress was considered to be spatially constant but able to vary temporally. The mechanical properties were governed by a flow law of the form ε ˙ = A σ n ⁡ exp ( − Q + P V R T ) , [S5]where ε ˙ is the strain rate, A is a material constant, σ is the applied stress, n is the stress exponent, Q is the activation energy, P is the confining pressure, V is the activation volume, R is the gas constant, and T is the temperature. Values for A, n, and Q were taken from Hirth and Kohlstedt (40) for dry dislocation creep of olivine. The flow law above is formulated for triaxial compression. To facilitate application of this flow law to stresses and strain rates in simple shear, engineering shear strain rates were calculated as γ ˙ = 3 ε ˙ , and shear stresses were calculated as τ = σ / 3 (e.g., ref. 41). A value of 14 × 10−6 m3/mol was used for the activation volume (42). Although the upper mantle contains some dissociated water incorporated into the olivine structure, we do not implement a wet flow law. In the context of the evolution of seismic anisotropy, the effect of water on viscosity is an unnecessary complication. Using the wet flow law (40) with 300 ppm H/Si in olivine reduces the applied stress (and therefore the viscosity) by a factor of ∼2 (e.g., ref. 43), but the magnitude and distribution of strain are not affected. Thus, the viscosity color scale in Fig. 4 and Fig. S2 shifts slightly when water is present, but the plots are otherwise unchanged. This difference is within the experimental error associated with the extrapolation to geological rates. For the mantle melting calculation presented in Calculation of the Spatial Distribution of Anisotropy we do include 300 ppm H/Si in olivine, as this significantly influences the depth and extent of melting. Although several deformation mechanisms are likely operating simultaneously in the upper mantle (40), we suggest a reasonable first-order approximation of the flow pattern is obtained by only considering dislocation creep. Diffusion creep is potentially contributing significantly at the base of the upper mantle, but predictions of its depth extent depend strongly on the poorly constrained values of the activation volume for each creep mechanism. Recent estimates suggest that diffusion creep does not contribute significantly to the deformation over the upper 200 km of the mantle and potentially the upper 400 km (43). Because we are primarily concerned with structure in the uppermost mantle, we opt not to include diffusion creep. Recent work has also established that a grain-size-sensitive mechanism such as dislocation-accommodated grain-boundary sliding (disGBS) may be dominant throughout large parts of melt-free (44, 45) and melt-present (46) upper mantle. Current understanding of this mechanism suggests that it is largely similar to dislocation creep in terms of its microphysical processes (44, 45). Furthermore, the textures developed in olivine by dislocation creep and disGBS are virtually indistinguishable (22, 47). Thus, the primary potential influence of disGBS in our model is the addition of a grain size dependence of the viscosity. To fully incorporate grain-size-dependent flow would also require prediction of the grain size evolution. As we are interested in a first-order description of the flow field, we do not include disGBS, although it can have important implications for the spatial distribution of grain sizes and, therefore, seismic properties such as attenuation (39). Preliminary tests of a model that includes grain size evolution along with flow laws for disGBS and diffusion creep did not demonstrably change the results presented here. Finally, we note that our model does not incorporate two pressure-related transitions in deformation behavior that may become important in the deep uppermost mantle. First, transitions in the activation volume at greater depths will potentially modify the strain rate (and therefore strain) distribution below 200 to 300 km depth (48). Although modifying the activation volume will change the magnitude of strain in the deep upper mantle, it will not affect the discontinuities in the uppermost (<200 km) mantle considered here. Second, there is increasing experimental evidence for a transition in textural evolution at pressures corresponding to depths below ∼200 km (48⇓–50). Although this transition in textural evolution would likely lead to a change in the orientation of azimuthal anisotropy, this effect is only relevant at depths below those considered in the current study. The time dependence of flow in the upper mantle was investigated by allowing the temperature to change according to a half-space cooling model (51) (Fig. S2). Material properties incorporated into this model include a mantle potential temperature of 1,350 °C, a specific heat of 1,350 J⋅kg−1⋅K−1, a coefficient of thermal expansion of 2.9 × 10−5 kg−1, and a thermal diffusivity of 10−6 m2/s. The stress at any given time was found through an iterative procedure. A stress was guessed and a strain rate profile was calculated using the current temperature profile. The velocity at the upper boundary was then calculated from the resulting strain rate profile. This procedure was iterated over a range of stresses until the stress yielding an upper velocity closest to 100 mm/y was found. Because the mean temperature of the system decreased throughout the simulation, the applied stress increased from 0.8 to 7.1 MPa over the course of 200 My. The finite strain at any depth and time was calculated using a central finite differencing scheme (Fig. S2).

Calculation of the Spatial Distribution of Anisotropy We used the results presented in Fig. 2 to predict seismic properties over the depth and age intervals included in our flow model. At depths below those at which melting occurs at any time, the predicted finite strain was used in conjunction with Eq. S3 to map the distribution or radial anisotropy. Fig. S3 presents a comparison between the radial anisotropy developed in our flow model at 200 km and that observed in the tomographic model of Beghein et al. (13). The good correspondence in timescale of the evolution of anisotropy between model and observation increases confidence in extrapolation of laboratory results to upper-mantle processes. We assumed melt is present at depths and ages at which melting is predicted to occur according to the peridotite solidus (yellow region in Fig. S2, Lower Right). We modified the solidus of Hirschmann (52) for a water content of 300 H/106 Si in olivine following the argument of Hirth and Kohlstedt (8). The evolution of radial anisotropy in partially molten regions was calculated with Eq. S4. We do not modify the flow law in the partially molten region, even though the presence of a partial melt reduces viscosity. A viscosity reduction due to the presence of melt has little to no effect on the evolution of anisotropy, because the magnitude and orientation of anisotropy in melt-present samples is not a function of strain (Fig. 2). After material has been transported away from the ridge and temperature cooled below the solidus, we begin to describe the evolution of radial anisotropy using Eq. S3. A difficulty arises in that Eq. S3 assumes that the entire strain history has been in the absence of melt. We make the first-order approximation that, once the system cools below the solidus, the magnitude of radial anisotropy will increase at a rate expected for a melt-free aggregate with similar anisotropy. In other words, the anisotropy will increase according to the curve for the melt-free case (red lines in Fig. 2), beginning at a magnitude of anisotropy similar to that at the end of the melt-rich deformation. As an example of the resulting output incorporating textures formed during melt-free and melt-present deformation, a 1D profile of radial anisotropy at a plate age of 100 My and its vertical derivative are provided in Fig. S4. These profiles are comparable to those presented by Beghein et al. (13) in their figure S9, with the main difference being the magnitude of the vertical gradient near the discontinuity at ∼80 km. We suggest that the smaller gradient in Beghein et al. (13) is due to the smoothed natured of the tomographic model. Fig. S4. Predicted distribution of seismic anisotropy in the Pacific upper mantle for a plate age of 100 My. Blue curves indicate the model output, and red curves indicate the smoothed output. Right panels illustrate the vertical derivatives of the corresponding Left panels. The gray dashed line indicates the location of the anisotropy discontinuity as determined by the (Upper) maximum gradient in radial anisotropy and (Lower) minimum gradient in the fast azimuth deviation. We predicted the distribution of orientations of azimuthal anisotropy in a similar manner to radial anisotropy. Because our flow model is only 1D (or pseudo-2D when vertical profiles at neighboring times are juxtaposed), we cannot intrinsically predict changes in the fast azimuth for wave propagation. Stratification of azimuthal anisotropy likely results from anisotropy being frozen in as the lithosphere thickens, followed by subsequent changes in plate motion direction. As an ad hoc solution, we suggest that the plate motion direction changes azimuth by 0.3°/My, a value consistent with the angular velocity of the Euler pole for the Pacific plate (53). We thus calculate the azimuthal anisotropy according to Eq. S1, but we modify the apparent strain to account for the misalignment between the plate motion direction and the fast azimuth that has occurred due to plate rotation in that time step. Regions deforming at faster strain rates allow for increased alignment of the predicted fast azimuth with the plate motion direction. We do not modify this scheme for regions that have undergone melting because, as illustrated in Fig. 2A, the timescale for realignment of azimuthal anisotropy is similar for both nominally melt-free samples (red squares) and our numerical mixture of textures from melt-free and melt-present samples (black triangles). We use the difference between the predicted fast azimuth and the plate motion direction to assess which regions are best aligned with the local plate motion direction. An example 1D profile of the resulting prediction of azimuthal anisotropy at a plate age of 100 My and its vertical derivative are provided in Fig. S4. As noted above for gradients in radial anisotropy, the magnitude of the vertical derivative of the deviation in the fast azimuth direction is larger in our study than presented by Beghein et al. (13), which we attribute to the smoothed nature of the tomographic model.

Acknowledgments The authors thank Thorsten Becker, Caroline Beghein, Jean-Paul Montagner, and Luiz Morales. We thank David Kohlstedt for input on several versions of this manuscript. This research was supported by funding awards, including John Fell Fund 123/718 and Natural Environment Research Council (NERC) NE/M000966/1 (to L.N.H.), National Science Foundation (NSF) EAR-1255620 (to J.M.W.), and NSF EAR-1214876 (to D. L. Kohlstedt).

Footnotes Author contributions: L.N.H. designed research; L.N.H. and C.Q. performed research; L.N.H. analyzed data; and L.N.H., C.Q., and J.M.W. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1608269113/-/DCSupplemental.