Research ArticleGEOPHYSICS

Seismic anisotropy reveals crustal flow driven by mantle vertical loading in the Pacific NW

See allHide authors and affiliations

Science Advances  08 Jul 2020:
Vol. 6, no. 28, eabb0476
DOI: 10.1126/sciadv.abb0476


Buoyancy anomalies within Earth’s mantle create large convective currents that are thought to control the evolution of the lithosphere. While tectonic plate motions provide evidence for this relation, the mechanism by which mantle processes influence near-surface tectonics remains elusive. Here, we present an azimuthal anisotropy model for the Pacific Northwest crust that strongly correlates with high-velocity structures in the underlying mantle but shows no association with the regional mantle flow field. We suggest that the crustal anisotropy is decoupled from horizontal basal tractions and, instead, created by upper mantle vertical loading, which generates pressure gradients that drive channelized flow in the mid-lower crust. We then demonstrate the interplay between mantle heterogeneities and lithosphere dynamics by predicting the viscous crustal flow that is driven by local buoyancy sources within the upper mantle. Our findings reveal how mantle vertical load distribution can actively control crustal deformation on a scale of several hundred kilometers.


Geodynamic models commonly describe the relation between buoyancy-driven mantle convection and plate tectonics with two components of traction applied to the base of the lithosphere—vertical tractions giving rise to dynamic topography (1, 2), and horizontal basal tractions driving plate motion and tectonic deformation (3). Although lithospheric stress field measurements (4, 5) and mantle flow patterns (6, 7) provide critical constraints on the dynamics of these interactions, attempts to isolate their relative influence on near-surface tectonics often yield ambiguous results. In most cases, the difficulty derives from our imperfect knowledge of the mantle density structure and the high variability in the material strength of the lithospheric rocks, which greatly influences the degree of mechanical coupling between the tectonic plates and the convective flow (8). While substantial advancements in seismic imaging have permitted the construction of high-resolution models of the mantle’s mass distribution, an ability to accurately quantify the degree of coupling between the mantle and the lithosphere remains underdeveloped. This limitation, in combination with a paucity of observational constraints, has prevented any reliable assessment of how mantle-based forces interact with plate-scale processes to give rise to the tectonic stresses that drive surface deformation. Here, we show that, under certain rheological conditions, crustal anisotropy is transparent to the structural complications of the crust and can reveal a crustal flow driven by the vertical coupling of the mantle and the lithosphere.

The process for detecting mantle-induced vertical deformation generally involves identifying regions that have experienced rapid surface uplift or subsidence (9, 10) or areas with sharp elevation contrasts that are difficult to explain with simple isostatic models (11). The Wallowa Mountain block in northeastern (NE) Oregon, for example, represents a remarkable regime where both mantle- and crustal-based stresses appear to have played an essential role in lithosphere dynamics. These mountains are composed of a sizable granitic batholith that rapidly rose ~2 km above the surrounding area shortly after the deposition of the Columbia River flood basalts (CRB) ~16 million years (Ma) ago (12, 13), creating an impressive topographic bullseye centered on the Wallowa Mountains (Fig. 1). The compact and isolated uplift of the granitic Wallowa batholith suggests the foundering (13) of a dense garnet-rich (14) pluton root during or shortly after the CRB eruptions. However, the regional post-CRB uplift of the entire topographic bullseye region (14) indicates the existence of a larger-scale mechanism that dynamically drives crustal deformation around the site of the inferred foundering event. In the mantle beneath the Wallowa Mountains, high-resolution tomographic images persistently reveal the presence of a major high-velocity anomaly (the Wallowa anomaly) that is circular in map view and extends to a depth of ~350 km (15, 16). This structure appears to be part of a system of ancient slab fragments that are dangling beneath the Pacific Northwest (NW) and, together with buoyant plumes of rising asthenosphere, is hypothesized to drive the small-scale mantle convection that actively modifies the western U.S. lithosphere (16). Although the precise role of the Wallowa anomaly and other nearby mantle heterogeneities in shaping the topography of NE Oregon remains unclear, recent seismic imaging studies have revealed that the crust just north of the Wallowas is about 20 km thicker than the surrounding area (17, 18). The correlation between these two puzzling features seems to suggest that the negatively buoyant Wallowa anomaly is responsible for the localized pull-down on the Moho.

Fig. 1 Regional and location maps for northeastern (NE) Oregon.

(A) Regional map showing the broadband seismic stations used in this study (black inverted triangles). The dashed blue line depicts the Snake River Plain (SRP). WA, Washington; OR, Oregon; ID, Idaho. (B) Global map centered in NE Oregon. The thick black line encloses the region shown in (A). (C) Elevation map of the topographic bullseye region [red area in (A)]. The dashed larger ellipse is the outer limit of the bullseye, whereas the inner ellipse locates the Wallowa batholith.

In this study, we developed an azimuthal anisotropy model for the crust of NE Oregon and its surrounding regions using short-period (3 to 17 s) Rayleigh waves extracted from ambient seismic noise cross-correlations (Fig. 2). With the concentration of broadband stations in the area and wide azimuthal interstation path coverage, we can reliably resolve the lateral variations of seismic anisotropy for the uppermost ~35 km of the crust. To measure the anisotropy, we implemented a beamforming scheme that allows us to characterize the seismic wavefield’s velocity dependence with propagation direction beneath each station (19, 20). The reliability of our model was then verified by comparing our surface wave anisotropy measurements to those that were obtained by characterizing the azimuthal dependence of receiver functions (RFs) at stations surrounding the Wallowa Mountains (fig. S1).

Fig. 2 Example of beamformer outputs and final azimuthal anisotropy model.

(A) Two-dimensional histograms over the azimuth-velocity space for the 3- to 17-s period band with the best-fitting anisotropy model (red dashed lines). The green bars on top of each panel indicate the number of noise cross-correlations available for each azimuth. (B) Azimuthal anisotropy model for the crust of the Pacific NW. Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color represents the intersection density of the anisotropy vectors assuming that they are of infinite length (i.e., projected to the bounds of the study region). The green and blue dots indicate the location of the two stations beamformed in (A).


The azimuthal anisotropy model for the crust underlying this region does not correlate with the surface geology, the structural trends, or the mapped crustal stress field (21). However, it instead holds a remarkable connection to the upper mantle velocity distribution (Fig. 3A). The fast directions of anisotropy, which are generally thought to reflect the coherent deformation of small-scale structures and preferred alignment of anisotropic minerals (22), show a simple and well-defined radial pattern that strongly correlates with the Wallowa anomaly. Moreover, the northern and easternmost anisotropy vectors display a subtle fan-like pattern that correlates notably well with the geographic extent of the Siletzia slab curtain beneath Idaho (15, 16). The amplitude of the azimuthal anisotropy also decreases to near-zero values for seismic stations above the Wallowa anomaly and slab curtain kink, where the geometry of anisotropy may be transitioning into one that is null in the horizontal plane. The connection between crustal anisotropy and upper mantle velocity structure suggests that mantle gravitational loads actively induce vertical stresses on the overlying material and, through this relation, control crustal deformation in NE Oregon and its adjacent regions.

Fig. 3 Comparison of seismic and geodynamic results with crustal anisotropy.

(A) Azimuthal anisotropy model for the crust of the Pacific NW overlying a depth slice through the Vp tomography model at 250 km (16). The red dashed lines depict the Wallowa anomaly and the Siletzia slab curtain. (B) Modeled Moho stress and mid-lower crustal flow velocity for the Pacific NW. The colored contours represent the vertical stress at the Moho based on a global geodynamic model driven by density anomalies derived from the P-wave velocity structure. The black arrows denote the predicted mid-lower crustal flow velocity that results from the application of the modeled Moho stress to a viscously heterogeneous crust. The red bars represent the anisotropy measurements derived from this study.

On the basis of the spatial coherence of our measurements, we hypothesize that the azimuthal anisotropy in NE Oregon results from the lattice preferred orientation (LPO) of anisotropic minerals with the subhorizontal flow of the mid-lower crust (23, 24). Recent numerical studies show that stresses transmitted upward from the underlying mantle can induce large amounts of intraplate deformation through Poiseuille and Couette flow due to lateral pressure variations and basal shear (25). This style of deformation requires the lithospheric rocks to have low viscous strength to form a channeled ductile flow system in the mid-lower crust that decouples the upper crustal and upper mantle stress fields. Because of the relatively recent magmatic activity in NE Oregon, the crust beneath this region can achieve the adequate thermal conditions (700° to 1000°C) to create such a ductile and mobile environment, especially along the Snake River Plain and beneath the Wallowa Mountains (26, 27). The existence of a mid-lower crustal weak channel would then allow the crust to flow in response to lateral pressure gradients and accommodate the vertical stresses exerted by the underlying mantle. Because of the high local Moho temperature, the mapped crustal anisotropy is most likely to be dominated by the LPO of type II and III fabrics in amphibole, for which the fast direction of anisotropy is subparallel to the flow direction (28). The alignment of mica crystals may also contribute to the overall observed anisotropy; however, global compilations of the structure of the continental crust suggest that its deep portion contains rather little mica and that amphibole takes a larger fraction of its composition (29).

At sublithospheric depths, shear wave splitting observations reveal that the asthenospheric flow in the Pacific NW is primarily controlled by a combination of North American plate motion and the sinking of the Juan de Fuca and Farallon slab systems (30, 31). These measurements also reveal that there is little, if any, mantle deformation caused by the downward movement of the Wallowa anomaly. As a matter of fact, the most recent shear wave splitting observations in NE Oregon indicate that mantle materials flow smoothly around the lateral boundaries of the Wallowa anomaly rather than converging on the site of lithospheric load (Fig. 4) (32). The lack of a strong disturbance in the mantle flow field beneath this area suggests that the asthenospheric strain that is created by the downwelling velocity of the Wallowa anomaly is not strong enough to perturb the current LPO that has been established by the long-term movement of the tectonic plates. This observation leads us to the notion that whatever vertical forces are being exerted by the upper mantle and driving the crustal flow are almost entirely derived from the negative buoyancy of its dense structures rather than the weak vertical asthenosphere flow that is excited by their vertical movement. Furthermore, the absence of correlation between the crustal and the upper mantle deformation fields suggests that the mantle lithospheric strength isolates horizontal asthenospheric flow from that in the crust, such that there is insignificant basal shearing by the underlying mantle (i.e., mantle flow is not driving Couette flow in the viscous crust). This type of decoupling is consistent with previous tomographic findings that show a weak correlation between crustal and upper mantle anisotropy in most regions of the western United States (33). Here, it is important to note that a key to the low strain rate of the dense mantle structures is their greater viscosity and that their low sinking rate is a result of them being attached to the North American lithosphere, as seismically imaged (15, 16).

Fig. 4 Station averaged shear wave splitting measurements for the Pacific NW.

(A) SKS splitting measurements for the entire western United States (31). The red arrow depicts the relative motion between the North American plate (NA) and the hot spot reference frame (HS) (55). (B) SKS splitting measurements for our study region [red area in (A)]. The thick blue vectors depict the measurements of Niday and Humphreys (32), and the black vectors are from the database of Becker et al. (31). In both panels, the orientation of the vectors gives the angle of the fast polarization, and the length of the bars is proportional to the magnitude of the shear wave splitting. The white trajectories through the anisotropy field lines in (B) are used to represent the streamlines of the mantle flow assuming an east-oriented flow (56). The red circle in the background marks the location of the Wallowa anomaly. Note how mantle materials appear to flow smoothly around the lateral boundaries of the Wallowa anomaly.

In the ductile regime, viscous strain rate preferentially orients minerals relative to one another, generating seismic anisotropy that is aligned with the flow direction (34). Within this framework, we model the crustal deformation induced by mantle loading through scaling the seismically imaged mantle velocity anomalies (16) to density structure (35) and predicting vertical tractions at Moho depths. Here, we exclude the subducting Juan de Fuca slab and North American craton because these structures have been in long-term steady state relative to the more recent mantle structures of the interior Pacific NW and are, therefore, unlikely to play a crucial role in shaping the present-day crustal strain field. We also impose a small-scale load, with moderate stress magnitude, to the predicted Moho traction in the Wallowas area to account for the foundering of the mountain’s pluton root. This last step is taken because this event is a rather short-lived transient phenomenon that is not captured by the seismic tomography and may also contribute to the observed anisotropy; the destabilization and subsequent removal of the root would cause the weak mid-lower crust around the Wallowas to flow toward the vacated root region (25). Moho traction calculations made without incorporating the localized Wallowa load reveal that the first-order vertical stress distribution is not significantly altered and are presented in fig. S2.

The final mantle-derived vertical tractions are then used to drive viscous Stokes flow in a thin crustal channel using surface heat flow as a modulator of crustal viscosity (36). Figure 3B shows a comparison between the observed crustal anisotropy, the estimated vertical stress at the Moho, and the predicted mid-lower crustal flow. Both the relative amplitude and orientation of the crustal flow velocity agree remarkably well with the measured crustal anisotropy within the main study area, displaying a dominant radial flow pattern centralized at the Wallowa Mountains site. Such flow would lead to crustal thickening in the low-pressure regions, which is consistent with the nearly circular ~20-km Moho depression that is observed above the Wallowa mantle anomaly (17, 18). Note, however, that there exists a disagreement between our modeled and measured anisotropy in the southwestern part of our study region. This discrepancy may be the result of either nonmantle processes that are not included in our numerical model (i.e., tectonic strain in the active Basin and Range) or due to the inherent shortcoming of the beamforming technique to resolve lateral sharp changes in the anisotropic structure. Nonetheless, the general agreement between the crustal flow predicted by our simple model and the seismic observables strongly suggests that mantle-induced stresses can, in some cases, have more control on intraplate deformation than those transmitted laterally from nearby active plate boundaries (2).


On the basis that the lithosphere is rheologically stratified, we propose an upside-down water bed model in which vertical mantle loads cause the ductile rocks inside the weak mechanical layer to migrate horizontally toward low-pressure regions through Poiseuille flow, involving little mantle deformation (Fig. 5). The mechanics of the channelized flow that is induced by this model are similar to the ones that are typically invoked to explain near-surface deformation in extreme tectonic environments such as the Tibetan Plateau (37) or the Altiplano in the Bolivian Andes (38). The difference, however, lies on the fact that crustal flow in such regimes is generally thought to be driven tectonically or gravitationally as a response of the buoyancy forces that arise from differential crustal densities (39) or the pressure differences caused by varying crustal thickness (40). Therefore, the evidence that mantle gravitational loads are capable of displacing weak crustal materials in a comparable manner not only refines our understanding of the interaction between crustal tectonics and mantle dynamics but also brings to light another source of deformation that might be necessary to explain the state of stress in the crust of other tectonically enigmatic regions. Figure 6 for instance shows the crustal anisotropy measurements of Lin et al. (33) around Southern and Central California, where other dense mantle anomalies have been imaged by different tomographic studies. Similar to the case of the Wallowas, the strong correlation between the crustal anisotropy and upper mantle velocity structure in this region suggests that the upside-down water bed model is playing a crucial role in driving the evolution of the lithosphere. However, because of the unique tectonic of California, deciphering the precise role and contribution of mantle-based stresses in its surficial processes would require a more complete modeling that incorporates the rather elevated horizontal strains that are exerted from the active plate margin.

Fig. 5 Schematic representation of the upside-down water bed model.

The load of the mantle lithosphere is a force creating vertical stresses (σzz) on the Moho. The lithospheric load pulls down on the crust, which creates a lateral pressure gradient that drives Poiseuille flow in the ductile mid-lower crust. The asthenosphere flows independently (as evidenced by its independent anisotropy field; Fig. 4), creating a local Couette flow that is decoupled from the mid-lower crust by the mantle lithosphere.

Fig. 6 Crustal anisotropy and upper mantle velocity structure of California.

The black vectors depict Lin et al.’s (33) surface wave anisotropy measurements for the 12-s period. Bar orientation gives the fast direction of azimuthal anisotropy, and bar length is proportional to anisotropy amplitude. The background color corresponds to a depth slice through the Vp tomography model at 195 km (16). The red dashed lines denote the two seismically fast and likely dense mantle anomalies.

In general, the crustal anisotropy that results from the upside-down water bed model is most sensitive to recent deformation and, hence, relevant for addressing young tectonic evolution. For the case of NE Oregon, the temporal sequence of the Siletzia slab curtain formation ~53 Ma ago and the Wallowa anomaly delamination ~16 Ma ago (15) may well explain the apparent dominance of the Wallowa anomaly in aligning the fast directions of anisotropy. This argument is supported by the fact that, although the additional load imposed at the site of the Wallowa batholith was initially designed to incorporate the effects of the foundering of its root, the localized Wallowa enhancement would still be required to represent the latest stage of upper mantle vertical forcing and achieve the remarkable centralized flow pattern that is illuminated by the crustal anisotropy. The flow that is predicted by mantle loading alone thus suggests that the amount of strain exerted by the delamination of the Wallowa mantle anomaly is enough to effectively align mid-crustal minerals and overprint on any prestrained fabric. An alternative and simple interpretation is that the viscous strength of the crustal rocks beneath and around the Wallowa Mountains is weaker because of the recent CRB eruptions (27) and have deformed easily, flowing toward the site of the mantle loading and observed crustal thickening. Regardless of the relative level of contribution of each of these mechanisms, our findings provide strong observational evidence of regional-scale mantle-crust vertical coupling and highlight the fundamental importance of upper mantle buoyancy in understanding near-surface tectonics.


Ambient seismic noise cross-correlation

The dataset used in this study consists of Rayleigh waves obtained from the cross-correlation of ambient noise recorded at over 350 broadband seismic stations (table S1). To extract the surface wave signals, we used a similar routine to the one described in Bensen et al. (41) and computed the vertical-to-vertical cross-correlations of continuous recordings between all synchronous station pairs. The preprocessing of the ambient noise data consisted in (i) downsampling the continuous records to one sample per second and dividing them into 1-day time windows, (ii) removing the mean and trend in each time window, (iii) band-pass filtering between the 3- and 100-s period band, (iv) whitening the spectra, and (v) normalizing in the time domain. After preprocessing, each time window was cross-correlated, normalized to unit peak amplitude, and averaged over time. To further enhance the surface wave signals, the causal and anticausal parts of the cross-correlation functions were stacked to obtain the so-called symmetric cross-correlations. The application of this technique yielded over 50,000 ambient noise cross-correlation functions with a signal-to-noise ratio (SNR) greater than 10. Here, we defined the SNR by the ratio of the peak amplitude of the surface wave signal to the root mean square of the noise trailing the ballistic arrival (41).

Azimuthal anisotropy measurements

We measured the azimuthal anisotropy at individual stations by finding the average group velocity of fundamental-mode Rayleigh waves traveling from all available backazimuth ranges and characterizing the wavefield’s azimuthal dependence. We choose this method over traditional tomographic inversions because the station distribution in NE Oregon provides us with sufficiently dense coverage to estimate the anisotropy at individual stations and still be able to interpret the regional crustal flow field as a whole. This advantage allowed us to directly quantify the variations of wave parameters with the propagation direction at discrete locations without the necessity of imposing any type of regularization, which tends to bias the final solution.

To characterize the wavefield’s azimuthal dependence at a given station, we adopted a beamform approach and used a plane wave approximation to find the best-fitting slowness for every possible azimuth range (19, 20). For this process, we searched for the maximum coherent output over velocities from 1 to 5 km/s in 5° bins from 0° to 360° azimuth with 70% overlap for the 3- to 17-s period band. We used this frequency range so that Rayleigh waves were exclusively sensitive to Earth’s crust and not to upper mantle properties (fig. S3). To ensure the robustness of our measurements, we only beamformed cross-correlations with an SNR higher than 10 and an interstation distance larger than one wavelength of the lowest period of the band-pass filter (41). We also only used stations at which the azimuth range of 180° was sampled by at least three ray paths in a five-bin range.

To measure the anisotropy, we used a minimization algorithm and fitted the first three parameters of Smith and Dahlen’s (42) anisotropy model for surface wave velocityv(T,θ)=a0(T)+a1(T)cos(2θ)+a2(T)sin(2θ)+a3(T)cos(4θ)+a4(T)sin(4θ)(1)where v is the velocity, T is the period, θ is the azimuth, a0 is the isotropic velocity, and a1−4 are the azimuthal coefficients (43) to every beamformer output. The reason why we only kept up to the 2θ terms is that the azimuthal dependence of the fundamental-mode Rayleigh waves is practically insensitive to the 4θ terms (44). After obtaining the azimuthal coefficients, we calculated the magnitude of the anisotropy, A, and its seismically fast direction, ∅, usingA=a12+a22(2)=12tan1a2a1(3)

To assess the uncertainty in our anisotropy parameters, as well as their statistical significance, we applied a bootstrapping method to calculate the 95% confidence limits using a total of 100 resamples. Figure S4 presents the bootstrap SE for the A and ∅ parameters for all stations.

Azimuthal anisotropy from RFs

To test the reliability of our surface wave measurements, we used RFs to measure the time variations of P-to-S converted phases at the Moho (Pms) beneath 16 seismic stations surrounding the Wallowa Mountains (4 stations for each geographic quadrant) and evaluated their similarity to the surface wave results. To this end, we made use of the EarthScope Automated Receiver Survey (45, 46) to obtain the RFs and quantified the crustal anisotropy at each site following the workflow of Shen et al. (47).

To measure the anisotropy, we first applied a moveout correction for a reference slowness of 6.4 S/° to all available RFs. Then, for a given station, we stacked the radial RFs in a 10° azimuth bin and characterized the systematic moveout of the Moho Pms converted phases using a model that assumes a single layer of anisotropy with a horizontal symmetry axistpms=t0δt2cos[2(ψc)](4)and a model that assumes a single layer of anisotropy with a tilted symmetry axis of anisotropy (48)tpms=t0+δt2cos(ψc)(5)

In both models, t0 represents the reference time of arrival of the Pms phase, δt the split time caused by the anisotropy, ψ the backazimuth of the RF, and ∅c the fast direction of anisotropy. For each seismic station, we conserved the anisotropy model that resulted in the largest stacking energy. In the few isolated cases where both models resulted in similarly good fit, we performed an azimuth weighted stack (AWST) to distinguish the periodicity of the converted signals in the tangential RFs. To build the AWST sections, we used the simplified functional form of Girardin and Farra (49)ST(k,t,ψ)=i=1nsink(ψi)Ti(t)(6)where Ti(t) is the individual tangential RF after moveout correction, ∅i is the backazimuth of the ith RF, n is the number of RFs, k is the harmonic order, and ψ is a variable parameter of azimuth. It can be shown that, in the presence of anisotropy, the amplitude of the AWST section will have a maximum amplitude around the Pms phase whenever ψ = ∅c. Here, we limited the harmonic order to k = 1 and k = 2, as the AWST is mostly sensitive to the anisotropy and dip of the symmetry axis for k = 2. We also only used seismic stations that had at least nine binned RFs and no 180° azimuthal gap. Figures S5 to S8 show the anisotropy characterization process of the 16 seismic stations.

Geodynamic modeling of crustal flow

We numerically simulate the water bed model in two steps. First, we calculate the expected radial normal stress at the Moho resulting from mantle density anomalies. We apply that stress field as a forcing term to the lower crust, modeled as channelized fluid flow in a heterogeneous viscous layer. We interpret the instantaneous strain rate in the lower crustal channel as an indicator of time-integrated finite strain and, therefore, as an indirect proxy for expected LPO and seismic anisotropy.

To calculate the expected Moho stress state, we apply a simple geodynamic model using the HC global mantle convection code (50). We use the tomographic Vp (P-wave velocity) image of Schmandt and Humphreys (16) as a prescribed body force, with depth-dependent density scaling against the PREM (Preliminary Reference Earth Model) (fig. S9) (51). We impose neutral density everywhere outside of the study area, and the structures within the study area are shifted to produce net-neutral buoyancy conditions within our model space. We use the reference mantle viscosity profile of Steinberger and Calderwood (52), as implemented by Steinberger and Holme (53). Variations from this reference model, including adjusting the thickness and viscosity of the lithosphere, have minimal effect on the qualitative results. From these geodynamic calculations, we extract radial normal stress predicted at the uppermost model layer. Last, to incorporate the Wallowa enhancement, we add a 55-km-radius Gaussian-shaped load with moderate stress magnitude to the predicted Moho traction in the Wallowas area. This composed stress field is then applied as a forcing term for our crustal deformation model.

We consider the lower crust as a viscous fluid undergoing parabolic Poiseuille flow in a heterogeneous narrow channel. Neglecting internal body forces, our system is governed by Stokes momentum balance·ηε̇¯(u)P=0(7)with velocity, u; viscosity, η; strain rate, ε̇¯=(u+(u)T)/2; and pressure, P. Boldface symbols represent vector quantities, and underlines indicate tensors of rank 2.

The lower crustal channel is modeled in a layer-centered coordinate system, where z = 0 at the channel’s midplane. The channel has half thickness, h, and lateral components of flow are assumed to have parabolic form, tapering to zero at the upper and lower channel boundariesux,y(z)=U(1z2h2)(8)where U is a two-dimensional vector field representing the flow velocity at z = 0. By isolating the channel’s midplane, combining Eqs. (7 and (8 yieldsη2U2ηh2UP=0(9)

Substituting the Moho stress field (calculated above) for P, the system becomes well posed and can be solved numerically.

To obtain numerical solutions to Eq. (9, we first nondimensionalize the system. We define nondimensional variables, indicated by a “′” superscript, in terms of representative magnitudes, indicated by a “0” subscript: η = η0 η′, Ui = U0 Ui′, P = P0 P′, h = h0 h′, xi = h0 x′i, and ∇ = h0 – 1 ∇′, and rewrite the system in terms of a nondimensional parameterΠ=U0η0h0P0(10)

Equation (9 then becomesΠη(2U2h2U)P=0(11)

Finite element implementation

We use the TerraFERMA modeling framework (54) to solve Eq. (11. We scale all empirically derived variables by representative magnitudes: h0 = 20 km, η0 = 1 × 1020 Pa/s, u0 = 3.17 × 10−9 m/s (10 cm/year), and P0 = 1 MPa, equivalent to Π ≈ 15. We assign h ′ = 1 everywhere. Varying crustal thickness in proportion to Moho topography has negligible effect on the form or magnitude of the solution. Viscosity perturbation (η′) varies exponentially with observed surface heat flow of Blackwell et al. (36). Applicable coefficients are scaled such that 0.1 < η′ < 10. All input files required to reproduce such models are included in data file S1.


Supplementary material for this article is available at

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 P. Ryan and R. Guy for help in deploying the Wallowa2 array. We thank S. King and an anonymous reviewer for the careful and constructive suggestions. We also thank F.-C. Lin for providing us with the crustal anisotropy measurements for California. Y.K. acknowledges support from the Creative Pioneering Researchers Program of Seoul National University (SNU SRnD 3345-20160014). The figures presented in this paper were made using the Generic Mapping Tools version 4.5.9 ( We used HC 1.0.5 (50, 5759) published under the GPL2 license. Funding: This project was supported by NSF/EAR 1546635 and 1547594. Author contributions: J.C.C. produced the anisotropy results and cowrote the paper. J.P.-H. modeled the coupling between mantle and crust and cowrote the paper. Y.K., A.C.S., and B.N. synthesized the receiver functions, seismic tomography, and SKS results. R.W.C. and E.H. conceptualized the project and edited the paper. All authors contributed to the interpretation of the results and the preparation of the manuscript for publication. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The project used data from High Lava Plains, IDOR, Wallowa, and Wallowa2 surveys, along with data from the T.A. (Transportable Array) from EarthScope. All data are available from the IRIS-DMC ( 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