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.

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

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.

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.


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

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.

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

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:
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:
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:
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
which can be recombined as
The value that is measured in our inversion is
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
Applying a simple Taylor expansion, where we expand around Dmean = 0 and ignore the Dmean2 term, gives us
If we define the average azimuthal anisotropy in a layer as T = (UFastUSlow)/Ur, then we get

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.

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.

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.

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.

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.

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.

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.

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

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.

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.

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

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

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.

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.

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.

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.

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

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.

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.

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.

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 or access the full-text article on to view the Supplemental Materials.
Science Editor: Shanaka de Silva
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data