Abstract

The Isabella anomaly, a prominent upper-mantle high-speed P-wave anomaly located within the southern Great Valley and southwestern foothills of the Sierra Nevada, has been interpreted either as foundering sub-Sierran lithosphere or as remnant oceanic lithosphere. We used Vp/Vs anisotropy tomography to distinguish among the probable origins of the Isabella anomaly. S waveforms were rotated into the Sierran SKSFast and SKSSlow directions determined from SKS-splitting studies. Teleseismic P-, SFast-, SSlow-, SKSFast-, and SKSSlow-wave arrival times were then inverted to obtain three-dimensional (3-D) perturbations in Vp, Vp/VsMean, and percent azimuthal anisotropy using three surface wave 3-D starting models and one one-dimensional (1-D) model. We observed the highest Vp/Vs anomalies associated with slower velocities in regions marked by young volcanism, with the largest of these anomalies being the Mono anomaly under the Long Valley region, which extends to depths of at least 75 km. Peak Vp/Vs perturbations of +4% were found at 40 km depth. The low velocities and high Vp/Vs values of this anomaly could be related to partial melt.

The high wave speeds of the Isabella anomaly coincide with low Vp/Vs values with peak perturbations of −2%, yet they do not covary spatially. The P-wave inversion imaged the Isabella anomaly as a unimodal eastward-plunging body. However, the volume of that Isabella anomaly contains three separate bodies as defined by varying Vp/Vs values. High speeds, regionally average Vp/Vs values (higher than the other two anomalies), and lower anisotropy characterize the core of the Isabella anomaly. The western and shallowest part has high wave speeds and lower Vp/Vs values than the surrounding mantle. The eastern and deepest part of the anomaly also contains high speeds and lower Vp/Vs values but exhibits higher anisotropy. We considered combinations of varying temperature, Mg content (melt depletion), or modal garnet to reproduce our observations. Our results suggest that the displaced garnet-rich mafic root of the Mesozoic Sierra Nevada batholith is found in the core of the Isabella anomaly. If remnant oceanic lithosphere exists within the Isabella anomaly, it most likely resides in the shallow, westernmost feature.

Within the Sierra Nevada, the highest upper-mantle anisotropy is largely contained within the central portion of the range and the adjacent Great Valley. Anisotropy along the Sierra crest is shallow and confined to the lithosphere between 20 and 40 km depth. Directly below, there is a zone of low anisotropy (from 170 to 220 km depth), low velocities, and high Vp/Vs values. These features suggest the presence of vertically upwelling asthenosphere and consequent horizontal flow at shallower depths. High anisotropy beneath the adjacent western foothills and Great Valley is found at ∼120 km depth and could represent localized mantle deformation produced as asthenosphere filled in a slab gap.

INTRODUCTION

The seismic velocity structure beneath the Sierra Nevada (Fig. 1) has been an important constraint in decoding the uplift history of this range and the geodynamic consequences of uplift on surrounding areas. One widely discussed hypothesis for uplift of the Sierra Nevada is through removal of mafic mantle lithosphere and/or dense lower crust. While the exact mechanisms are still debated, removal of dense, garnet-rich material beneath the Sierra Nevada would result in convective upwelling of asthenosphere directly beneath and east of the range (e.g., Ducea and Saleeby, 1998). Upward migration of hot, buoyant asthenosphere would lead to crustal thinning, range uplift, and volcanism—all of which have been inferred from xenolith, geologic, seismic, and magnetotelluric studies (e.g., Ducea and Saleeby, 1996; Jones et al., 1994; Park, 2004).

Figure 1.

Overview map of study region. Physiographic provinces are labeled in blue. Black lines denote plate boundaries, including the San Andreas fault (thick line), Mendocino fracture (thin line), and the Cascadia subduction zone (teeth). Blue polygons outline the high-velocity seismic anomalies at 170 km depth using the P-wave tomographic results from Jones et al. (2014) with the “Moschetti-free” starting model. Abbreviations: ML—Mono Lake; LVC—Long Valley caldera; SAF—San Andreas fault; IA—Isabella anomaly; RA—Redding anomaly; G—Gorda plate; JdF—Juan de Fuca plate; MM—Monterey microplate.

Figure 1.

Overview map of study region. Physiographic provinces are labeled in blue. Black lines denote plate boundaries, including the San Andreas fault (thick line), Mendocino fracture (thin line), and the Cascadia subduction zone (teeth). Blue polygons outline the high-velocity seismic anomalies at 170 km depth using the P-wave tomographic results from Jones et al. (2014) with the “Moschetti-free” starting model. Abbreviations: ML—Mono Lake; LVC—Long Valley caldera; SAF—San Andreas fault; IA—Isabella anomaly; RA—Redding anomaly; G—Gorda plate; JdF—Juan de Fuca plate; MM—Monterey microplate.

The distribution of upper-mantle seismic velocities is somewhat consistent with this hypothesis. The western foothills of the Sierra Nevada exhibit fast lower-crust/uppermost-mantle seismic velocities relative to surrounding material (e.g., Biasi and Humphreys, 1992; Fliedner et al., 1996; Jones et al., 1994) and thick crust, up to 55 km (Frassetto et al., 2011). The crust below the Sierra crest has slower velocities and thins to 35 km, which is comparable to crustal thicknesses in the Basin and Range Province. West of the southern Sierra Nevada, there lies an isolated, high-velocity, upper-mantle anomaly, termed the Isabella anomaly (Fig. 1). The magnitude of the Isabella anomaly is large, up to ∼10% faster in P-wave speed than surrounding material (e.g., Jones et al., 2014). The Isabella anomaly is an important feature in unraveling the evolutionary history of the Sierra Nevada.

There has been debate over the Isabella anomaly’s origin, geometry, and composition since it was first identified by Raikes (1980). In east-west cross sections, the inferred geometry has varied from being almost vertical (Benz and Zandt, 1993; Biasi and Humphreys, 1992; Zandt and Carrigan, 1993; Jones et al., 1994) to steeply east-plunging (Jones et al., 2014) to more gently east-plunging (e.g., Boyd et al., 2004; Yang and Forsyth, 2006). The Isabella anomaly has also been inferred to have originated as sub-Tehachapi Mountains lithosphere (Jones et al., 1994), Sierran lithosphere (Biasi and Humphreys, 1992), and a fragment of Farallon oceanic lithosphere called the Monterey microplate (Benz and Zandt, 1993; see also Fig. 1, MM, herein). Thus, the composition of the Isabella anomaly has variously been interpreted as being continental lithospheric mantle, solely Sierran lower crust and upper mantle, or dehydrated, melt-depleted oceanic lithosphere of the Farallon slab (e.g., Wang et al., 2013; Pikser et al., 2012; Cox et al., 2016). A clearer understanding of the seismic velocity heterogeneities in this region is needed in order to distinguish among the upper-mantle processes that have been proposed.

Because velocity heterogeneities can arise from variations in temperature, composition, melt fraction, mineral phase, and hydration, identifying the physical state of the subsurface demands more than just observations of a single kind of wave speed. In addition, spurious isotropic velocity variations can emerge from anisotropy when raypaths are not uniformly distributed, as is typically the case for teleseismic tomography. For example, Bezada et al. (2016) observed significant anisotropy-induced artifacts in isotropic P-wave tomogram simulations of subduction zones. These artifacts were more prominent for slow-velocity anomalies. However, these authors also observed that using a uniform distribution of events could exacerbate these effects. Sobolev et al. (1999) observed the largest anisotropy-induced artifacts coming from anisotropic material with a plunging olivine a-axis. Because temperature has a more pronounced effect on S waves than P waves (e.g., Goes et al., 2000), shear-wave speeds help to isolate thermal effects from other processes.

Seismic stations that are near the Isabella anomaly have generally exhibited low SKS splitting delay times (Fig. 2B; Yang et al., 2016; Liu et al., 2014; Table S1 in the Supplemental Materials1) that might reflect compositional variations or structural fabrics. Boyd et al. (2004) suggested that the low-shear-wave splitting can either be attributed to the presence of garnet, which is isotropic, or to a vertically oriented fast axis, which reduces splitting for near-vertical incoming S waves. If the low SKS splits are due to the presence of garnet, then it might be possible to resolve differences between garnet and ultramafic lithologies using Vp/Vs and anisotropy.

Figure 2.

(A) Seismic stations used in this study. Red circles refer to the Sierra Nevada Earthscope Project (SNEP) network. Blue circles refer to the Sierran Paradox Experiment (SPE) network. White circles are other stations used in this study from a variety of networks including Berkeley, Caltech, Transportable Array, USArray/Terrascope, University of Oregon, University of Nevada, and Leo Brady. (B) Map of SKS splits from Liu et al. (2014; black), Polet and Kanamori (2002; magenta), and new measurements (Table S2 [footnote 1]; blue).

Figure 2.

(A) Seismic stations used in this study. Red circles refer to the Sierra Nevada Earthscope Project (SNEP) network. Blue circles refer to the Sierran Paradox Experiment (SPE) network. White circles are other stations used in this study from a variety of networks including Berkeley, Caltech, Transportable Array, USArray/Terrascope, University of Oregon, University of Nevada, and Leo Brady. (B) Map of SKS splits from Liu et al. (2014; black), Polet and Kanamori (2002; magenta), and new measurements (Table S2 [footnote 1]; blue).

In this study, we utilized Vp/Vs anisotropy tomography to better separate the roles of temperature and composition in generating the high-velocity Isabella anomaly and other features. This approach allowed us to address unresolved questions, such as: Is the Isabella anomaly seismically fast because it is composed of melt-depleted, oceanic material or because it represents cold, possibly garnet-rich lithospheric material? Progressive extraction of partial melt will result in an increase in Mg# [Mg/(Mg + Fe) × 100] as iron is lost to melt. This effect, termed melt depletion, will decrease density, Vp/Vs, and increase seismic velocities, but an increase in garnet content will increase density, seismic velocities, and Vp/Vs (e.g., Boyd et al., 2004).

Our approach to estimating Vp/Vs differs from most other similar studies in this region (e.g., Schmandt and Humphreys, 2010a, 2010b). Teleseismic shear waves can be partitioned into two orthogonal pulses traveling with different speeds that are determined by the mineralogical orientation and/or structural fabric of a region. This shear-wave splitting can introduce errors in traveltime measurements if splitting varies within a region. We accounted for anisotropy in our study by measuring both S- and SKS-wave traveltimes on components rotated parallel and perpendicular to the regionally consistent SKS-fast direction (Fig. 2B; e.g., Boyd et al., 2004; Schutt and Humphreys, 2004).

We inverted for lateral variations in the magnitude of azimuthal anisotropy, defined as hexagonal symmetry with a horizontal symmetry axis, to further constrain the tectonics beneath the Sierra Nevada region. Most teleseismic anisotropy is introduced in the upper mantle, where it is dominantly controlled by the presence of olivine and is a result of strain-induced lattice preferred mineral orientation. In contrast, more mafic lithologies such as eclogites typically exhibit very little anisotropy (e.g., Fountain and Christensen, 1989). By using measurements on both teleseismic S and SKS arrival times, we gained depth resolution from raypaths with different incidence angles, compared to the uniformly steep incidence angles of SKS waves used in shear-wave splitting studies (e.g., Silver and Chan, 1991). The shear waves used in this study are well within the shear-wave window, with incidence angles less than 35° (i.e., Savage, 1999). A discussion of the limitations in our anisotropy tomography approach can be found in the Results section.

DATA AND METHODS

Data

This study expands upon an earlier P-wave study of the Sierra Nevada (Jones et al., 2014). The P-wave arrival times used here were described by Jones et al. (2014) and do not take into account anisotropic effects on P-wave traveltimes. S- and SKS-wave arrival times were measured from earthquakes recorded on 389 seismometers deployed within California, western Nevada, southern Oregon, and southwestern Idaho as part of nine temporary and permanent seismic networks from the time periods of June to October 1997 and May 2005 to September 2007 (Fig. 2A; Table S2 [footnote 1]). Of particular focus were the temporary deployments within the Sierra Nevada. The Sierra Nevada Earthscope Project (SNEP; Owens et al., 2005) network spanned the northern to central portion of the range and was deployed from May 2005 to September 2007. The Sierran Paradox Experiment (SPE; Jones and Phinney, 1997) covered the southern portion of the range and was deployed from June to October 1997. The SNEP and SPE network stations were closely spaced, at roughly 25 km, compared to the more broadly spaced Transportable Array (TA) and other permanent network stations in the Sierra Nevada. This close spacing allowed for finer imaging of upper-mantle features beneath the range.

Teleseismic S and SKS phases used in this study came from 104 M >5.5 earthquakes with epicentral distances ranging from 30° to 80° and from 96° to 120° (Fig. 3; Table S3).

Figure 3.

Map of seismic events used for the tomographic inversion. SFast arrival times (left) were picked independently of SSlow (right) and thus have a different event distribution.

Figure 3.

Map of seismic events used for the tomographic inversion. SFast arrival times (left) were picked independently of SSlow (right) and thus have a different event distribution.

Typical S-wave tomography studies transform the three component waveforms from the vertical-north-east (ZNE) coordinate system to vertical-radial-transverse (ZRT) system to maximize S-wave energy arriving in the transverse component, i.e., the first-arriving SH waves. Waveform coherence breaks down in the presence of lateral variations in azimuthal anisotropy, which can lead to the transverse component containing variously shifted combinations of the original S wave (e.g., Schutt and Humphreys, 2004). We reduced traveltime uncertainty due to anisotropic effects by rotating all S and SKS waveforms into the Sierran SKSFast (N75°E) and SKSSlow (N15°W) directions observed within the Sierra Nevada (Fig. 2B; Table S1). We note that rotating all shear waveforms across the study area to one coordinate system unique to the Sierra Nevada may produce artifacts in regions that do not exhibit the same fast orientation. Accordingly, our anisotropy tomography interpretations are restricted to the Sierra Nevada.

For events arriving from the east-northeast and west-southwest, S- and SKS-wave arrivals had too low of a signal-to-noise ratio for one of the two (fast and slow) components. We postulate that this arose from a combination of factors: (1) fewer M 5.5+ teleseisms exist in the east-northeast back-azimuth range, with arrivals typically obscured by noise; and (2) shear waves initially polarized into the fast and slow directions will not split into two components (e.g., Savage, 1999).

Data Processing

S-waveform processing maximized arrival signal over noise. S waveforms were first filtered with a Butterworth band-pass filter from 0.03 to 0.5 Hz. Events with noisier waveforms were additionally band-pass filtered from 0.01 to 0.1 Hz. SFast and SSlow arrival times were separately and independently picked using dbxcor, a waveform cross-correlation, beam-forming program (Fig. 4; Pavlis and Vernon, 2010). Seismograms for one event were cross-correlated to create a beamed waveform. An arrival time was then picked on the beamed waveform and appropriately adjusted in time to each individual waveform. This yielded absolute arrival times for each station-event pair. All residuals for one event were de-meaned before the inversion, removing any systematic time shifts arising from poor picks on the beamed waveform. SFast residuals were de-meaned separately from SSlow residuals. P-wave arrival times from Jones et al. (2014) were used in this study. In total, 29,182 P-, 7169 SFast-, and 5255 SSlow-wave arrival time measurements were inverted to obtain three-dimensional (3-D) perturbations in Vp, Vp/VsMean, and anisotropy.

Figure 4.

SFast and SSlow waveforms for a Mw 7.5 event originating from northern Peru. White lines indicate the shear-wave arrival time determined from dbxcor.

Figure 4.

SFast and SSlow waveforms for a Mw 7.5 event originating from northern Peru. White lines indicate the shear-wave arrival time determined from dbxcor.

A separate set of events was analyzed for the single-station shear-wave splitting measurements we report in Table S1 (footnote 1). SKS phases for earthquakes of mb ≥ 5.5 occurring at epicentral distances of ≥88° were inspected for the duration of the SNEP deployment, along with a few selected events from the 1997 SPE experiment. Seismograms were filtered prior to splitting analysis using a zero-phase Butterworth band-pass filter with corner frequencies of 0.04 and 0.3 Hz. Splitting parameters were constrained using the semi-automated approach of Teanby et al. (2004), which is based on the Silver and Chan (1991) approach. Horizontal components were rotated and time-shifted to minimize the second eigenvalue of the covariance matrix for particle motion within a time window around the SKS pulse. This is equivalent to linearizing the particle motion and minimizing the tangential component of shear-wave energy. The Silver and Chan (1991) approach takes a single, manually picked, shear-wave analysis window. In the cluster analysis approach we adopted here, however, the splitting analysis is performed for a range of window lengths, and cluster analysis is utilized to find measurements that are stable over many different windows. All splitting parameters were determined after analysis of 100 different windows. An example of the analysis is shown in Figure S1.

Traveltime Residuals

Early Arrivals Related to High-Velocity Bodies

The pattern and amplitude of early SFast- and SSlow-wave residuals varied with back azimuth and incidence angle consistent with an upper-mantle source for most of these variations. Similar to P-wave residuals from Jones et al. (2014), the earliest S-wave arrivals were those that traversed the Isabella anomaly. SFast and SSlow residuals from rays traversing the Isabella anomaly were earliest from south-southeastern back azimuths (120°–180°; Fig. 5I). Arrivals in this back-azimuth range were up to +3 s earlier than predicted and thus could travel within the Isabella anomaly longer, an effect that has also been inferred for early arrivals traversing the Gorda slab in northern California (Benz et al., 1992).

Figure 5.

Map of SFast and SSlow residuals with respect to a predicted arrival time using the one-dimensional (1-D) IASP91 seismic model (Kennett and Engdahl, 1991; interactive). Early (late) arrivals are denoted in cooler (warmer) colors. To interact with figure, click on figure to activate. Move mouse across hemisphere map (in the lower right) to display results from a particular back azimuth (Baz). If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.5i to interact with Figure 5.

Figure 5.

Map of SFast and SSlow residuals with respect to a predicted arrival time using the one-dimensional (1-D) IASP91 seismic model (Kennett and Engdahl, 1991; interactive). Early (late) arrivals are denoted in cooler (warmer) colors. To interact with figure, click on figure to activate. Move mouse across hemisphere map (in the lower right) to display results from a particular back azimuth (Baz). If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.5i to interact with Figure 5.

Delayed Arrivals

Delayed arrivals were most prominent along the California-Nevada border, east of the Sierra Nevada. These arrivals extended as far north as 40°N in some back azimuths (Fig. 5I). However, the largest cluster of delays was found east of the Sierra Nevada around 38°N, near the Long Valley caldera (Fig. 1, LVC) and Mono-Inyo Craters region (Fig. 5I). While arrivals at stations within the Long Valley region were consistently late, the magnitudes varied with incidence angle and back azimuth. Delays were most significant and coherent for SKS phases (epicentral distances greater than 96°; Fig. 5I). Variations in these delayed arrivals with back azimuth and distance were not as prominent as those related to the Isabella anomaly, suggesting a shallower source. However, the slight azimuthal variations are sufficient to preclude a purely crustal origin for this feature.

Inversion Parameterization

We modified an existing iterative, spherical, finite-difference (FD) tomography code, sphfd (Roecker et al., 2006), to include separate fast and slow shear-wave arrival times. This code performs a weighted, damped least squares inversion designed to optimize the trade-off between data and model misfit through appropriate damping and smoothing parameters. Our model space was parameterized by a set of nodes laid on a rectangular grid at specified depths. This inversion scheme incorporated the use of coarse and fine grids, where the coarse grid was used for the traveltime inversion, and the finer grid was used for calculating raypaths and traveltimes. The coarse grid spacing was an integer multiple of the fine grid’s spacing. Similar to Jones et al. (2014), our coarse interval was set to ∼25 km within our region of interest. Grid spacing increased toward the edges of our model to allow for teleseisms to enter at the base of our model and not on the sides. Nodes were placed at depths of −30, 0, 20, 40, 70, 120, 170, 220, 270, 320, and 430 km below sea level. P-, SFast-, and SSlow-wave arrival times were inverted to obtain 3-D perturbations in Vp, Vp/VsMean, and anisotropy, defined in this study as: 
graphic
 
graphic
 
graphic
 
graphic
where A is anisotropy (though half of a more common definition of anisotropy). SFast and SSlow velocities were determined from Equations 3 and 4. Note that we inverted for Vp/VsMean instead of Vs, a choice that avoids spurious Vp/Vs anomalies where either Vp or Vs observations are insufficient to constrain that wave-speed model.

Iterative Scheme

Traveltimes were computed using a finite difference solution to the eikonal equation modified from Hole and Zelt (1995). Two sets of traveltimes were calculated. First, times from each station to the bottom edge of the model space were calculated. Next, traveltimes from each earthquake to the base of the model were computed and then added to the previous times. The earliest arrivals were used as the predicted P- and S-wave arrival times. The partial derivatives of traveltime with respect to slowness were determined separately for SFast and SSlow and related changes in slowness, Vp/VsMean, and anisotropy to traveltime residuals. SFast and SSlow partial derivatives relating our traveltime observables to our model parameters at each grid point, i, along a raypath were calculated as: 
graphic
 
graphic
where ψP, ψSFast, and ψSSlow are, respectively, the slowness values for P, SFast, and SSlow.

The partial derivatives (Equations 5 and 6) that were calculated for the fine grid were then accumulated for the coarse grid by summing the contributions of the fine grid elements.

The system of linear equations was inverted using the LSQR algorithm of Paige and Saunders (1982). Damping factors were used in the inversion to limit the amount of changes in a single iteration (with respect to the previous iteration’s model) by accepting a poorer fit to the data. One global damper was applied to all three variables. This damper was set to 250 during the first eight iterations and then was subsequently ramped down to 150 by the fourteenth iteration. These damping values were determined from Jones et al. (2014). However, using one damper for three different variables (with three different units) suppresses the results for Vp/Vs and anisotropy. Therefore, Vp/Vs and anisotropy factors were rescaled by factors of 4 and 3, respectively, in the relevant columns in the partial derivative matrix. These values are in proportion to the expected difference in standard deviation of those terms relative to the P-wave results.

For each iteration, the perturbations from the inversion were smoothed twice using a moving-window average of three nodes in the x, y, and z directions. The smoothing operator limited variations with respect to adjacent nodes. The perturbations were then added to that iteration’s starting model. Each inversion ran for 14 iterations.

Similar to Jones et al. (2014), we employed both one-dimensional (1-D) and 3-D starting models to separate features that were common across all models from those that were artifacts arising from specific starting models, which exposed the uncertainty from the limited resolution of the teleseismic inversion in the uppermost parts of the model. For each starting model, anisotropy was set to zero. We used IASP91 (Kennett and Engdahl, 1991) as our 1-D starting model and three regional surface wave models for our 3-D models (Moschetti et al., 2010; Gilbert et al., 2012; Shen et al., 2013). Gilbert et al. (2012) inverted Rayleigh waves for shear-wave structure using teleseisms recorded at SNEP, SPE, TA, and other permanent stations within California. Moschetti et al. (2010) produced a 3-D isotropic shear-wave model of the western United States by inverting both Rayleigh and Love wave dispersion measurements from ambient noise and teleseisms recorded by TA stations. Shen et al. (2013) created a 3-D shear-wave model of the central and western United States by inverting Rayleigh wave phase velocity maps from ambient noise and earthquakes with receiver functions using TA stations. All 3-D starting models were resampled to match the grid dimensions used in this study. The starting P-wave velocities were obtained by converting the 3-D surface wave shear models using the empirical relations between Vp and Vs from Brocher (2005). Because there is a variation of Vp/Vs with Vs in this work, our starting models do have variations in Vp/Vs. Wave speeds were averaged to keep the vertical transit time of both P and S waves consistent between the original model and our model parameterization. Given that the lateral resolution from surface waves decreases with depth, we limited the use of the P- and S-wave velocities from these models to a depth of 170 km. P- and S-wave velocities from IASP91 were used to fill out the rest of our model space from 220 to 430 km depth.

Scaled Anisotropy

Our inverted parameter A only captures the relative variations in azimuthal anisotropy between our predetermined fast and slow axes. That is, A is insensitive to the absolute average azimuthal anisotropy. We therefore present a modified version of this parameter to more closely match the absolute magnitude of azimuthal anisotropy. Instead of showing anisotropy perturbations centered on 0, a constant is added to A that reflects values normally used to describe upper-mantle anisotropy. Percent anisotropy is more commonly expressed as the difference between SFast and SSlow velocity over the mean: 
graphic
Our inversion is not sensitive to the actual anisotropy but rather to the deviations (D) from the mean (U) fast and slow velocities for each layer, where 
graphic
Thus, 
graphic
which can be recombined as 
graphic
The value that is measured in our inversion is 
graphic
where Ur is the mean S-wave velocity in the layer, which should equal (UFast+ USlow)/2, and Dmean is the mean deviation at that node in both fast and slow velocities, so Dmean = (DFast+ DSlow)/2.
Assuming this to be true, we find 
graphic
Applying a simple Taylor expansion, where we expand around Dmean = 0 and ignore the Dmean2 term, gives us 
graphic
If we define the average azimuthal anisotropy in a layer as T = (UFastUSlow)/Ur, then we get 
graphic

For simplicity, T was assigned an arbitrary value of 4.5% for all depths. Previous SKS splitting studies within California have assumed an upper-mantle anisotropy of 4% when determining the depth extent of anisotropic layers (Polet and Kanamori, 2002; Özalaybey and Savage, 1995). Xie et al. (2015) observed intrinsic (called “inherent” in their study) S-wave anisotropy from 4% to 5% within the upper mantle of the western United States.

Only the real inversion results (in map view and cross section) are displayed using the scaled values (C). Synthetic tests were performed and are displayed using the original, unscaled values (A).

Model Resolution

To test the resolving power of the inversion code to image variations independent of anisotropy, we created a checkerboard model with nonzero values confined to nodes at 170 km depth. Square anomalies 75 km wide (∼4 nodes) were assigned a peak P-wave perturbation of ±5%, a peak Vp/Vsmean perturbation of ±5%, and no anisotropy. A noiseless synthetic traveltime data set was made using this checkerboard model and the real source-receiver pairs. We then inverted the synthetic traveltime data using the same grid geometry and parameters as in our real inversion. After fifteen iterations, the peak values of the recovered anomalies were 63% of the starting model for P and ∼50% for Vp/Vsmean, and anisotropy values remained under 1% (Fig. 6). The traveltime variance of our synthetic data set was reduced by ∼95% after the fifteenth iteration. Although anisotropy should be zero for this model, anomalies as large as 0.7% were observed in the eastern portions of our model space. The underestimation of peak magnitudes resulted mainly from the vertical smearing of compact anomalies into adjacent depth layers. However, vertically integrating these anomalies largely recovered the starting magnitudes, especially more toward the center of the anomalies and less near the edges. For P-wave slowness perturbations, recovered magnitudes were similar to Jones et al. (2014), ranging from 50% to 140%. The lowest integral was found in the corner nodes of the anomalies, and the highest integral was found in the center of the anomalies. Recovered integrated perturbations in Vp/Vs ranged from ∼17% to 140%. The high integral was also found in the core of the anomaly. The corner nodes varied from 17% (at the lowest) in the southwest-northeast nodes to ∼60% for the northwest-southeast corner nodes.

Figure 6.

Checkerboard synthetic tests at 170 km depth. Initial model only varies Vp and Vp/Vs in this layer. Final model recovers pattern of Vp and Vp/Vs anomalies but underestimates magnitudes. The presence of anisotropy anomalies in the final model was found where station sampling was poorer.

Figure 6.

Checkerboard synthetic tests at 170 km depth. Initial model only varies Vp and Vp/Vs in this layer. Final model recovers pattern of Vp and Vp/Vs anomalies but underestimates magnitudes. The presence of anisotropy anomalies in the final model was found where station sampling was poorer.

A similar test was conducted to test the stability of the code to separate anisotropy from isotropic variations. For this second test, another checkerboard test was created, and variations were again confined to 170 km depth, with no P-wave perturbation, no Vp/Vsmean perturbation, and ±5% perturbation in anisotropy. The final inversion recovered up to 54% of the peak starting anisotropy anomaly (Fig. 7). For this model, P-wave and Vp/VsMean perturbations should be 0. However, our final model included up to 1% P-wave and 0.3% Vp/VsMean perturbations. The pattern of the Vp/VsMean variations resembled the anisotropy checkerboard, implying that some cross-talk does exist between the two parameters. The large errant signals in P-wave perturbations were only observed in the easternmost portion of our model, where ray coverage was minimal owing to fewer stations. Traveltime variance was reduced by 94%–95% after the final iteration. Vertically integrated values for anisotropy indicated that recovered magnitudes ranged from 30% (corner nodes) to 140% (center nodes). Similar to Vp/Vs, though not as dramatic, northwest-southeast corner nodes recovered generally ∼30% of the starting anomalies, and southwest-northeast corner nodes recovered ∼70% of the starting anomalies.

Figure 7.

Checkerboard synthetic tests at 170 km depth. Initial model only varies anisotropy. Final model recovers the pattern of anisotropy anomalies but underestimates their magnitudes. The presence of positive polarity anomalies (green/blue) in Vp indicates less stability toward the edge of our model where sampling resolution is reduced. Vp/Vs results exhibit up to 0.3% anomalies mapped into the model space.

Figure 7.

Checkerboard synthetic tests at 170 km depth. Initial model only varies anisotropy. Final model recovers the pattern of anisotropy anomalies but underestimates their magnitudes. The presence of positive polarity anomalies (green/blue) in Vp indicates less stability toward the edge of our model where sampling resolution is reduced. Vp/Vs results exhibit up to 0.3% anomalies mapped into the model space.

Synthetic Vp, Vp/VsMean, and anisotropy anomalies were better recovered and spurious anomalies were minimized within and near the Sierra Nevada. The closer station spacing in this region allowed for better ray coverage. Vertical resolution for both tests was poor compared to lateral resolution, owing to the steep raypaths of teleseismic P and S phases, especially SKS. Consequently, vertical smearing of compact anomalies originally at 170 km extended from ∼100 to ∼200 km. However, the largest recovered amplitudes were still located at 170 km depth. As was the case with checkerboard synthetic tests, our runs revealed that our values for Vp, Vp/VsMean, and anisotropy can underestimate the true magnitude of anomalies. In our tests, underestimations were at most a factor of 2 within the core region of interest for this study. For a more detailed discussion of resolution and smearing effects using the SNEP and SPE networks, refer to Jones et al. (2014).

We also tested whether excluding anisotropy from the inversion reduced the fit to the data. For this test, we reran our inversion on the observed arrival times using the Moschetti starting model. Anisotropy was damped by a factor of 30 for the first seven iterations and then released on the eighth iteration to levels used in the real inversion. Although the results for Vp, Vsmean, and Vp/Vsmean were similar to our normally damped results, the SFast- and SSlow-wave data variances over the first seven iterations were higher than those produced by the inversion including anisotropy. After seven iterations, the test SFast variance was 0.415 s2, compared to the SFast variance (after seven iterations) of 0.387 s2 when anisotropy was allowed. The test SSlow variance was 0.374 s2, while the normally damped SSlow variance was 0.290 s2. Once we returned the anisotropy damper to values from the regular inversion at iteration 8, the difference in data variance between this test and the normally damped inversion became negligible (Fig. 8). The pattern of high and low azimuthal anisotropy remained consistent between the final test and real models. However, the final test model generally underestimated peak magnitudes. These results confirm our earlier synthetic test inversion of a checkerboard without anisotropy; both indicated that inferred variations in anisotropy are not artifacts from misfitting isotropic structures but are features required by observations that separately reduce the data variance by ∼7% for SFast and over 20% for SSlow.

Figure 8.

Data variance comparison between full inversion and isotropic test inversion. “Full inversion” includes anisotropy during all fourteen iterations. “Test inversion” only introduces anisotropy at the eighth iteration.

Figure 8.

Data variance comparison between full inversion and isotropic test inversion. “Full inversion” includes anisotropy during all fourteen iterations. “Test inversion” only introduces anisotropy at the eighth iteration.

RESULTS

Broadly speaking, inversions initiated from each of the four starting models fit the arrival times nearly equally well (Table 1). For all of our starting models, after 14 iterations, the variance in traveltime residuals was reduced by ∼70% for P, ∼60% for SFast, and ∼70% for SSlow from starting values.

TABLE 1.

DATA VARIANCE

Comparison of Results from Different Starting Models

The four different starting models yielded somewhat different final solutions, with the differences reflecting nonunique aspects of the inversion. While the pattern of anomalous features was consistent across the different models, there were variations in the shape and magnitude of these features. However, these variations were more profound for P- and S-wave velocities than for Vp/VsMean and anisotropy. For example, the shape of the high-velocity Isabella anomaly varied substantially between the Shen and Moschetti models. For both inversions, however, the Vp/VsMean and anisotropy results exhibited minor differences in shape (Fig. 9). The minimal differences in anisotropy between the various starting models is presumably because of the identical starting conditions (no anisotropy) for all the inversions; thus, the absence of differences in anisotropy between different inversions does not guarantee robustness. In contrast, because Vp/Vs does vary between the starting models, the similarity of results below 20 km depth probably does provide an indication of the robustness of our results.

Figure 9.

Comparison of east-west cross sections across the Isabella anomaly using different starting models as labeled at the top of each column. Anisotropy results are from the scaled version. High velocity and anisotropy are denoted by cooler colors. High Vp/Vs is denoted by warmer colors. Outlines of subanomalies IA1–IA3 are shown; in interactive version, button at lower-left toggles the overlay on and off. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.9i to interact with Figure 9.

Figure 9.

Comparison of east-west cross sections across the Isabella anomaly using different starting models as labeled at the top of each column. Anisotropy results are from the scaled version. High velocity and anisotropy are denoted by cooler colors. High Vp/Vs is denoted by warmer colors. Outlines of subanomalies IA1–IA3 are shown; in interactive version, button at lower-left toggles the overlay on and off. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.9i to interact with Figure 9.

Five Vp/VsMean anomalies were consistent across the different starting models. Broadly speaking, low Vp/VsMean anomalies were found in regions characterized by fast velocities related to the Isabella anomaly and the Redding anomaly, a high-velocity feature in northern California (RA in Fig. 1; see also Fig. 10I). This observation agrees with other Vp/Vs tomographic studies (Boyd et al., 2004; Schmandt and Humphreys, 2010a, 2010b). There were also three prominent high Vp/VsMean features that were coincident with regions characterized by lower velocities (Fig. 10I). These Vp/VsMean highs coincided with regions of known Holocene and latest Quaternary volcanism: the Clear Lake volcanic field, the Coso volcanic field, and the Long Valley caldera and Mono-Inyo Craters region.

Figure 10.

Tomographic results in map view by starting model (interactive). P-wave perturbations are in left panel, Vp/VsMean perturbations are in middle panel, and scaled anisotropy is in the right panel. To interact with figure, first click on figure to activate and then click on depth and model to display results. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.10i to interact with Figure 10.

Figure 10.

Tomographic results in map view by starting model (interactive). P-wave perturbations are in left panel, Vp/VsMean perturbations are in middle panel, and scaled anisotropy is in the right panel. To interact with figure, first click on figure to activate and then click on depth and model to display results. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.10i to interact with Figure 10.

Isabella Anomaly

The previously defined P-wave Isabella anomaly encloses variations in both Vp/VsMean and azimuthal anisotropy, which suggest it is a polygenetic body. While the Isabella anomaly still retains a fairly compact, unimodal higher-than-average P- and S-wave speed anomaly, east-west cross sections through it surprisingly reveal two low Vp/VsMean anomalies separated by a region with higher, albeit background, Vp/Vs values (Fig. 9). This division resembles the variations in Vp/Vs observed within the Isabella anomaly by Boyd et al. (2004) but differs considerably in the patterns of anisotropy.

The Vp/VsMean lows were located at different depths. The shallower and westernmost anomaly, defined as IA1, exhibited peak Vp/VsMean perturbations ranging from −1% to −3% between 40 and 70 km depth. The strength of this signal varied with starting model and was strongest using the Shen starting model (Fig. 9). This anomaly was weakest using the IASP91 starting model. Because the IASP91 1-D model does not include information about Great Valley sediments, low wave speeds that should be located closer to the surface are smeared into deeper depths by the limited vertical resolution of the tomography at shallow depths, which distorts the shallower IA1 Vp/Vs feature.

The deepest low Vp/VsMean anomaly, defined as IA3, was located beneath the eastern portion of the Sierra Nevada between 170 and 220 km depth. Vp/VsMean perturbations were as low as ∼–2%. This anomaly is within the easternmost location of the Isabella anomaly (Fig. 9). While still faster than surrounding mantle, this portion of the Isabella anomaly was not as fast as IA1 or the central feature found at ∼200 km depth. Instead, peak P- and S-wave velocities ranged from 1% to 4% faster than surroundings. There was minimal variation in the shape and magnitude of this Vp/VsMean feature across the starting models, an indication that ray coverage was sufficient at this depth and location. There were also no variations between starting models because surface wave models did not resolve anomalies at these depths, and the ray coverage prevented shallower differences from bleeding into these depths (Fig. 9).

Finally, the core of the Isabella anomaly observed in the P-wave tomograms (Jones et al., 2014), defined as IA2, corresponded to a region with a low-to-average Vp/VsMean ratio (Fig. 9) that was higher than the eastern (deep) and western (shallow) anomalies.

Anisotropy also varied between the different Vp/VsMean bodies (Fig. 9). The easternmost and deepest Vp/VsMean anomaly, IA3, exhibited greater anisotropy than the other two anomalies, but was comparable to background levels. IA2 (core) exhibited the lowest anisotropy among the three features. IA1 (westernmost) coincided with a questionable region of low anisotropy that could result from streaking of IA2 structures. IA1 was the least resolved of the three anomalies owing to a lack of seismic stations in the San Joaquin Valley and its lower crust–shallow mantle location, where results from differing starting models diverge. Vertical smearing effects were observed in our anisotropic checkerboard synthetic tests, so we do not interpret the low anisotropy for the shallowest part of the Isabella anomaly. The descriptions of the Isabella anomalies (IA1–IA3) are summarized in Table 2.

TABLE 2.

ISABELLA ANOMALIES FROM THE DIFFERENT STARTING MODELS

Because the Vp/Vs contrasts suggest that the Isabella anomaly exhibits variations in physical properties, we tested whether the multicomponent Vp/VsMean anomaly was an artifact introduced into our model by limited ray coverage or limited resolution (Figs. 11 and 12). First, we tested whether a unimodal, plunging velocity anomaly like that seen in the P-wave tomography could create two Vp/VsMean low anomalies (Fig. 11). Synthetic traveltime data were created from a starting model characterized by a rectangular plunging anomaly with peak Vp/Vs perturbations of −4% and peak P-wave perturbations of 3.7%. This anomaly was placed between 90 and 200 km, plunging 45° to the east. The final results did not reproduce the multifeature Vp/Vs anomaly seen in our real data. Instead, we recovered a single plunging anomaly with up to 59% of the peak starting Vp/Vs magnitude. Second, we tested whether two Vp/VsMean anomalies could create a singular velocity anomaly (Fig. 12). Two Vp/Vs anomalies were given a peak Vp/Vs perturbation of −4%. The shallower anomaly was placed between 50 and 90 km depth, while the deeper anomaly was placed from 160 to 200 km depth. Although we assigned the same magnitude perturbation to both starting anomalies, final inversions revealed that the deeper Vp/Vs anomaly was better recovered, up to 41% of the starting magnitude was recovered at the peak depth. The shallower anomaly only recovered up to 25% of the starting magnitude at its peak depth. Results from this test reveal that if two high-velocity anomalies were separated by slower ambient mantle, we would resolve these two bodies. Neither test reproduced the observed relationship between P- and S-wave velocity and Vp/VsMean in our results, indicating that such features are not artifacts. Additional synthetic tests revealed that the Vp-Vp/Vs pattern observed within the Isabella anomaly can arise from a starting model where IA2 is characterized by higher Vp and higher Vp/Vs than IA1 and IA3 (Fig. S2 [footnote 1]). While the spatial distribution of anomalies was not a result of ray geometry artifacts, some signal from a P-wave velocity anomaly might have influenced the Vp/Vs magnitude (Fig. S3).

Figure 11.

Testing Isabella anomaly structures: One body. Synthetic anomaly test determining whether station spacing is sufficient to properly recover a single, plunging Isabella anomaly. The final results, although exhibiting vertical smearing, indicate that our resolution is sufficient to prevent a uniform body from being imaged as two separate bodies. Cross-section results using the Moschetti starting model have been added to the third column for comparison.

Figure 11.

Testing Isabella anomaly structures: One body. Synthetic anomaly test determining whether station spacing is sufficient to properly recover a single, plunging Isabella anomaly. The final results, although exhibiting vertical smearing, indicate that our resolution is sufficient to prevent a uniform body from being imaged as two separate bodies. Cross-section results using the Moschetti starting model have been added to the third column for comparison.

Figure 12.

Testing Isabella anomaly structures: Multiple bodies. Synthetic anomaly test determining whether our sampling resolution is sufficient to resolve two high-velocity anomalies separated by background mantle. Final results recover the two distinct anomalies in the final inversion. Cross-section results using the Moschetti starting model have been added to the third column for comparison.

Figure 12.

Testing Isabella anomaly structures: Multiple bodies. Synthetic anomaly test determining whether our sampling resolution is sufficient to resolve two high-velocity anomalies separated by background mantle. Final results recover the two distinct anomalies in the final inversion. Cross-section results using the Moschetti starting model have been added to the third column for comparison.

Mono, Clear Lake, and Coso Anomalies

Our results also contained three high Vp/Vs regions that were all located in areas marked by young volcanism. These locations are the Clear Lake volcanic field, the Coso volcanic field, and the Mono-Inyo magmatic region (Fig. 10I). The Coso Vp/VsMean anomaly, located east of the Sierra Nevada near 36°N, was up to 3% higher than average (Fig. 13). This feature corresponds with a zone of low velocities, down to −6% for P waves. While the velocity anomaly was more compact, mainly present above 50 km, the Vp/VsMean anomaly extended below 50 km. This Vp/VsMean anomaly also exhibited low anisotropy.

Figure 13.

East-west cross section across the high Vp/Vs anomaly near the Coso volcanic field.

Figure 13.

East-west cross section across the high Vp/Vs anomaly near the Coso volcanic field.

The Clear Lake Vp/VsMean anomaly was located away from the SNEP/SPE networks, in a poorly sampled part of the study area. Consequently, only the broadest aspects of this feature will be discussed. The peak Vp/VsMean signal was consistently located at 40 km depth, across all starting models, with amplitudes up to 3% (Fig. 14). The shape of the Vp/VsMean was also consistent across the starting models.

Figure 14.

East-west cross section across the high Vp/Vs anomaly near the Clear Lake volcanic field.

Figure 14.

East-west cross section across the high Vp/Vs anomaly near the Clear Lake volcanic field.

The Vp/VsMean anomaly within the Mono-Inyo region was the strongest of the three volcanic anomalies, with amplitudes as high as ∼3.9% at depths from 40 to 70 km depth. Many magmatic features are found in this region, including the Long Valley caldera and the Mono-Inyo craters. We term this the Mono anomaly for simplicity. Perturbations of 3% extended to depths of 100 km (Fig. 15). The shape of the Vp/VsMean feature was more equant than the corresponding P-wave anomaly. Within the Vp/VsMean anomaly, P-wave velocities were ∼3% slower than surrounding material. The Mono Vp/VsMean anomaly was large, exceeding 75 km in diameter. In east-west cross sections, the anomaly extended to depths of 170 km, though features below 100 km could represent smearing artifacts. Anisotropy related to the Mono anomaly was comparable to background values (Fig. 15).

Figure 15.

Southwest-northeast cross section through the Mono high Vp/Vs anomaly, located on the eastern side of the central Sierra Nevada. High Vp/Vs is denoted in warmer colors, whereas high velocity and anisotropy are denoted in cooler colors.

Figure 15.

Southwest-northeast cross section through the Mono high Vp/Vs anomaly, located on the eastern side of the central Sierra Nevada. High Vp/Vs is denoted in warmer colors, whereas high velocity and anisotropy are denoted in cooler colors.

The pattern of high Vp/VsMean within the Mono-Inyo region, if real, should be observed in the raw P and S traveltime residuals. To create a high Vp/Vs anomaly, the S-wave residual needs to be significantly more delayed than appropriately scaled P-wave raw residuals. We tested the null hypothesis that anomalies in the area have a constant Vp/Vs ratio with the data by rescaling the P residuals and comparing them with the observed S-wave residuals. We chose the 10 May 2006 Mw 6.4 event from the Aleutian Islands to perform our test, since Mono delays were particularly strong for west-northwest back azimuths. P-wave residuals at each station were multiplied by an upper-mantle Vp/Vs value of 1.81, which should mimic S-wave residuals in the absence of a Vp/Vs anomaly under this region. We then used station SNEP11 as a reference station and subtracted its value from all other stations to emphasize local variations in Vp/Vs (Fig. 16). This hypothetical S residual was then compared to the actual SFast and SSlow traveltime residuals. While the Mono-Inyo region generally did exhibit more delays in P and S traveltime residuals, the hypothetical S residual (Px1.81) delays were smaller than actually observed in either SFast or SSlow. The largest hypothetical delay, relative to SNEP11, in the Mono-Inyo region was ∼2 s. Within the real S residuals, the largest relative delays (to SNEP11) were up to 5 s for SFast and up to 6 s for SSlow (Fig. 16).

Figure 16.

P-, SFast-, and SSlow-wave residuals for the M6.4 Fox Island (Aleutian Islands) earthquake on 10 May 2006 with a back azimuth (BAZ) of 309°. P residuals were multiplied by a Vp/Vs value of 1.81 to provide the null hypothesis that there is no significant Vp/Vs anomaly. P residual × 1.81, SFast, and SSlow residuals were compared to that of a reference station (SNP11, black star). P residuals × 1.81 are smaller than the magnitudes exhibited by the S residuals, thus demonstrating the need for higher Vp/Vs ratios along rays penetrating under the Mono volcanic region.

Figure 16.

P-, SFast-, and SSlow-wave residuals for the M6.4 Fox Island (Aleutian Islands) earthquake on 10 May 2006 with a back azimuth (BAZ) of 309°. P residuals were multiplied by a Vp/Vs value of 1.81 to provide the null hypothesis that there is no significant Vp/Vs anomaly. P residual × 1.81, SFast, and SSlow residuals were compared to that of a reference station (SNP11, black star). P residuals × 1.81 are smaller than the magnitudes exhibited by the S residuals, thus demonstrating the need for higher Vp/Vs ratios along rays penetrating under the Mono volcanic region.

We modeled the Mono Vp/VsMean anomaly to determine the minimal extent of this anomaly. The Mono anomaly was modeled as a cube 50 km on a side extending from 40 to 90 km depth. The background Vp/Vs was set to 1.75, and the cube was given a Vp/Vs of 1.79. This resulted in a peak Vp/Vs perturbation of 2.3% between the anomaly and the surroundings. After fourteen iterations, our inversion recovered 40% of the peak amplitude at 40 km depth. We found that ∼95% of the integrated anomaly was recovered between 0 and 170 km depth (Fig. 17). The peak amplitude was preserved between 40 and 90 km depth, and the overall pattern appeared similar to our actual results (cf. Fig. 17 with Fig. 15). Not surprisingly, the lateral resolution was better than the vertical resolution.

Figure 17.

Synthetic test of the Mono Vp/Vs anomaly across a southwest-northeast cross section. Top: Starting anomaly Vp/Vs perturbation of +2%. Bottom: Recovered model underestimates the true amplitudes and vertically smears past 100 km. However, the position of the peak anomaly is conserved.

Figure 17.

Synthetic test of the Mono Vp/Vs anomaly across a southwest-northeast cross section. Top: Starting anomaly Vp/Vs perturbation of +2%. Bottom: Recovered model underestimates the true amplitudes and vertically smears past 100 km. However, the position of the peak anomaly is conserved.

Anisotropy in the Sierra Nevada

Our estimates of variations in upper-mantle anisotropy across our region of interest bear some similarity to published SKS splitting delay magnitudes (Özalaybey and Savage, 1995; Polet and Kanamori, 2002; Hongsresawat et al., 2015), irrespective of the orientation of the fast direction (Fig. 18). Predicted splitting delay times, made from our anisotropy results, exhibit some agreement with SKS delay times. Predicted delay times for the Sierra Nevada region were shifted by 1 s to account for systematic time shifts in our calculations and to match magnitudes from other studies. Regions of high anisotropy (from our study) that coincided with high SKS splitting delay times up to 2 s included (1) the central Great Valley and western foothills of the Sierra Nevada around 38°N; (2) an east-west band extending from northeastern California into northwestern Nevada at ∼41°N; and (3) a localized region east of the Sierran crest at ∼37°N and ∼118°W (Figs. 18 and 19). Three prominent low-anisotropy regions that were associated with low SKS splitting delay times included (1) near 40°N, 122°W, (2) east of the Sierran crest at 36°N, and (3) within the Isabella anomaly. While the anisotropy results showed some vertical smearing, anisotropy anomalies peaked at well-defined depths.

Figure 18.

Left: Comparison of integrated anisotropy tomography results (from 0 to 270 km depth) to existing SKS splitting measurements. Anisotropy tomogram uses Moschetti et al. (2010) as the starting model. Circles represent stations with SKS splitting measurements. Measurements from Polet and Kanamori (2002), Liu et al. (2014), Yang et al. (2016), and this paper (Table S1 [text footnote 1]) were averaged for each station. Circle color indicates whether averaged SKS-fast direction is closely parallel to Sierran SKS-fast direction (N75°E). Size of circle denotes magnitude of average splitting delay- time. Checker pattern highlights regions where SKS fast directions are “misoriented” (station fast-direction is more than 20° away from Sierran SKS fast direction). Right: Comparison of predicted splitting delay time (using our results) to existing SKS splits. Negative delay times are generally found where our choice of projected coordinate system falls apart. Nonetheless, the magnitudes of the delays are comparable to those found in SKS splitting measurements. We added a 1 s time shift to the model to align model results (which are insensitive to absolute splitting magnitude) with observed splitting delay magnitudes.

Figure 18.

Left: Comparison of integrated anisotropy tomography results (from 0 to 270 km depth) to existing SKS splitting measurements. Anisotropy tomogram uses Moschetti et al. (2010) as the starting model. Circles represent stations with SKS splitting measurements. Measurements from Polet and Kanamori (2002), Liu et al. (2014), Yang et al. (2016), and this paper (Table S1 [text footnote 1]) were averaged for each station. Circle color indicates whether averaged SKS-fast direction is closely parallel to Sierran SKS-fast direction (N75°E). Size of circle denotes magnitude of average splitting delay- time. Checker pattern highlights regions where SKS fast directions are “misoriented” (station fast-direction is more than 20° away from Sierran SKS fast direction). Right: Comparison of predicted splitting delay time (using our results) to existing SKS splits. Negative delay times are generally found where our choice of projected coordinate system falls apart. Nonetheless, the magnitudes of the delays are comparable to those found in SKS splitting measurements. We added a 1 s time shift to the model to align model results (which are insensitive to absolute splitting magnitude) with observed splitting delay magnitudes.

Figure 19.

Anisotropy tomography results in map view and in cross section. Map-view tomograms are from 40 and 170 km depth using the Moschetti et al. (2010) starting model. Cross sections run along (A-A′) and across (B-B′) the Sierra Nevada.

Figure 19.

Anisotropy tomography results in map view and in cross section. Map-view tomograms are from 40 and 170 km depth using the Moschetti et al. (2010) starting model. Cross sections run along (A-A′) and across (B-B′) the Sierra Nevada.

Within the vicinity of the Sierra Nevada, anisotropy was highest in the central portion of the range and westward into the Great Valley. The central Sierran anisotropic high was found at shallow depths from 20 to 40 km. Below this layer, there was a zone of low anisotropy from 120 to 220 km depth that peaked at 170 km depth (Fig. 19). This anisotropic low coincided with high Vp/Vs and low Vp.

The high anisotropy observed along the Great Valley and western Sierra foothills was generally deeper than the central Sierran feature. However, peak values related to this feature were shallower to the south. Around 38°N, high anisotropy was found at 120 km depth, while at 37°N, high anisotropy extended from 20 to 70 km depth (Fig. 19).

SKS fast orientations were not consistent across our study region (Fig. 18), which means that anisotropy values are ambiguous where SKS fast directions are misoriented with respect to our chosen reference. SKS fast directions in the Coast Ranges are oriented northwest to southeast (∼N60°W), where San Andreas fault motion dominates. The northern Sierra Nevada–southern Cascade region around 39.5°N also hosts SKS fast directions oriented northwest to southeast (Fig. 18). We thus limited our interpretations to the Sierra Nevada and Great Valley regions.

We tested the resolvability of an anisotropic feature in the upper mantle near the position of the Isabella anomaly. A rectangular anomaly 100 km wide and 80 km thick was centered at 120 km depth and was given a peak anisotropic anomaly of 5% with no variation in Vp or Vsmean (Fig. 20). This initial rectangular feature smeared vertically to become more elongated in the vertical direction (Fig. 20). However, the lateral dimensions of the initial anomaly were well preserved. Although we were not able to recover the limited vertical extent of the anomaly, we observed the highest amplitudes at the center of the starting anomaly, at 120 km. At this depth, 35% of the starting anomaly was recovered. If we included magnitudes from the adjacent depth layers, 70 and 170 km, we recovered 91% of the starting anomaly. Our synthetic tests indicate that while the vertical extent of the anomaly is poorly constrained, and while our amplitudes represent underestimations of real azimuthal anisotropy, we can recover the mean depth and lateral extent of anisotropic features.

Figure 20.

Testing the robustness of anisotropic material at depth. A high (+5%) anisotropy anomaly at 120 km depth exhibits vertical smearing and amplitude underestimation, but the depth of the peak anomaly is still recovered.

Figure 20.

Testing the robustness of anisotropic material at depth. A high (+5%) anisotropy anomaly at 120 km depth exhibits vertical smearing and amplitude underestimation, but the depth of the peak anomaly is still recovered.

DISCUSSION

Multicomponent Isabella Anomaly

We identified three distinct Vp/Vs anomalies that coincide with the Isabella P-velocity anomaly, which synthetic tests have shown are not artifacts. Therefore, the Isabella anomaly is not homogeneous from top to bottom and thus cannot be explained by variation in one physical attribute or one composition. A segmented Vp/Vs feature associated with the Isabella anomaly is similar to the results of Boyd et al. (2004), whose tomographic study employed similar methods but used a far smaller data set and focused only on the region near the Isabella anomaly. Vp/Vs cross sections through the Isabella anomaly from Schmandt and Humphreys (2010a, 2010b) did not resolve multiple bodies. This discrepancy might result from assuming isotropy both in measuring S-wave times and inverting them, from differing smoothing procedures, or from a different inversion approach.

Creating a High-Velocity Isabella Anomaly

Elastic Calculations

Many physical processes can create a high-velocity Isabella anomaly, faster than ambient mantle: colder temperatures, dehydration, phase transitions (specifically from basalt to eclogite), Mg enrichment (also referred to as melt depletion), orthopyroxene enrichment, and garnet enrichment. We compared our Vp, Vp/VsMean, and anisotropy results with rock physical properties calculated using the Abers and Hacker (2016) workbook to narrow down the probable causes of the Isabella anomaly.

Elastic properties for various ultramafic compositions were calculated at 3 GPa. Compositions included oceanic upper mantle as well as continental lower crust and upper mantle. For oceanic upper mantle, we used the starting modal proportions for harzburgite, enriched harzburgite, and lherzolite from Hacker et al. (2003). Depleted harzburgite was also included. For depleted compositions, elastic properties were calculated using only the Mg-rich mineralogical end members (forsterite and diopside, for example). Temperature varied from 300 °C to 1300 °C to simulate thermal variations that would be observed at any point within the oceanic slab. However, it is unlikely that temperatures colder than 700 °C would be observed in oceanic lithosphere at ∼90 km depth; these values are included purely for completeness.

For Sierran lower-crust and upper-mantle compositions, we varied temperature from 900 °C to 1300 °C, similar to Boyd et al. (2004). We used garnet-pyroxenite (eclogite), garnet peridotite, and spinel peridotite mineral compositions defined in Boyd et al. (2004). Spinel peridotite calculations were performed at 1.5 GPa, since higher pressures would result in the transformation from spinel peridotite to garnet peridotite. We also included four garnet-pyroxenite (eclogite/“arclogite”) and garnet websterite modal compositions from central Sierra Nevada xenoliths collected at Big Creek (Lee et al., 2006) to investigate how much garnet, if present, could realistically be within the Isabella anomaly. Lee et al. (2006) divided garnet-pyroxenite xenoliths into two groups based on MgO content. The low-MgO xenoliths were composed of more than 50% garnet, clinopyroxene, and lesser amounts of orthopyroxene. One sample contained up to 77% garnet. The high-MgO xenoliths were dominantly composed of clinopyroxene (>50%), less garnet (<50%), and some orthopyroxene. Rock compositions used in these calculations are summarized in Table 3 and Figure 21.

TABLE 3.

MINERALOGICAL COMPOSITIONS FOR LITHOLOGIES USED IN THE ROCK PHYSICAL PROPERTIES CALCULATIONS

Figure 21.

(A–B) Variations in Vp/Vs and Vp of different mantle lithologies from (A) temperature and (B) depletion of iron, holding other factors constant. (C) Effect of increasing garnet on eclogitic (arclogitic) lower-crustal lithologies. Letters are keyed to Table 3. Variations were derived for 3 GPa, except the spinel peridotite (E) in panel A, which was calculated at 1.5 GPa. (D) All variations derived from Abers and Hacker (2016), applying the relationships of A–C as basis functions from the mean mantle (red star), assuming no variation in melt. Additional melt in the Isabella anomaly area moves up and left from the red star; less melt than average mantle in the region would move the reference state down toward the blue star. The deviations of the three different portions of the Isabella anomaly from three different inversions are indicated by the circles, which are colored by inversion and numbered by the subanomaly of the Isabella body that they represent.

Figure 21.

(A–B) Variations in Vp/Vs and Vp of different mantle lithologies from (A) temperature and (B) depletion of iron, holding other factors constant. (C) Effect of increasing garnet on eclogitic (arclogitic) lower-crustal lithologies. Letters are keyed to Table 3. Variations were derived for 3 GPa, except the spinel peridotite (E) in panel A, which was calculated at 1.5 GPa. (D) All variations derived from Abers and Hacker (2016), applying the relationships of A–C as basis functions from the mean mantle (red star), assuming no variation in melt. Additional melt in the Isabella anomaly area moves up and left from the red star; less melt than average mantle in the region would move the reference state down toward the blue star. The deviations of the three different portions of the Isabella anomaly from three different inversions are indicated by the circles, which are colored by inversion and numbered by the subanomaly of the Isabella body that they represent.

These results defined three trends relative to melt-free but otherwise typical asthenosphere. A temperature drop of 200 °C should increase Vp by ∼1.5% while decreasing Vp/Vs by ∼0.25%. Adding garnet up to the greatest values found in xenoliths will increase Vp by up to ∼5.5% while increasing Vp/Vs by ∼0.5%, though these values depend on other characteristics of the mineralogy. Depleting fertile mantle through progressive partial melt extraction will tend to increase Vp up to ∼2.5% while decreasing Vp/Vs by ∼1.3%. Thus, we investigated the combination of these variations that would reproduce the features found in the Isabella anomaly.

To use this information to interpret our tomographic results, we recognize that tomography constrains lateral variations in seismic properties and is fairly insensitive to the mean 1-D structure of Earth. Thus, to use those variations with temperature, iron content, and garnet content, we needed to know the reference state for our tomography. Although the area immediately around the Isabella anomaly is almost certainly melt-free (Park, 2004), the broader footprint of our study is almost certainly not. The presence of 1% melt fraction would reduce P-wave seismic velocities by 3.6% and increase Vp/Vs by 4.5% (Hammond and Humphreys, 2000) for a typical melt geometry. Scaling these values down for a presumed average melt fraction of 0.25% would mean that the tomographic reference mantle would be 0.9% slower than a melt-free reference and have a Vp/Vs ratio 1.2% higher than a melt-free reference.

We summarized this information as vectors originating from one possible melt-free reference point with the peak Vp/Vs perturbations and corresponding P-wave velocity perturbations for the three Isabella Vp/Vs anomalies using all starting models (labeled as 1, 2, 3), and these results are plotted in Figures 21 and 22I. We also estimated the values for melt-free mantle, assuming the average ambient mantle across our region contains 0.25% melt (Fig. 23).

Figure 22.

Interactive figure showing approximate seismic values produced from user-entered arbitrary values of regional melt percent and local temperature, depletion, and garnet differences (Fig. 21). Melt percent is allowed to range from −0.2% to 0.7%, local temperatures range from –200 °C to +900 °C, Mg# ranges from 84 to 100, and % garnet ranges from 0 to 100%. Peak anomalies from the three Isabella subanomalies for each of the three 3-D starting models can be shown using the checkboxes. The reader is invited to try and identify a simple set of variations that would explain all three subanomalies. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.22i to interact with Figure 22.

Figure 22.

Interactive figure showing approximate seismic values produced from user-entered arbitrary values of regional melt percent and local temperature, depletion, and garnet differences (Fig. 21). Melt percent is allowed to range from −0.2% to 0.7%, local temperatures range from –200 °C to +900 °C, Mg# ranges from 84 to 100, and % garnet ranges from 0 to 100%. Peak anomalies from the three Isabella subanomalies for each of the three 3-D starting models can be shown using the checkboxes. The reader is invited to try and identify a simple set of variations that would explain all three subanomalies. If you are reading the full-text html version of the paper, please view the PDF of this paper or visit https://doi.org/10.1130/GES02093.22i to interact with Figure 22.

Figure 23.

Isabella anomaly interpretation. (A) Summary of results of this study. (B) Relation of the three Isabella Vp/Vs anomalies (IA1–IA3 = 1–3) to physical causes. (C) Interpretation for this study. Hz—harzburgite, GtPy—garnet pyroxenite, GtPe—garnet peridotite, MM—possible Monterey microplate source, SN—possible Sierra Nevada source.

Figure 23.

Isabella anomaly interpretation. (A) Summary of results of this study. (B) Relation of the three Isabella Vp/Vs anomalies (IA1–IA3 = 1–3) to physical causes. (C) Interpretation for this study. Hz—harzburgite, GtPy—garnet pyroxenite, GtPe—garnet peridotite, MM—possible Monterey microplate source, SN—possible Sierra Nevada source.

The goal of these calculations was to determine which compositional and/or physical variations could reproduce the Vp and Vp/Vs values that were observed in the Isabella anomaly. We also sought to determine the processes that would generate high P velocity and low Vp/Vs (IA1, IA3), and high P velocity and high Vp/Vs (IA2). Plausible processes would also have to make sense with our anisotropy observations. Since our tomographic inversions tended to underestimate the true magnitude of seismic anomalies, the compositions and processes that best match our observations represent minimal deviations from ambient mantle.

Elastic Calculations Compared to Inversion Results

The first observation from this compilation was that the differences between the three subbodies of the Isabella anomaly are not parallel to any of the individual parameters affecting seismic velocity. Put another way, you cannot choose, say, IA1 as a reference and then decrease the temperature, or extract some melt, or add garnet to get to IA2. The differences between different parts of the Isabella anomaly require more than one variation.

Even considering the total range of possible reference states (i.e., varying the average amount of melt present in the region at different depths), the three parts of the Isabella anomaly seem discordant. For instance, if the regionally averaged mantle contains 0.25% melt at shallow depths, then IA1 could represent (1) fully depleted mantle that is 100 °C colder than average; (2) slightly depleted mantle that is 300 °C colder than average, or (3) garnet peridotite (up to 14% garnet) with no temperature contrast against surrounding mantle containing 0.35% melt. For IA3, the deepest anomaly, colder temperatures cannot be used as an explanation. The only interpretation solely relying on a thermal contrast is if the average mantle at this depth had ∼0.5% melt, and that this part of the Isabella anomaly was actually warmer than average mantle by 100–200 °C, which is entirely inconsistent with the melt-free gradients illustrated. Instead, the anomaly can be explained by a combination of a partially molten surrounding mantle and an increase in Mg# (depletion) within the anomaly. If the surrounding mantle contains 0.35% melt, then IA3 only needs to be melt-free. However, if the surrounding mantle contains less melt (say 0.25%), then IA3 could be melt-free and depleted (Fig. 21). Clearly, a simple interpretation of the Isabella anomaly cannot be made in a vacuum. Instead, we consider the two main hypotheses for this body and their implications after summarizing the main variations capable of affecting our observables.

Physical Processes and Related Vp-Vp/Vs Trends

While thermal effects have often been considered to be the dominant influence on seismic velocity variations in the upper mantle (e.g., Goes et al., 2000), we ruled out temperature variations as the sole contributor to the variations in Vp/Vs and velocity observed within the western foothills of the southern Sierra Nevada. A decrease in temperature would consistently lead to an increase in P velocities and a decrease in Vp/Vs, as seen in Figure 21. The discrepancy between Vp/Vs and velocity within the Isabella anomaly cannot be explained only by variations in temperature. If the Isabella velocity anomaly were solely a thermal anomaly, it would result in covarying Vp and Vp/Vs, which is not observed.

There is still debate as to whether seismic velocities and Vp/Vs are sensitive to melt depletion in the upper mantle. Melt depletion, the result of partial melt extraction, has been observed to increase Vs more than Vp, decrease density, and decrease Vp/Vs with increasing Mg# [Mg/(Mg + Fe)] for peridotites (Lee, 2003). However, Schutt and Lesher (2006) suggested that melt depletion has minimal effect on velocity. Afonso et al. (2010) argued that any correlation between Vp/Vs and Mg# is limited to garnet-bearing assemblages with temperatures less than 900 °C and breaks down when spinel is the stable Al-rich phase. Melt depletion is often used to explain the neutral buoyancy required to stabilize cratonic lithosphere, but it has also been used to explain the assumed neutral buoyancy and high velocities of a proposed remnant Monterey microplate.

Dehydration has also been mentioned as a possible explanation for the high velocities observed within the Isabella anomaly, specifically, that the anomaly represents remnant oceanic lithosphere (i.e., Wang et al., 2013). For this interpretation, the slab is considered to have been dehydrated during the formation of oceanic lithosphere at a mid-ocean ridge. Dehydration of oceanic lithosphere would result in an increase in seismic velocities and a decrease in Vp/Vs (e.g., Faccenda, 2014). On a related note, a decrease in serpentinization would also result in an increase in seismic velocities and a decrease in Vp/Vs. Both of these would trend in approximately the same direction as changes in temperature and iron content.

Orthopyroxene enrichment has been invoked to explain the high shear-wave velocities and low Vp/Vs values observed in regions above subduction zones (e.g., Wagner et al., 2005, 2008). However, these results trend similarly to the other mechanisms. Of all the effects mentioned, only an increase in garnet will lead to both an increase in velocities and Vp/Vs. Worthington et al. (2013) observed that eclogite formed from high-pressure metamorphism generally exhibited higher Vp/Vs than peridotite. Results from this exercise are summarized in Table 4.

TABLE 4.

Vp AND Vp/Vs TRENDS

Isabella Anomaly Derived from the Sierra Nevada

The strongest case for the presence of sub-Sierran lithosphere within the Isabella anomaly is within IA2. High velocities, high Vp/Vs, and low anisotropy values characterize this feature. As mentioned above, the seismic velocity and Vp/Vs trend can be explained by an increase in garnet. The presence of garnet would also reduce or preclude the development of azimuthal anisotropy. However, we cannot rule out the possibility that the low amount of azimuthal anisotropy coincides with strong radial anisotropy. We suggest that IA2 represents the mafic “eclogitic” root of the Sierra Nevada, sometimes termed “arclogite,” due to its origin under an arc and differences from eclogites produced in different environments.

Because IA3 exhibits higher anisotropy than IA2, we posit that IA3 is peridotitic and could be derived from beneath the Sierra Nevada. IA3 could represent garnet peridotite, if the surrounding mantle contained ∼0.35% partial melt. Geometrically, it is possible that IA3 represents Sierran garnet peridotitic upper mantle. IA3 is deeper and to the east of IA2. A conceptual stratigraphic column through the Sierra Nevada prior to foundering would place garnet peridotite somewhere beneath the garnet pyroxenitic root. The plunge of what originally resembled a vertical section most likely reflects the process of foundering, perhaps due to asymmetry in initial instabilities (e.g., Saleeby et al., 2013), tilting during downwelling (Saleeby et al., 2003) or convective motion in the mantle (e.g., Zandt, 2003).

We also suggest that IA1 is peridotitic given its lower Vp/Vs values than IA2, which are of similar magnitude to IA3. The seismic anisotropy for this feature is unclear, because smearing effects might contribute to the lower anisotropy produced by our inversions. Consequently, there is less confidence that this feature is related to the Sierra Nevada. Other relevant interpretations for IA1 include Great Valley ophiolite, in situ oceanic upper mantle from a Jurassic oblique rifting event, Early Cretaceous arc material, or remnant oceanic material from the Monterey microplate.

Isabella Anomaly Derived from Oceanic Lithosphere (Monterey Microplate)

Using Abers and Hacker’s (2016) computational tools, we investigated the plausibility that a portion of the Isabella anomaly is derived from remnant oceanic lithosphere of the Monterey microplate. This interpretation of the Isabella anomaly attributes the high velocities to a combination of cooler temperatures (than surrounding mantle), depletion, and dehydration (i.e., Pikser et al., 2012; Wang et al., 2013; Cox et al., 2016). We did not include eclogite in our oceanic compositions. While eclogite is produced in oceanic lithosphere as basaltic crust descends below ∼100 km depth, the density increase associated with this transition would make the slab negatively buoyant (i.e., Pikser et al., 2012) and prone to foundering. The temperature contrast between the slab and asthenospheric (ambient) mantle ranges only from 100° to 300 °C, as calculated from thermal models (Van Wijk et al., 2001). The slab cannot be too much colder than surrounding mantle, owing to its young age at the time of subduction and ensuing thermal reequilibration with surrounding asthenosphere (e.g., Bohannon and Parsons, 1995).

The fertile (enriched?) and depleted harzburgite calculations, for temperature contrasts of 100 °C to 300 °C, define the seismic perturbations that can be explained by the presence of an oceanic slab (Fig. 22). Of the three Isabella anomalies, only IA1 displays perturbations consistent with a slab origin. While there was considerable variability in seismic values for IA1 between the starting models (as mentioned in the previous section), results from two of our starting models plotted within this “slab anomaly” region. Thus, it is seismically possible that depleted, slightly cooler oceanic material can produce a sufficiently high-velocity seismic anomaly to reproduce observations of IA1. The most plausible location for remnant oceanic material within the Isabella anomaly is in IA1, the shallowest and westernmost feature.

It is also possible that IA3 is remnant oceanic material. If this were the case, IA3 would represent an older portion of the slab lacking a cold thermal anomaly to explain the Vp and Vp/Vs values. Accordingly, if this were older Monterey material, it would represent material that has thermally reequilibrated with surrounding mantle. However, because a fossil slab should have fairly uniform compositional and seismic properties, it is difficult to reconcile IA2’s position between IA1 and IA3 as part of the Monterey microplate, complicating an interpretation of IA3 as a fossil slab. This is especially true if the slab is intact as suggested by Pikser et al. (2012).

Isabella Anomaly as Both Sierran and Monterey Lithosphere?

It is seismically possible that the Isabella Vp/Vs anomaly is both Sierran and Monterey lithosphere. For this interpretation, IA1 is Monterey upper mantle and IA2 + IA3 are Sierran lower crust and upper mantle (Fig. 23).

It is also possible that both features are currently present beneath the southwestern foothills of the Sierra Nevada. This requires either happenstance as the Monterey plate arrived as Sierran lithosphere descended to its east, or some interference where either the arrival of the Monterey slab attracted the Sierran downwelling, or the presence of the downwelling captured the slab and slowed or stopped its northward trajectory. As it is unlikely that the Sierran drip would migrate northward at rates comparable to the motion of the Pacific plate, given its position near the southern edge of the Sierra, it would seem that this juxtaposition would have to have occurred fairly recently as the offshore remnant of the Monterey plate only came opposite to the Isabella anomaly ca. 3 Ma, and interaction at that time would require the subducted part of the plate to lie to the northeast of the offshore remnant instead of to the east-southeast, as at present. The removal of Sierran lithosphere could have begun by 12 Ma, but it seems certain to have been largely complete by 3.5 Ma, when unusual potassic magmas erupted across much of the southern Sierra (Farmer et al., 2002; Manley et al., 2000) and the first mantle xenoliths containing no evidence of older deep lithosphere appear (Ducea and Saleeby, 1996, 1998). Similarly, the lithospheric downwelling might have appeared as recently as 3–2 Ma (Le Pourhiet et al., 2006) or before 5 Ma (Saleeby et al., 2013). Short of other evidence or a set of geodynamic experiments to explore how two such different bodies might interact, we cannot rule out such scenarios, but they seem to require special pleading.

Partial Melt at Upper-Mantle Depths beneath the Mono-Inyo Region?

While results from this study image high Vp/Vs within the Clear Lake, Coso, and Mono-Inyo volcanic fields, we limit our discussion to the Mono anomaly, since Clear Lake and Coso were more poorly resolved. Given considerable Pleistocene and Recent volcanism in the area (e.g., Hildreth, 2004), the presence of both elevated temperatures and partial melt in the upper mantle of the Mono-Inyo region is likely. Anomalies in Vp and Vp/Vs within the Mono anomaly are roughly −2% and +3.5%, respectively.

We first consider a purely thermal origin for these anomalies. Using the temperature derivatives from Karato (1993), Schmandt and Humphreys (2010a) calculated that a 500 K temperature contrast would produce Vp and Vp/Vs anomalies of −5% and +3%, respectively. While their Vp/Vs variation is comparable to our results, our Vp anomaly is too small. If we instead use a smaller temperature contrast of 300 K, Qs of 100, and Qp of 200, the Vp and Vp/Vs variation becomes −3% and +1.7%, respectively, and we fail to match the Vp/Vs anomaly. Thus, to match the Vp variations, our temperature contrast needs to be smaller, which will not fully account for our Vp/Vs variations.

Partial melt appears to be necessary to reconcile our results with relationships described in the literature. Hammond and Humphreys (2000) found that 1% partial melt should cause Vp to decline by no less than 3.9% and Vs to decline by no less than 7.9%. The minimum amount of partial melt that could explain the Mono anomaly (in the absence of anelastic and thermal effects) ranges from 0.56% and 0.78%. A melt anomaly alone fits our average observations fairly well, but the difference in shape of the Vp and Vp/Vs anomalies (Fig. 15) suggests a more complex variation of melt and temperature in this magmatic system.

Our interpretation is in agreement with previous geophysical studies beneath the Sierra Nevada–eastern California and southern California regions. Inversion of a magnetotelluric profile across Yosemite National Park, and extending as far east as Mono Lake, produced a resistive upper-mantle feature beneath the Mono Lake region confined to depths of ∼40–100 km (Ostos and Park, 2012). However, sensitivity tests revealed that this feature does not need to be resistive and could be conductive. If conductive, then this feature might be related to conductive anomalies found beneath the Sierra Nevada upper mantle that have been attributed to upwelling asthenosphere containing <1% partial melt (e.g., Ostos and Park, 2012; Park, 2004). Yang and Forsyth (2008) also suggested that the low shear velocities and related attenuation found in the upper mantle of southern California could be attributed to damp asthenosphere containing partial melt. For Schmandt and Humphreys (2010b), the presence of <1% partial melt in the asthenosphere extending from 150 km to 200 km depth could explain the magnitudes in their Vp, Vs, and Vp/Vs perturbations.

Anisotropy Tomography Limitations

Our anisotropy tomography method represents a simple approach for providing depth constraints on anisotropic upper-mantle features. However, this technique has its limitations. First, this method assumes that delays or advances of S waves only depend on the anisotropy of the medium and the length of the raypath in that medium. An assumption inherent in our approach is that the anisotropic fabric exhibits hexagonal symmetry with a horizontal symmetry axis (azimuthal anisotropy). Shear-wave splitting patterns for azimuthally anisotropic material depend on both the azimuth relative to the fast orientation and the angle of incidence (e.g., Savage, 1999). For example, Schulte-Pelkum and Blackman (2003) observed that the fast direction can vary by up to 25° from the fast symmetry axis for S waves at shallow incidence angles. The splitting delay time can also vary with increasing incidence angle. However, because we are not measuring split times, this gets more complicated. Using S and SKS phases with varying incidence angles will plausibly average out the effects of anisotropy strength.

The use of steeply arriving teleseismic P, S, and SKS waves makes our approach less sensitive to plunging anisotropy or radial anisotropy. Having said that, incorporating teleseismic S with varying incidence angles (all of which are within the shear-wave window of Savage, 1999) might introduce slight sensitivity to radial anisotropy.

Second, by rotating the horizontal components to one fast and slow orientation, we are restricted from interpreting results in regions where SKS fast orientations deviate from those of the Sierra Nevada. Third, by using shear-wave fast and slow traveltime residuals to constrain anisotropy, we assume any anisotropy originates from anisotropic layers with a single orientation, thus ignoring effects from multiple layers of anisotropy with different orientations. Also, the presence of multiple layers would be difficult to resolve using S and SKS phases, since our tomograms do exhibit vertical smearing. Therefore, our tomography represents apparent azimuthal anisotropy in these cases. These points are obviously an issue within the vicinity of the San Andreas fault in central California. Here, not only do the fast orientations for stations in this region deviate from those in the Sierra Nevada, but back-azimuthal variations in SKS splitting parameters are best explained by a two-layer anisotropic model (Savage and Silver, 1993; Özalaybey and Savage, 1995). SKS splitting measurements within the Sierra Nevada exhibit less variation with back azimuth (this study; Özalaybey and Savage, 1995; Polet and Kanamori, 2002), making our interpretations and focus on this range less problematic.

Finally, by not accounting for anisotropy in our P-wave traveltimes, the high-velocity Isabella anomaly (should it represent vertically deformed fabric, which our teleseismic phases cannot resolve) could partly be explained by anisotropic artifacts mapping into velocity. Although there is a correlation between the prominent high-velocity upper-mantle anomalies (Isabella and Redding) and low anisotropy, this coincidence is imperfect. The low anisotropy found in the northern Sierra Nevada is associated with ambient P-wave velocity.

Anisotropy within the Sierra Nevada

The high anisotropy observed in the central Sierra crest and Great Valley–western foothills (Figs. 10 and 18) is consistent with other studies suggesting the presence of anisotropic fabrics in this region (Özalaybey and Savage, 1995; Polet and Kanamori, 2002; Zandt et al., 2004; Frassetto et al., 2011). The foothills feature (38°N) found at 120 km depth is adjacent to station CMB (Fig. 19), which has been characterized by high SKS splitting delays of up to 2 s (Özalaybey and Savage, 1995; Polet and Kanamori, 2002). Özalaybey and Savage (1995) suggested that the large delays indicate that this anisotropic feature resides in the asthenosphere (which is consistent with our peak depth of 120 km) and possibly results from asthenospheric infill related to the slabless window caused by the passage of the Farallon plate.

In contrast, the central Sierra Nevada anisotropic feature to the east is much shallower at 20–40 km depth. Receiver function studies in the southern Sierra Nevada observed mid- to lower-crust negative polarity arrivals around ∼20–30 km, indicative of a possible sheared/deformed layer representing a remnant detachment zone related to the foundering of the Sierran mafic root (Zandt et al., 2004) or top-to-the-west low-angle normal faulting associated with extension to the east (Jones and Phinney, 1998). The zone observed in this study could be related. Below this anisotropic layer, there is a zone of low anisotropy minimized from 170 to 220 km depth. This anisotropic low is associated with slightly higher Vp/Vs values than surrounding mantle, possibly reflecting the presence of some additional melt. When considering the overall pattern of anisotropy, the combination of low-azimuthal anisotropy lying under a thinner high-anisotropy zone could reflect vertical upwelling of asthenosphere that drives a horizontal component of mantle flow at shallower depths (Fig. 19). Shear on subhorizontal planes would tend to align minerals to produce the observed azimuthal anisotropy.

CONCLUSIONS

Vp/Vs anisotropy tomography of the Sierra Nevada presents new insights on the compositional, physical, and deformational structure of the upper mantle in this region. Results from this study present an unexpected image of a compound Isabella anomaly. Instead of being one unimodal feature, as previously identified in P-wave tomographic studies, the Isabella anomaly is composed of three distinct bodies differing in Vp/Vs and anisotropy values. Synthetic tests indicate that the data set used here should be able to resolve these separate anomalies and that they are not artifacts. The core of the P-velocity anomaly (IA2) is characterized by fast velocities, regionally average Vp/Vs amplitudes, and low anisotropy. The westernmost (and shallowest) anomaly (IA1) has lower Vp/Vs and ambiguous anisotropy. The easternmost (and deepest) anomaly (IA3) also contains lower Vp/Vs amplitudes but the highest anisotropy of the three. The discordant nature of the Vp/Vs anisotropy anomaly suggests that these features cannot be explained by one composition or process alone. Rock physics calculations indicate that the Isabella anomaly reflects two to three separate mafic to ultramafic lithologies. Anomaly IA2 is most likely garnet-rich “eclogitic” material derived from the Sierra Nevada. The anomalies with lower Vp/Vs amplitudes (IA1 and IA3) point to a more peridotitic and depleted material that is less diagnostic of the feature’s origins. A temperature contrast is not necessary to explain IA3, but it may contribute to the high velocities and low Vp/Vs observed in IA1. Based on these calculations, the seismic constraints presented here allow for the possibility that portions of the Isabella anomaly that appear to be peridotitic could be either Sierran or Monterey in origin.

The largest Vp/Vs anomalies in our study area, which are associated with slow velocities, correlate to regions of known magmatism. The largest of these features is the Mono anomaly, located within the eastern Sierra Nevada, which includes the Long Valley caldera, Mammoth Mountain, the Mono-Inyo Craters, and Mono Lake. Low anisotropy is observed within the Mono anomaly, possibly due to upwelling of partially molten asthenospheric material.

Anisotropy tomography reveals patterns that are broadly consistent with SKS splitting delay measurements and previous receiver function studies within the Sierra Nevada region. High anisotropy is localized to the central Sierran crest and Great Valley–western Sierra foothills. The central crest feature is found between 20 and 40 km depth, a depth range that is consistent with shear fabrics inferred from receiver function studies. Below this feature, there lies a zone of lower anisotropy between 170 and 220 km. Together, these results suggests that upwelling asthenospheric material in the upper mantle might be causing horizontal flow fabrics at shallower depths. The Great Valley–western foothills feature is much deeper at 120 km depth and could represent mantle flowing in as the Gorda slab moved to the north.

ACKNOWLEDGMENTS

Steve Roecker provided and helped suggest the changes we made to his sphfd code. G. Pavlis helped guide us to modify his dbxcor to make our unconventional S-wave picks. Vera Schulte-Pelkum aided our discussion on anisotropy. Finally, Don Forsyth and an anonymous reviewer provided insightful critiques of our approach and presentation that helped to clarify our study. This project was funded from the National Science Foundation (NSF) grant no. DGE 1144083, with some additional work supported by NSF EAR-1547123. The original data acquisition and some of the early work on this project were supported by NSF EAR grants 0454535, 0454524, and 0454554 to the University of Colorado, University of South Carolina, and University of Arizona, respectively.

1Supplemental Materials. Details of SKS measurements, two additional synthetic tests, and tables of station, event, and arrival data. Please visit https://doi.org/10.1130/GES02093.S1 or access the full-text article on www.gsapubs.org to view the Supplemental Materials.
Science Editor: Shanaka de Silva

REFERENCES CITED

Abers
,
G.A.
, and
Hacker
,
B.R.
,
2016
,
A MATLAB toolbox and Excel workbook for calculating the densities, seismic wave speeds, and major element composition of minerals and rocks at pressure and temperature
:
Geochemistry Geophysics Geosystems
 , v.
17
, no.
2
, p.
616
624
, https://doi.org/10.1002/2015GC006171.
Afonso
,
J.C.
,
Ranalli
,
G.
,
Fernàndez
,
M.
,
Griffin
,
W.L.
,
O’Reilly
,
S.Y.
, and
Faul
,
U.
,
2010
,
On the Vp/Vs–Mg# correlation in mantle peridotites: Implications for the identification of thermal and compositional anomalies in the upper mantle
:
Earth and Planetary Science Letters
 , v.
289
, p.
606
618
, https://doi.org/10.1016/j.epsl.2009.12.005.
Benz
,
H.M.
, and
Zandt
,
G.
,
1993
,
Teleseismic tomography: Lithospheric structure of the San Andreas fault system in northern and central California
, in
Iyer
,
H.M.
, and
Hirahara
,
K.
, eds.,
Seismic Tomography: Theory and Practice
 :
New York
,
Chapman and Hall
, p.
440
465
.
Benz
,
H.M.
,
Zandt
,
G.
, and
Oppenheimer
,
D.H.
,
1992
,
Lithospheric structure of Northern California from teleseismic images of the upper mantle
:
Journal of Geophysical Research
 , v.
97
, p.
4791
4807
, https://doi.org/10.1029/92JB00067.
Bezada
,
M.J.
,
Faccenda
,
M.
, and
Toomey
,
D.R.
,
2016
,
Representing anisotropic subduction zones with isotropic velocity models: A characterization of the problem and some steps on a possible path forward
:
Geochemistry Geophysics Geosystems
 , v.
17
, p.
3164
3189
, https://doi.org/10.1002/2016GC006507.
Biasi
,
G.P.
, and
Humphreys
,
E.D.
,
1992
,
P-wave image of the upper mantle structure of central California and southern Nevada
:
Geophysical Research Letters
 , v.
19
, no.
11
, p.
1161
1164
, https://doi.org/10.1029/92GL00439.
Bohannon
,
R.G.
, and
Parsons
,
T.
,
1995
,
Tectonic implications of post–30 Ma Pacific and North American relative plate motions
:
Geological Society of America Bulletin
 , v.
107
, p.
937
959
, https://doi.org/10.1130/0016-7606(1995)107<0937:TIOPMP>2.3.CO;2.
Boyd
,
O.S.
,
Jones
,
C.H.
, and
Sheehan
,
A.F.
,
2004
,
Foundering lithosphere imaged beneath the southern Sierra Nevada, California, USA
:
Science
 , v.
305
, p.
660
662
, https://doi.org/10.1126/science.1099181.
Brocher
,
T.M.
,
2005
,
Empirical relations between elastic wavespeeds and density in the Earth’s crust
:
Bulletin of the Seismological Society of America
 , v.
95
, p.
2081
2092
, https://doi.org/10.1785/0120050077.
Cox
,
P.
,
Stubailo
,
I.
, and
Davis
,
P.
,
2016
,
Receiver function and geometric tomography along the Monterey microplate to test slab delamination or lithospheric drip models of the Isabella anomaly, California
:
Bulletin of the Seismological Society of America
, v.
106
, p.
267
280
, https://doi.org/10.1785/0120140339.
Ducea
,
M.N.
, and
Saleeby
,
J.B.
,
1996
,
Buoyancy sources for a large, unrooted mountain range, the Sierra Nevada, California: Evidence from xenolith thermobarometry
:
Journal of Geophysical Research
 , v.
101
, p.
8229
8244
, https://doi.org/10.1029/95JB03452.
Ducea
,
M.N.
, and
Saleeby
,
J.B.
,
1998
,
A case for delamination of the deep batholithic crust beneath the Sierra Nevada, California
:
International Geology Review
 , v.
40
, p.
78
93
, https://doi.org/10.1080/00206819809465199.
Faccenda
,
M.
,
2014
,
Water in the slab: A trilogy
:
Tectonophysics
 , v.
614
, p.
1
30
, https://doi.org/10.1016/j.tecto.2013.12.020.
Farmer
,
G.L.
,
Glazner
,
A.F.
, and
Manley
,
C.R.
,
2002
,
Did lithospheric delamination trigger late Cenozoic potassic volcanism in the southern Sierra Nevada, California?
:
Geological Society of America Bulletin
 , v.
114
, p.
754
768
, https://doi.org/10.1130/0016-7606(2002)114<0754:DLDTLC>2.0.CO;2.
Fliedner
,
M.M.
,
Ruppert
,
S.
, and
SSCD Working Group
,
1996
,
Three-dimensional crustal structure of the southern Sierra Nevada from seismic fan profiles and gravity modeling
:
Geology
 , v.
24
, p.
367
370
, https://doi.org/10.1130/0091-7613(1996)024<0367:TDCSOT>2.3.CO;2.
Fountain
,
D.M.
, and
Christensen
,
N.I.
,
1989
,
Composition of the continental crust and upper mantle: A review
, in
Pakiser
,
L.C.
, and
Mooney
,
W.D.
, eds.,
Geophysical Framework of the Continental United States
 :
Geological Society of America Memoir
172
, p.
711
742
, https://doi.org/10.1130/MEM172-p711.
Frassetto
,
A.
,
Zandt
,
G.
,
Gilbert
,
H.
,
Owens
,
T.J.
, and
Jones
,
C.
,
2011
,
Lithospheric structure of the Sierra Nevada from receiver functions and implications for lithospheric foundering
:
Geosphere
 , v.
7
, no.
4
, p.
898
921
, https://doi.org/10.1130/GES00570.1.
Gilbert
,
H.
,
Yang
,
Y.
,
Forsyth
,
D.W.
,
Jones
,
C.H.
,
Owens
,
T.J.
,
Zandt
,
G.
, and
Stachnik
,
J.C.
,
2012
,
Imaging lithospheric foundering in the structure of the Sierra Nevada
:
Geosphere
 , v.
8
, p.
1310
1330
, https://doi.org/10.1130/GES00790.1.
Goes
,
S.
,
Govers
,
R.
, and
Vacher
,
P.
,
2000
,
Shallow upper mantle temperatures under Europe from P and S wave tomography
:
Journal of Geophysical Research
 , v.
105
, p.
11,153
11,169
, https://doi.org/10.1029/1999JB900300.
Hacker
,
B.R.
,
Abers
,
G.A.
, and
Peacock
,
S.M.
,
2003
,
Subduction factory 1: Theoretical mineralogy, density, seismic wavespeeds, and H2O content
:
Journal of Geophysical Research
 , v.
108
,
2029
, https://doi.org/10.1029/2001JB001127.
Hammond
,
W.C.
, and
Humphreys
,
E.D.
,
2000
,
Upper mantle seismic wave velocity: Effects of realistic partial melt geometries
:
Journal of Geophysical Research–Solid Earth
 , v.
105
, p.
10,975
10,986
, https://doi.org/10.1029/2000JB900041.
Hildreth
,
W.
,
2004
,
Volcanological perspectives on Long Valley, Mammoth Mountain, and Mono Craters: Several contiguous but discrete systems
:
Journal of Volcanology and Geothermal Research
 , v.
136
, p.
169
198
, https://doi.org/10.1016/j.jvolgeores.2004.05.019.
Hole
,
J.A.
, and
Zelt
,
B.C.
,
1995
,
3-D finite-difference reflection traveltimes
:
Geophysical Journal International
 , v.
121
, p.
427
434
, https://doi.org/10.1111/j.1365-246X.1995.tb05723.x.
Hongsresawat
,
S.
,
Panning
,
M.P.
,
Russo
,
R.M.
,
Foster
,
D.A.
,
Monteiller
,
V.
, and
Chevrot
,
S.
,
2015
,
USArray shear wave splitting shows seismic anisotropy from both lithosphere and asthenosphere
:
Geology
 , v.
43
, p.
667
670
, https://doi.org/10.1130/G36610.1.
Jones
,
C.H.
, and
Phinney
,
R.
,
1997
,
Continental Mountains in Extensional Environments: The Sierran Paradox
:
International Federation of Digital Seismograph Networks, Other/Seismic Network
, https://doi.org/10.7914/SN/XJ_1997.
Jones
,
C.H.
, and
Phinney
,
R.A.
,
1998
,
Seismic structure of the lithosphere from teleseismic converted arrivals observed at small arrays in the southern Sierra Nevada and vicinity, California
:
Journal of Geophysical Research
 , v.
103
, p.
10,065
10,090
, https://doi.org/10.1029/97JB03540.
Jones
,
C.H.
,
Kanamori
,
H.
, and
Roecker
,
S.W.
,
1994
,
Missing roots and mantle “drips”: Regional Pn and teleseismic arrival times in the southern Sierra Nevada and vicinity, California
:
Journal of Geophysical Research
 , v.
99
, p.
4567
4601
, https://doi.org/10.1029/93JB01232.
Jones
,
C.H.
,
Reeg
,
H.
,
Zandt
,
G.
,
Gilbert
,
H.
,
Owens
,
T.J.
, and
Stachnik
,
J.
,
2014
,
P-wave tomography of potential convective downwellings and their source regions, Sierra Nevada, California
:
Geosphere
 , v.
10
, p.
505
533
, https://doi.org/10.1130/GES00961.1.
Karato
,
S.
,
1993
,
Importance of anelasticity in the interpretation of seismic tomography
:
Geophysical Research Letters
 , v.
20
, p.
1623
1626
, https://doi.org/10.1029/93GL01767.
Kennett
,
B.L.N.
, and
Engdahl
,
E.R.
,
1991
,
Traveltimes for global earthquake location and phase identification
:
Geophysical Journal International
 , v.
105
, p.
429
465
, https://doi.org/10.1111/j.1365-246X.1991.tb06724.x.
Lee
,
C.-T.A.
,
2003
,
Compositional variation of density and seismic velocities in natural peridotites at STP conditions: Implications for seismic imaging of compositional heterogeneities in the upper mantle
:
Journal of Geophysical Research–Solid Earth
 , v.
108
,
2441
, https://doi.org/10.1029/2003JB002413.
Lee
,
C.-T.A.
,
Cheng
,
X.
, and
Horodyskyj
,
U.
,
2006
,
The development and refinement of continental arcs by primary basaltic magmatism, garnet pyroxenite accumulation, basaltic recharge and delamination: Insights from the Sierra Nevada, California
:
Contributions to Mineralogy and Petrology
 , v.
151
, p.
222
242
, https://doi.org/10.1007/s00410-005-0056-1.
Le Pourhiet
,
L.
,
Gurnis
,
M.
, and
Saleeby
,
J.
,
2006
,
Mantle instability beneath the Sierra Nevada mountains in California and Death Valley extension
:
Earth and Planetary Science Letters
 , v.
251
, no.
1–2
, p.
104
119
, https://doi.org/10.1016/j.epsl.2006.08.028.
Liu
,
K.H.
,
Elsheikh
,
A.
,
Lemnifi
,
A.
,
Purevsuren
,
U.
,
Ray
,
M.
,
Refayee
,
H.
,
Yang
,
B.
,
Yu
,
Y.
, and
Gao
,
S.S.
,
2014
,
A uniform database of teleseismic shear wave splitting measurements for the western and central United States
:
Geochemistry Geophysics Geosystems
 , v.
15
, p.
2075
2085
, https://doi.org/10.1002/2014GC005267.
Manley
,
C.R.
,
Glazner
,
A.F.
, and
Farmer
,
G.L.
,
2000
,
Timing of volcanism in the Sierra Nevada of California: Evidence for Pliocene delamination of the batholithic root?
:
Geology
 , v.
28
, no.
9
, p.
811
814
, https://doi.org/10.1130/0091-7613(2000)28<811:TOVITS>2.0.CO;2.
Moschetti
,
M.P.
,
Ritzwoller
,
M.H.
,
Lin
,
F.C.
, and
Yang
,
Y.
,
2010
,
Crustal shear wave velocity structure of the western United States inferred from ambient seismic noise and earthquake data
:
Journal of Geophysical Research–Solid Earth
 , v.
115
,
B10306
, https://doi.org/10.1029/2010JB007448.
Ostos
,
L.
, and
Park
,
S.K.
,
2012
,
Foundering lithosphere imaged with magnetotelluric data beneath Yosemite National Park, California
:
Geosphere
 , v.
8
, p.
98
104
, https://doi.org/10.1130/GES00657.1.
Owens
,
T.
,
Zandt
,
G.
, and
Jones
,
C.
,
2005
,
Sierra Nevada EarthScope Project
:
International Federation of Digital Seismograph Networks, Other/Seismic Network
, https://doi.org/10.7914/SN/XE_2005.
Özalaybey
,
S.
, and
Savage
,
M.K.
,
1995
,
Shear-wave splitting beneath western United States in relation to plate tectonics
:
Journal of Geophysical Research–Solid Earth
 , v.
100
, p.
18,135
18,149
, https://doi.org/10.1029/95JB00715.
Paige
,
C.C.
, and
Saunders
,
M.A.
,
1982
,
LSQR: An algorithm for sparse linear equations and sparse least squares
:
ACM Transactions on Mathematical Software
 , v.
8
, p.
43
71
, https://doi.org/10.1145/355984.355989.
Park
,
S.K.
,
2004
,
Mantle heterogeneity beneath eastern California from magnetotelluric measurements
:
Journal of Geophysical Research–Solid Earth
 , v.
109
,
B09406
, https://doi.org/10.1029/2003jb002948.
Pavlis
,
G.L.
, and
Vernon
,
F.L.
,
2010
,
Array processing of teleseismic body waves with the USArray
:
Computers & Geosciences
 , v.
36
, p.
910
920
, https://doi.org/10.1016/j.cageo.2009.10.008.
Pikser
,
J.E.
,
Forsyth
,
D.W.
, and
Hirth
,
G.
,
2012
,
Along-strike translation of a fossil slab
:
Earth and Planetary Science Letters
 , v.
331
, p.
315
321
, https://doi.org/10.1016/j.epsl.2012.03.027.
Polet
,
J.
, and
Kanamori
,
H.
,
2002
,
Anisotropy beneath California: Shear wave splitting measurements using a dense broadband array
:
Geophysical Journal International
 , v.
149
, p.
313
327
, https://doi.org/10.1046/j.1365-246X.2002.01630.x.
Raikes
,
S.A.
,
1980
,
Regional variations in upper mantle structure beneath Southern California
:
Geophysical Journal International
 , v.
63
, p.
187
216
, https://doi.org/10.1111/j.1365-246X.1980.tb02616.x.
Roecker
,
S.
,
Thurber
,
C.
,
Roberts
,
K.
, and
Powell
,
L.
,
2006
,
Refining the image of the San Andreas fault near Parkfield, California, using a finite difference travel time computation technique
:
Tectonophysics
 , v.
426
, p.
189
205
, https://doi.org/10.1016/j.tecto.2006.02.026.
Saleeby
,
J.
,
Ducea
,
M.
, and
Clemens-Knott
,
D.
,
2003
,
Production and loss of high-density batholithic root, southern Sierra Nevada, California
:
Tectonics
 , v.
22
, p.
1604
, https://doi.org/10.1029/2002TC001374.
Saleeby
,
J.B.
,
Saleeby
,
Z.
, and
Le Pourhiet
,
L.
,
2013
,
Epeirogenic transients related to mantle lithosphere removal in the southern Sierra Nevada region, California: Part II. Implications of rock uplift and basin subsidence relations
:
Geosphere
 , v.
9
, no.
3
, p.
394
425
, https://doi.org/10.1130/GES00816.1.
Savage
,
M.K.
,
1999
,
Seismic anisotropy and mantle deformation: What have we learned from shear wave splitting?
:
Reviews of Geophysics
 , v.
37
, no.
1
, p.
65
106
, https://doi.org/10.1029/98RG02075.
Savage
,
M.K.
, and
Silver
,
P.G.
,
1993
,
Mantle deformation and tectonics: Constraints from seismic anisotropy in western United States
:
Physics of the Earth and Planetary Interiors
 , v.
78
, p.
207
227
, https://doi.org/10.1016/0031-9201(93)90156-4.
Schmandt
,
B.
, and
Humphreys
,
E.
,
2010a
,
Complex subduction and small-scale convection revealed by body-wave tomography of the western United States upper mantle
:
Earth and Planetary Science Letters
 , v.
297
, p.
435
445
, https://doi.org/10.1016/j.epsl.2010.06.047.
Schmandt
,
B.
, and
Humphreys
,
E.
,
2010b
,
Seismic heterogeneity and small-scale convection in the southern California upper mantle
:
Geochemistry Geophysics Geosystems
 , v.
11
,
Q05004
, https://doi.org/10.1029/2010GC003042.
Schulte-Pelkum
,
V.
, and
Blackman
,
D.K.
,
2003
,
A synthesis of seismic P and S anisotropy
:
Geophysical Journal International
 , v.
154
, p.
166
178
, https://doi.org/10.1046/j.1365-246X.2003.01951.x.
Schutt
,
D.L.
, and
Humphreys
,
E.D.
,
2004
,
P and S wave velocity and VP/VS in the wake of the Yellowstone hot spot
:
Journal of Geophysical Research–Solid Earth
 , v.
109
, no. B1, B01305, https://doi.org/10.1029/2003JB002442.
Schutt
,
D.L.
, and
Lesher
,
C.E.
,
2006
,
Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite
:
Journal of Geophysical Research–Solid Earth
 , v.
111
,
B05401
, https://doi.org/10.1029/2003JB002950.
Shen
,
W.
,
Ritzwoller
,
M.H.
, and
Schulte-Pelkum
,
V.
,
2013
,
A 3-D model of the crust and uppermost mantle beneath the central and western US by joint inversion of receiver functions and surface wave dispersion
:
Journal of Geophysical Research–Solid Earth
 , v.
118
, p.
262
276
, https://doi.org/10.1029/2012JB009602.
Silver
,
P.G.
, and
Chan
,
W.W.
,
1991
,
Shear wave splitting and subcontinental mantle deformation
:
Journal of Geophysical Research
 , v.
96
, p.
16,429
16,454
, https://doi.org/10.1029/91JB00899.
Sobolev
,
S.V.
,
Gresillaud
,
A.
, and
Cara
,
M.
,
1999
,
How robust is isotropic delay time tomography for anisotropic mantle?
:
Geophysical Research Letters
 , v.
26
, p.
509
512
, https://doi.org/10.1029/1998GL900206.
Teanby
,
N.A.
,
Kendall
,
J.-M.
, and
van der Baan
,
M.
,
2004
,
Automation of shear-wave splitting measurements using cluster analysis
:
Bulletin of the Seismological Society of America
, v,
94
, no.
2
, p.
453
463
, https://doi.org/10.1785/0120030123.
Van Wijk
,
J.W.
,
Govers
,
R.
, and
Furlong
,
K.P.
,
2001
,
Three-dimensional thermal modeling of the California upper mantle: A slab window vs. stalled slab
:
Earth and Planetary Science Letters
 , v.
186
, p.
175
186
, https://doi.org/10.1016/S0012-821X(01)00243-6.
Wagner
,
L.S.
,
Beck
,
S.
, and
Zandt
,
G.
,
2005
,
Upper mantle structure in the south central Chilean subduction zone (30° to 36°S)
:
Journal of Geophysical Research
 , v.
110
,
B01308
, https://doi.org/10.1029/2004JB003238.
Wagner
,
L.S.
,
Anderson
,
M.L.
,
Jackson
,
J.M.
,
Beck
,
S.L.
, and
Zandt
,
G.
,
2008
,
Seismic evidence for orthopyroxene enrichment in the continental lithosphere
:
Geology
 , v.
36
, p.
935
938
, https://doi.org/10.1130/G25108A.1.
Wang
,
Y.
,
Forsyth
,
D.W.
,
Rau
,
C.J.
,
Carriero
,
N.
,
Schmandt
,
B.
,
Gaherty
,
J.B.
, and
Savage
,
B.
,
2013
,
Fossil slabs attached to unsubducted fragments of the Farallon plate
:
Proceedings of the National Academy of Sciences of the United States of America
, p.
5342
5346
, https://doi.org/10.1073/pnas.1214880110.
Worthington
,
J.R.
,
Hacker
,
B.R.
, and
Zandt
,
G.
,
2013
,
Distinguishing eclogite from peridotite: EBSD-based calculations of seismic velocities
:
Geophysical Journal International
 , v.
193
, no.
1
, p.
489
505
, https://doi.org/10.1093/gji/ggt004.
Xie
,
J.
,
Ritzwoller
,
M.H.
,
Brownlee
,
S.J.
, and
Hacker
,
B.R.
,
2015
,
Inferring the oriented elastic tensor form surface wave observations: Preliminary application across the western United States
:
Geophysical Journal International
 , v.
201
, p.
996
1021
, https://doi.org/10.1093/gji/ggv054.
Yang
,
B.B.
,
Liu
,
K.H.
,
Dahm
,
H.H.
, and
Gao
,
S.S.
,
2016
,
A uniform database of teleseismic shear-wave splitting measurements for the western and central United States: December 2014 update
:
Seismological Research Letters
 , v.
87
, p.
295
300
, https://doi.org/10.1785/0220150213.
Yang
,
Y.
, and
Forsyth
,
D.W.
,
2006
,
Rayleigh wave phase velocities, small-scale convection, and azimuthal anisotropy beneath southern California
:
Journal of Geophysical Research–Solid Earth
 , v.
111
,
B07306
, https://doi.org/10.1029/2005jb004180.
Yang
,
Y.
, and
Forsyth
,
D.W.
,
2008
,
Attenuation in the upper mantle beneath Southern California: Physical state of the lithosphere and asthenosphere
:
Journal of Geophysical Research. Solid Earth
 , v.
113
,
B03308
, https://doi.org/10.1029/2007JB005118.
Zandt
,
G.
,
2003
,
The Southern Sierra Nevada drip and the mantle wind direction beneath the Southwestern United States
:
International Geology Review
 , v.
45
, p.
213
224
, https://doi.org/10.2747/0020-6814.45.3.213.
Zandt
,
G.
, and
Carrigan
,
C.R.
,
1993
,
Small-scale convective instability and upper mantle viscosity under California
:
Science
 , v.
261
, no.
5120
, p.
460
463
, https://doi.org/10.1126/science.261.5120.460.
Zandt
,
G.
,
Gilbert
,
H.
,
Owens
,
T.J.
,
Ducea
,
M.
,
Saleeby
,
J.
, and
Jones
,
C.H.
,
2004
,
Active foundering of a continental arc root beneath the southern Sierra Nevada in California
:
Nature
 , v.
431
, p.
41
46
, https://doi.org/10.1038/nature02847.
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data