Many studies have shown that typical oceanic lithosphere is underlain by a well-developed asthenosphere characterized by slow seismic velocities from ~100 to 250 km depth. However, the fate of the oceanic asthenosphere at subduction zones is poorly understood. I show here using shear-wave splitting of S waves emanating from earthquakes in the Juan de Fuca slab that upper mantle asthenospheric anisotropy beneath the slab is consistent with the presence of two distinct subducted asthenospheric layers, one with fast shear trends parallel to the subduction trench, and a second, deeper layer with fast upper mantle fabrics parallel to the motion of the Juan de Fuca plate with respect to the deeper mantle. The consistent orientation of unsubducted Pacific asthenospheric anisotropy in the direction of current plate motion implies that the trench-parallel, subslab anisotropy develops when the lithosphere subducts.
SEISMIC ANISOTROPY IN THE OCEANIC ASTHENOSPHERE
Many studies have shown that typical oceanic lithosphere is underlain by a well-developed asthenosphere characterized by slow seismic velocities from ~100 to 250 km depth (Helmberger, 1973; Dziewonski and Anderson, 1981; Kennett et al., 1995). Seismic anisotropy of the upper mantle is nearly ubiquitous, and, resulting from linear preferred orientation of olivine and pyroxene aggregates (Mainprice et al., 2000), is often viewed as a gauge of the sense and degree of shear coupling between moving lithospheric plates and underlying mantle mesosphere (Vinnik et al., 1989; Russo and Okal, 1998; Wolfe and Silver, 1998; Garcia and Russo, 2005). In the Pacific Ocean basin, the upper mantle at asthenospheric depths is well characterized by studies of surface-wave azimuthal anisotropy (Nishimura and Forsyth, 1988, 1989; Montagner and Tanimoto, 1990, 1991; Montagner, 2002; Maggi et al., 2006) and to a certain degree by shear-wave splitting observations at island seismic stations (Russo and Okal, 1998; Wolfe and Silver, 1998; Fontaine et al., 2005; Garcia and Russo, 2005): away from localized perturbations due to recent hotspot activity (e.g., Montagner, 2002; Maggi et al., 2006), fast shear trends of upper mantle anisotropy generally parallel current plate motions, indicating that the Pacific oceanic asthenosphere is a zone of simple shear decoupling the lithospheric plates from underlying mantle (Vinnik et al., 1989).
However, the fate of the oceanic asthenosphere at subduction zones is poorly understood. I show here using shear-wave splitting of S waves emanating from earthquakes in the Juan de Fuca slab that upper mantle asthenospheric anisotropy beneath the slab is consistent with the presence of two distinct subducted asthenospheric layers, one with fast shear trends parallel to the subduction trench, and a second, deeper layer with fast upper mantle fabrics parallel to the motion of the Juan de Fuca plate with respect to the deeper mantle. The consistent orientation of unsubducted Pacific asthenospheric anisotropy in the direction of current plate motion (Nishimura and Forsyth, 1988, 1989; Montagner and Tanimoto, 1990, 1991; Montagner, 2002; Maggi et al., 2006) implies that the trench-parallel, subslab anisotropy develops when the lithosphere subducts.
SUBDUCTION OF OCEANIC ASTHENOSPHERE BENEATH THE CASCADES
A straightforward test of the form of subslab asthenospheric anisotropy—a direct measure of subducted asthenospheric shear fabrics—can be applied where the Juan de Fuca plate subducts beneath North America (Fig. 1). The subduction zone is suitable for our purpose since it is relatively short along strike but is significantly curved, changing strike from nearly N-S south of 46°N in Oregon to NNE, then N, and finally to NW beneath Vancouver Island (Crosson and Owens, 1987; Bostock and VanDecar, 1995; Xue and Allen, 2007; Sigloch et al., 2008), and thus the potential for three-dimensional upper mantle flow can be assessed in a relatively small volume.
The major limitation of such a study is the scarcity of suitable Juan de Fuca slab earthquakes since the advent of the modern global seismic networks (i.e., Global Seismic Network [GSN], GEOSCOPE, and GEOFON) in the 1990s. These events must be large enough to be well recorded teleseismically (Mw of 5.5 at a minimum) and deeper than 40 km to ensure that near-source surface reflections pS and sS do not contaminate the direct S wave. In fact, only two events meet these criteria (Fig. 1 and Table 1), but given the wide coverage—azimuthally and in distance—afforded by the current global networks, and the compact, strongly three-dimensional nature of the subduction zone, the sampling of the subslab region by shear waves emanating from these two events is reasonably complete.
ISOLATION OF SOURCE-SIDE, SHEAR-WAVE SPLITTING
In order to isolate shear-wave splitting that occurs beneath the Juan de Fuca slab, an important consideration is where along path splitting may occur. Teleseismic S waves originating in subducting slabs turn in the lower mantle, so potential contributions to splitting could conceivably occur in the upper or lower mantle during downward propagation, in the core mantle boundary region, if the source-station distance is great enough (Δ > 80°), and in the lower or upper mantle on the upgoing path. Shear-wave splitting at the core-mantle boundary (CMB) region appears to be small (e.g., Hall et al., 2004), and the lower mantle itself is devoid of demonstrable splitting (Meade et al., 1995). Therefore, S waves at these distant stations are unlikely to be contaminated by splitting at the CMB, and the lower mantle contributes effectively nothing to source-side splitting thus isolated. The upper mantle below the olivine stability field (i.e., deeper than 410 km) has been shown to be a source of shear-wave splitting in several instances (e.g., Fouch and Fischer, 1996; Fischer and Wiens, 1996; Wookey and Kendall, 2004), attributed to preferred alignment of the higher pressure olivine polymorph wadsleyite by Fouch and Fischer (1996), but appears to be the exception rather than the rule. However, shear-wave splitting is nearly ubiquitous at stations recording core-traversing waves (e.g., SK(K)S), and it is necessary to correct for splitting that occurs in the upper mantle beneath the receiver station (Fig. 2).
Upper mantle anisotropy beneath the Juan de Fuca slab can thus be evaluated by isolating shear-wave splitting that occurs in the event source region as waves travel downward from the slab to distant stations (Fig. 2). Contributions to observed splitting in the source region are then determined by correcting for known splitting due to upper mantle anisotropy beneath the receiver stations (Russo and Silver, 1994). I used published shear-wave splitting parameters for a group of broadband seismic stations at the requisite teleseismic distances (30° < Δ < 84°) from the two Juan de Fuca slab events detailed in Table 1. Station locations and splitting parameters, and citations of the source publications from whence they are drawn, can be found in GSA Data Repository Table DR1.1 Station splitting corrections were applied to data used in this study in all cases when clear splitting was observed at the receiver stations; however, if no splitting was observed, I did not perform a receiver-side correction, as a splitting null after passage through the upper mantle on both the downgoing and upgoing legs indicates such correction is unwarranted and more than likely to result in spurious splitting. Note that only one receiver station used in this study is considered to be isotropic (PAYG; Fontaine et al., 2005), and thus never required correction for receiver-side splitting.
Numerical modeling studies of anisotropic wave propagation indicate that a systematic reorientation of apparent shear-wave splitting parameters occurs during passage of shear waves through regions with variable anisotropy (Saltzer et al., 2000), and the observed splitting fast shear azimuth, Φ, closely parallels the fast propagation axis of the last anisotropic layer traversed. Thus, to isolate splitting from beneath the Juan de Fuca slab, we need waves that sample this anisotropic region last on their way down through the isotropic lower mantle. Downgoing S waves propagating from the two slab events are therefore appropriate. Receiver-side splitting should affect upgoing S waves from a given back azimuth similarly, so corrections made using apparent splitting parameters for stations with no clear back-azimuthal splitting dependence isolate the subslab splitting contributions we seek (Figs. 2 and 3,302303304). I used the method of Russo and Silver (1994) to correct for known station splitting; details of the method can be found in the supplementary information (see footnote 1).
Application of the receiver-station splitting corrections detailed in Table DR1 (see footnote 1) yielded 22 clear source-side splitting measurements for event 99184 (3 July 1999) and 30 null splitting observations, where either the S-wave polarization parallels one of two possible anisotropic symmetry directions or the medium is isotropic. For the Nisqually earthquake (event 01059, 28 February 2001), I made 19 clear source-side splitting measurements and 33 null splitting observations. Details of the splitting nulls are shown in Figures DR2 and DR3 and Table DR2 (see footnote 1). Measurements of source-side, shear-wave splitting made with receiver corrections are shown in Figure 1 (see also Fig. DR2 and Table DR2 [footnote 1]). Resulting split S waves from the two events show a systematic variation in receiver-corrected apparent splitting parameters, depending on source-receiver takeoff azimuth (Fig. 1). For paths that sample the region immediately below the Juan de Fuca slab in the direction of slab dip (eastward azimuths) and along the slab strike (azimuths to NW and S-SE), apparent splitting parameters are systematically parallel to the subduction trench in the direction of wave propagation. For paths that sample the subslab region along W and WSW azimuths (i.e., at high angle to the slab strike and ~180° from slab dip), the apparent splitting parameters are close to the current velocity azimuth of the Juan de Fuca plate with respect to assumed-fixed hotspots (Fig. 1). Furthermore, one or the other of the allowed fast shear azimuths in the large number of splitting nulls observed for the two events is consistent with the observed splitting for source-receiver paths that sample similar subslab regions (see Figs. DR2 and DR3 [see footnote 1]).
The number of null-splitting observations for the two events is significant but not unusual. Note that strong splitting is frequently observed (mean delay time is 2.9 s) at a given station, but the S wave for a given event at a neighboring station is not apparently split (null). The S-wave travel paths for these stations are very similar, so it is unlikely that the subslab region includes large-scale isotropic pockets. Note also that in some instances stations recording null splitting for one event show clear splits for the other. Splitting nulls result from a combination of focal mechanism initial S-wave polarization (nodal for S radiation); source-receiver and anisotropic symmetry geometries; and event source complexity (the larger magnitude Nisqually event yielded significantly more nulls in part due to apparent rupture complexity [Ichinose et al., 2004], which complicates waveforms and may change S polarization during rupture); potentially because of subtle errors in receiver-station corrections; because of suppression of the split signal by attenuation or high station noise; and possibly due to complex heterogeneous anisotropy somewhere along the path. Nonetheless, the observed source-side splitting is largely consistent between the two events.
Subslab Sampling Differences
In order to demonstrate the distinct sampling of the source-side, shear-wave splitting, I traced rays through a Cartesian volume. Shear-wave velocities within the volume were assigned based on published variations from standard reference models for the area (e.g., Crosson and Owens, 1987; Bostock and VanDecar, 1995), which allowed construction of a Cascades “slab” anomaly 1%–2% higher in velocity than surroundings. No slow-velocity anomalies were included. Rays were traced along event-station azimuths from the event hypocenters down to 900 km depth. The principal result of the ray tracing was to show that the interpretation of two anisotropic layers is tenable, and to explain the apparent inconsistency of splitting results for rays taking off along NW azimuths. As shown in Figures 4,–,6, rays departing from the two earthquakes do sample significantly different volumes of the subslab region in this direction: blue rays from the 3 July 1999 event, with trench-parallel splitting, sample a region just below the slab that red rays from the 28 February 2001 event do not traverse. The latter rays cross a region more internal to the core of the plunging fold formed by the subduction zone. The interpretation is then that these red rays sample the Juan de Fuca absolute plate motion (APM)–parallel anisotropy layer.
Source-Side Splitting Versus Cascades Region Station Splitting
Note that a large majority of the source-side splitting results are inconsistent with observations of splitting of SK(K)S and PKS phases made at stations deployed above the Juan de Fuca slab, which consistently parallel the hotspot-defined absolute plate motions of the Juan de Fuca and North American plates (Silver and Chan, 1991; Bostock and Cassidy, 1995; Currie et al., 2004), ~N70°E; see Figure 1. As the SK(K)S and PKS waves traverse the entire anisotropic region of the upper mantle, from below slab to surface, it is likely that the observed differences with source-side, S-wave splitting are due to differences in seismic anisotropy above the slab (SKS splits), where seemingly strong anisotropy has developed parallel to the shear on the interplate interface in the upper mantle wedge, and below the slab (S splits), where anisotropy is apparently three-dimensional. Note that the source-side S waves do not traverse—and therefore cannot sample—anisotropy in the region above the Juan de Fuca slab. Despite the differences in the splitting of these two types of shear waves, we do observe a small number of source-side splits that parallel the Juan de Fuca APM direction (Fig. 1). These results can be explained as a consequence of three-dimensional flow below the curved slab, as detailed below.
Shear-Wave Splitting and Upper Mantle Anisotropy
Once made, splitting measurements require consideration of upper mantle anisotropy for interpretation. Teleseismic shear-wave splitting is commonly associated with development of linear preferred orientation (LPO) of primarily upper mantle olivine, with a tendency for aggregates to align in the shear plane parallel to the direction of tectonic extension (Gueguen and Nicolas, 1980; Christensen, 1984; Nicolas and Christensen, 1987; Ribe, 1989a, 1989b; Ribe and Yu, 1991; Zhang et al., 2000). Laboratory experiments on natural and synthetic aggregates representative of dry upper mantle compositions indicate that when such rocks yield shear-wave splitting, observed fast polarization directions parallel the  crystallographic axes of the aligned minerals (e.g., Hess, 1964; Carter et al., 1972; Kaminski and Ribe, 2001), the type-A fabric of Jung et al. (2006). Compilations of natural upper mantle samples from locales around the globe indicate that the predominant form of upper mantle anisotropy follows this basic type-A fabric (Mainprice and Silver, 1992; Ben Ismail and Mainprice, 1998). Jung et al. (2006) also note the existence of type-C and type-E fabrics that are petrographically distinct but yield shear-wave splitting with fast polarizations in the material flow direction, similar to that produced by the A-type fabric (see their Fig. 12). The C- and E-type fabrics are found in hydrous synthetic olivine aggregates.
Factors thought to complicate anisotropic fabrics at the scale of mineral grains or hand-sample aggregates of upper mantle minerals include presence of water under high stress conditions, which yield a distinctive B-type fabric (Jung and Karato, 2001; Karato, 2003; Jung et al., 2006), and presence of melt (Holtzman et al., 2003). Elevated water content and partial melt consistent with conditions in the mantle wedge above subducted slabs may modify LPO fabrics, yielding b-axis concentrations or girdles of crystallographic axes instead of the usual a-axis clustering. Numerical models meant to represent deformation processes on a larger scale (tens of kilometers) have indicated that episodes of noncoaxial finite strain may yield more complicated anisotropic fabrics (Tommasi et al., 1999; Kaminski and Ribe, 2001; Blackman and Kendall, 2002). However, direct observation of Pacific basin shear-wave splitting is inconsistent with the strongly variable asthenospheric anisotropy predicted by these models (e.g., Hess, 1964; Nishimura and Forsyth, 1988, 1989; Montagner and Tanimoto, 1991; Russo and Okal, 1998; Garcia and Russo, 2005). The asthenosphere beneath Pacific basin plates appears to be largely a zone of simple shear with fast shear trends in the local absolute plate motion direction (simple asthenospheric flow model of Vinnik et al. ).
Given that the asthenosphere beneath Pacific plates likely was thoroughly devolatilized by ridge processes, the asthenospheric channel beneath the Juan de Fuca slab is unlikely to include significant water. Although slab dewatering is an important process at subduction zones (Peacock et al., 2005; Screaton and Saffer, 2005), fluid migrates from the slab upward into the suprasubduction mantle wedge following the pressure gradient decrease. Thus, although this fluid hydrates the mantle wedge and likely pervasively modifies LPO (Jung and Karato, 2001; Karato, 2003; Jung et al., 2006; Bostock et al., 2002), there is no evidence that mantle beneath the slab is affected at all. Furthermore, given their cooling histories, a melt fraction at the base of subducted lithosphere or below—in the absence of any breach in the slab, or flow around a slab edge—is highly unlikely (e.g., Currie et al., 2004; Minshull et al., 2005; Peacock et al., 2005). Available tomographic studies have not been interpreted to show partial melt beneath the Juan de Fuca slab (e.g., Xue and Allen, 2007; Sigloch et al., 2008). Thus, shear-wave splitting results at this subduction zone can be linked to A-type anisotropy except in the mantle wedge region, where the increased likelihood of hydration and partial melt fraction make the B-type fabric a potential alternative (e.g., Abt and Fischer, 2008). Since the downgoing S waves used in this study do not sample the suprasubduction mantle wedge, we adopt the a-axis olivine LPO model for all interpretations.
Two Layers of Subslab Anisotropy
We rule out relating the source-side, S-wave splitting results to an anisotropic symmetry that plunges in the slab dip direction (see Fig. DR4 [see footnote 1]), as might be invoked for flow at subduction zones, because the potential patterns of 2π periodicity in apparent splitting parameters versus source-station azimuth for shape-preferred–orientation anisotropies (Crampin and Booth, 1985) and for olivine LPO (Chevrot and van der Hilst, 2003) differ markedly from the observations for fast shear azimuths and delay times of both events. Instead, as shown in Figure 7, the variation of apparent splitting with source-receiver geometry can be explained by a two-layer model of upper mantle fabrics beneath the Juan de Fuca slab: a layer immediately below the slab has predominant anisotropy with fast shear parallel to the local trench strike; and a deeper layer beneath the western, as-yet unsubducted portion of the Juan de Fuca plate, is characterized by fast shear in the Juan de Fuca absolute motion direction, i.e., ENE (Gripp and Gordon, 2002) (Fig. 7). The two layers are sampled heterogeneously by S waves emanating from the two slab source events, such that waves propagating to the west will sample an increasing portion of even a uniform-thickness, deeper, APM anisotropy layer, given the steeply ENE-plunging geometry of the curved Juan de Fuca slab. In contrast, waves propagating in other takeoff directions preferentially sample the trench-parallel anisotropy immediately beneath the slab. Three-dimensional views of the ray sampling are shown in Figures 4–6, which demonstrate that rays from the two events traverse the subslab region comprehensively, and that apparent inconsistency of results for takeoff azimuths to the NW is due to the differences in ray geometry resulting from the ~56 km distance between the event hypocenters.
Thus, observations of the asthenosphere beneath oceanic plates are consistent with simple asthenospheric flow, but, near the Juan de Fuca subduction zone, flow apparently becomes three-dimensional and a trench-parallel flow layer develops (Fig. 7). Given the geometry of sampling by S waves leaving the study region, the subslab trench-parallel flow layer cannot be more than ~120 km thick, although the oblique ray paths traveled by the downgoing waves are much longer, consistent with a nominal 3%–4% anisotropy to yield the long delay times observed (mean of 2.9 s). This trench-parallel layer may be a manifestation of strain partitioning in the subslab upper mantle due to interaction of the subducted slab with the overriding plate: compressional stresses across the subduction zone apparently yield a reorientation of subslab LPO. That the trench-parallel anisotropy layer is proximal to the slab, whereas the more distal APM anisotropic layer lies deeper, would be consistent with flow around the southern edge of the Juan de Fuca–Gorda slab, as has recently been proposed to account for regional-scale, shear-wave splitting measurements in western North America (Zandt and Humphreys, 2008). A global survey of subduction zone anisotropy showing ubiquitous, strong subslab fabrics (Long and Silver, 2008) is consistent with our interpretation, and our result is similar to those of recent studies detailing the variation in upper mantle flow and anisotropic fabrics above slabs in subduction-zone mantle wedges, which also show development of adjacent domains of nearly orthogonal anisotropy (Abt and Fischer, 2008; Hoernle et al., 2008).
Subduction of oceanic lithosphere appears to have a strong effect on upper mantle fabrics of the underlying asthenosphere. Whereas the oceanic asthenosphere beneath the unsubducted Pacific basin plates appears to be largely a zone of simple shear with fast shear trends parallel to current APM directions, the upper mantle beneath the Juan de Fuca slab is more complex: two layers of strong seismic anisotropy exist beneath the slab, one immediately beneath the slab characterized by strong trench-parallel anisotropic fabrics, and a second, deeper layer that maintains the Juan de Fuca APM anisotropic fabric presumably formed prior to subduction. The existence of the trench-parallel anisotropic layer at the subduction zone, but not in the Pacific basin to the west, implies that this layer formed as a consequence of subduction. Strain partitioning in the subducting asthenospheric upper mantle due to Juan de Fuca–North America convergence and compression across the subduction zone is clearly one possible mechanism for developing such a layer.
The shear-wave splitting database published online by the group at the Laboratoire de Tectonophysique, Université de Montpellier II, was indispensable in finding suitable stations for receiver-side correction (see www.gm.univ-montp2.fr/splitting/DB/). I am grateful to Mark Panning, David Foster, Karin Sigloch, Gene Humphreys, Matt Fouch, and David James for helpful discussions, and to an anonymous reviewer whose thorough, probing commentary resulted in significant improvement to the manuscript. Thanks also to Editor Jon Pelletier for seeing the work through review and publication. All data used in this study were downloaded from the Incorporated Research Institutions for Seismology (IRIS) Data Management Center.