JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 113, B03304, doi:10.1029/2007JB005092, 2008
Testing five of the simplest upper mantle anisotropic velocity
parameterizations using teleseismic S and SKS data from the Billings,
Montana PASSCAL array
Huaiyu Yuan,1,2 Ken Dueker,1 and Derek L. Schutt1
Received 4 April 2007; revised 15 October 2007; accepted 6 December 2007; published 13 March 2008.
[1] Five of the simplest parameterizations of upper mantle anisotropy are tested and ranked
for a data set collected on a dense temporary PASSCAL seismic array located 100-kmNE of
Yellowstone. These hexagonal symmetry anisotropymodels possess either one or two layers
with either flat or dipping fast velocity axis (FVA). Recordings from fifteen high quality
direct-S and SKS arrivals are stacked to provide accurate waveform and error estimates.
Source normalization is accomplished using the cross-convolution technique. A direct
Monte Carlo Neighborhood Algorithm is used to map the posteriori model probability
density (PPD) volume. Using the F test, we find that models with purely flat FVA can be
rejected at >97% probability. Our best model (P5) is a two layer dipping FVA
parameterization, albeit the two layer model with one flat and one dipping FVA can only
be rejected at 80% probability. The best model has a lower layer with a N65E FVA strike
and a12 dip (down to the southwest), and an upper layer with a N20WFVA strike and a
47 FVA dip (down to the southeast). The bottom asthenospheric layer FVA strikes
parallel to North America’s absolute plate motion direction and dips opposite to what
passive plate shear of the asthenosphere would predict. The upper lithospheric layer is
consistent with LPO accretion associated with north directed drift of the North American
plate during theMesozoic. Comparison between the SKS-and direct S-wave data sets shows
that the direct S waves improve resolution of the double layer anisotropic model parameters.
Citation: Yuan, H., K. Dueker, and D. L. Schutt (2008), Testing five of the simplest upper mantle anisotropic velocity
parameterizations using teleseismic S and SKS data from the Billings, Montana PASSCAL array, J. Geophys. Res., 113, B03304,
doi:10.1029/2007JB005092.
1. Introduction LPO symmetry development in the crust, mantle lithosphere
and upper mantle can be complex [e.g., Jung and Karato,
[2] Resolving anisotropic velocity variations between the
2001; Kaminski and Ribe, 2002; Holtzman et al., 2003;
crust, mantle lithosphere, and asthenosphere is important to
Katayame et al., 2004].
constrain a variety of processes that shape our Planet’s
[3] Many shear wave splitting methods use an asymptotic
evolution: e.g., lower crustal flow [Bank and Bostock,
low frequency ray theoretical normal-incidence approxima-
2003; Sherrington et al., 2004], lithospheric accretion
tion [Vinnik et al., 1989; Silver and Chan, 1991; Gledhill
[Babuska and Plomerova, 1993; Fox and Sheehan, 2005],
and Gubbins, 1996; Savage, 1999; Davis, 2003] that
slab-induced flow [Fischer et al., 2000; Anglin and Fouch,
extracts apparent splitting parameters (ASP)–either the fast
2005], and plume-lithosphere interaction [Rumpker and
or slow velocity axis strike and splitting time–via a grid
Silver, 2000; Behn et al., 2004]. Yet, issues remain that
search. In choosing whether to use a fast or slow velocity
strongly limit the constraints provided by shear wave
axis parameterization, the fast axis parameterization has
anisotropy imaging that utilizes shear wave birefringence
been chosen because pyrolitic composition upper mantle
[Fischer et al., 2005; Fouch and Rondenay, 2006]: e.g., the
typically displays fast axis LPO. The ASP uncertainties are
non-linearity of the data kernels makes model appraisal
approximated by fitting a 2-D Gaussian function, centered
difficult, the lack of available high-fold data sets from
on the model misfit minimum. If the ASP is invariant with
seismic arrays makes accurate isolation of anisotropic
respect to shear wave polarization, then a single layer of
signals from signal generated noise difficult, and olivine
anisotropy with a flat FVA is the simplest and presumably
best model. However, if multiple anisotropic layers with flat
FVA [Silver and Savage, 1994; Ozalaybey and Savage,
1Department of Geology and Geophysics, University of Wyoming,
1994; Levin et al., 1999] or a single layer of dipping FVA
USA.
2Now at Berkeley Seismological Lab, UC Berkeley, Berkeley, exist [Hartog and Schwartz, 2000; Levin et al., 2002], then
California, USA. the ASP are modulated by lower order sinusoidal terms.
With this methodology, the analysis devolves into assessing
Copyright 2008 by the American Geophysical Union. the statistical fits through the ASP error estimates.
0148-0227/08/2007JB005092$09.00
B03304 1 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 1. Billings array topography, seismic broadband station locations, and shear wave raypath
piercing points at 150 km depth. The SKS and direct-S wave piercing points are shaded as shown in the
legend. The gray dashed line denotes the strike of the Archeanaged Madison Mylonite shear zone
exposed in the Mountains just north of the Yellowstone Caldera.
[4] In this study, a new method is developed to discrim- graphic sub-setting of our data set, and comparison of
inate between multiple layer models of upper mantle an- models using the direct-S and SKS data subsets.
isotropy. We use the waveform data set from the high [5] Shear wave splitting measurements using Silver and
density Billings array (Figure 1) to test and statistically Chan [1988] eigenvector analysis show a nearly uniform
rank five of the simplest upper mantle anisotropy models. FVA strike (51) and split time (0.8 s) across the array
As a sub-project of the Continental Dynamics Yellowstone (Figure 2). Some systematic variations in the FVA orienta-
array project, the Billings array was designed to record tion at the northern end of the array are found, but the
teleseismic data in a region 100 km NE of the current maximum deviation from the mean FVA is only 12. On the
location of the Yellowstone hot spot located within Yellow- basis of the relative uniformity of the shear wave splitting
stone Park. The relatively small 150 km by 110 km aperture analysis, our further analysis of these ten SKS/SKKS waves
of the array permits coherent stacking of teleseismic shear and an additional five direct S-waves is performed by
wave signals from all 30 broadband stations to accurately stacking the entire array to one stack waveform pair for
isolate anisotropic S-wave signals with robust error esti- each of our 15 events.
mates. The stacking process also attenuates the signal
generated noise and small variations in the upper mantle 2. Data and Methods
anisotropy within the Billings array. To source normalize
2.1. Shear Wave Data Windowing and Polarization
the different events, the cross-convolution technique is used
[Menke and Levin, 2003]. A direct Monte-Carlo search [6] Our teleseismic data set consists of ten SKS/SKKS
method, the Neighborhood Algorithm (NA) [Sambridge, and five direct-S arrivals (Table 1) that are reasonably well
1999a, 1999b], is used to map the posterior model proba- distributed with respect to the events back-azimuth, inci-
bility density volumes (PPD) from which maximum likeli- dence angle, and wave polarization direction (Figure 3). The
hood models may be found. This permits the 1- and 2-D addition of the direct-S events significantly increases the
marginal distribution functions to be assessed quantitatively. diversity of incidence angle and polarization sampling
Model robustness is tested using synthetic data sets, geo- (Figure 3). The increased incidence angle of the direct-S
waves is particularly useful to improve the resolution of
2 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 3. SKS and direct-S wave event distribution with
respect to incidence angle (at 150 km depth) and back-
azimuth. The orientation of the line segments for each event
denotes the polarization of the principle component (S1) of
the wave. For the SKS events this observed polarization is
within 3 of the theoretical back-azimuth (radially polar-
ized). Most of the direct S-events are dominately SH
Figure 2. Billings array shear wave splitting measure- polarized.
ments using eigenvector analysis [Silver and Chan, 1988]
for the SKS/SKKS events in Table 1. The line segment
lengths are proportional to the split time and the line stations recording each event is used. The wavefield param-
segments are parallel to the fast velocity axis. The peak-to- eters constrained via the principle component analysis are:
peak split time is 0.8 s. The thick black arrow shows the the back-azimuth, incidence angle, and the dominant polar-
absolute plate motion direction of the North America Plate ization direction. The dominant polarization direction is
[Gripp and Gordon, 2002]. denoted as S1 with the secondary polarization direction
S2 defined as orthogonal to S1. Principle component
dipping FVA axis [Chevrot and van der Hilst, 2003]. The parameter uncertainties are estimated by bootstrap resam-
waveforms are zero-phase band-pass filtered with a two- pling [Efron and Tibshitani, 1986] of the 28–30 three-
pole Butterworth filter with corner frequencies of 0.015 and component stations that recorded each event (see an exam-
0.15 Hz. ple in Figure 4). The measured polarization of the S1
[7] To estimate the dominant polarization of the direct-S component of the SKS waves is found to have a mean
and SKS waves, a principle component analysis of the 30 value within 3 of the theoretically radially polarized SKS-
Table 1. Event Table
Event BAZ, () Event Year-Day Phase Latitude, () Longitude, () Depth, km Polarization, () Slowness, s/km
0.3 2000-133 SKS 35.98 70.66 108 357 0.046
141.4 1999-325 S 21.75 68.78 101 72 0.097
235.3 2000-228 SKS 31.51 179.73 358 231 0.064
235.3 2000-228 SKKS 31.51 179.73 358 232 0.043
235.3 2000-228 S 31.5 179.73 358 162 0.075
240.5 2000-166 SKS 25.52 178.05 605 238 0.045
240.5 2000-166 SKKS 25.52 178.05 605 238 0.065
240.5 2000-166 S 25.52 178.05 605 355 0.075
264.1 2000-008 SKS 9.81 159.81 33 265 0.046
273.3 2000-037 SKS 5.84 150.88 33 270 0.043
291.2 2000-057 SKS 13.8 144.78 132 296 0.052
307.1 2000-161 S 30.49 137.73 307 280 0.089
359.9 2000-199 SKS 36.28 70.82 141 350 0.046
359.9 1999-312 SKS 36.52 71.24 228 358 0.046
359.9 1999-312 S 36.52 71.24 228 248 0.076
3 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
This polarization anomaly increases as the angle between
the FVA and the raypath departs from perpendicularity. For
the mean 19 incidence angle of our direct S-rays at 150 km
depth (Figure 3), and the 12–28 FVA dip found by model
parameterization P4 and P5 (described in section 2.4) for the
bottom layer, the transverse wave polarization anomaly is
predicted to be 4–8. Thus this is a small error that we
have ignored.
[8] The waveforms are visually windowed to minimize
the influence of secondary converted and reverberated
phases from the crust-mantle boundary. For each individual
event, a common window length is found that starts at the
shear wave signal onset to avoid S-P precursors and the
length of the waveform segment is 10–20 s depending on
the magnitude of the earthquake. This window length
excludes all free surface Moho reflection except the Spms
phase which is predicted to arrive 16–19 s after the S-wave
arrival. This arrival time is estimated based on combined
receiver function [Yuan et al., 2006] and diffusive/ballistic
surface wave imaging beneath this array [Stachnik et al., in
review.]. To account for relative traveltime variations asso-
ciated with small isotropic velocity variations beneath the
array, the multichannel cross-correlation method [MC,
vanDecar and Crosson, 1990] is used to measure S-wave
arrival times. Relative traveltimes are measured on the S1
component, and the maximum relative traveltime variation
for the S1 component is 1.6 s. After applying the MC
calculated time-shifts to the S1 and S2 components recorded
Figure 4. SKS event 2000-133 (Table 1) recording and
polarization estimation. (a) Traveltime residual adjusted S1
component waveforms. (b) Residual adjusted S2 component
waveforms. (c) S1 polarization direction histogram derived
from all recording stations. (d) Incidence angle histogram.
(e) The S1 and S2 stack traces with one-sigma error bar
shaded in light gray. The thin dashed lines in (a), (b) and (e)
define the signal measurement window. The polarization
measurements are derived from 200 bootstraps of the
principle component analysis for all the stations. The thick
black line in (c) is the theoretical event back-azimuth and
the dashed line is the mean estimate of the S1 polarization
direction. For the SKS/SKKS arrival, the event back-
azimuth is very close to the S1 direction (Table 1). The
dashed line in (d) is the mean estimate of the incidence
angle from the histogram. Note that the S1 component has
high signal coherence while the split energy present on the
S2 components is less coherent. The relative incoherence of
the S2 component is attributed to signal generated noise
created by the crust and overlying sedimentary layer.
Figure 5. Direct-S and SKS waveform pairs. The five
direct-S waveform pairs are shown as dotted lines and the
waves (Table 1). Noteworthy is that for dipping FVA
ten SKS waveform pairs as solid lines. The thick and thin
anisotropy, the polarization directions of the fast and slow
lines are the S1 and S2 components, respectively. The
wave components are no longer orthogonal and a systematic
shaded region about each trace is one standard deviation
angular deviation (polarization anomaly) of the transverse
estimated from stacking all the stations recording each
wave is predicted [e.g., Chevrot and van der Hilst, 2003].
event.
4 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 6. Five simple anisotropy models. The P1 and P2 parameterizations are single layer anisotropy
models with P1 having a flat FVA and P2 having a potentially dipping FVA. Model P3 is a two layer
model with a flat FVA in both layers. Model P4 has the upper layer with flat FVA and a potentially
dipping FVA in the lower layer. Model P5 has a potentially dipping FVA in both layers. Note all
anisotropic layers are fixed to be 100 km thick. The crustal thickness is fixed to the array mean of 52 km
(Yuan et al., in prep.).
by each event, the waveforms are linearly stacked to obtain [11] Following Frederiksen et al. [2003], we assess the
15 event stack S1-S2 waveform-pairs that are subsequently misfit associated with equation (1) using cross correlation
modeled (Figure 5). and the L1 and L2 norms. We find that these different misfit
[9] Strategies to minimize potential source-side anisotro- norms produce very similar results. The L2 norm is used in
py signal in direct S-arrivals include: explicit correction for this research because this norm can be used as input to the F
known source-side anisotropy [e.g., Yang and Fischer, test to compare the significance of the variance reductions
1994; Anglin and Fouch, 2005]; or use of deep-focus events between model parameterizations that possess different
[Savage et al., 1990; Fischer and Yang, 1994; Hartog and number of degrees of freedom. Search for the best models
Schwartz, 2000; Long and van der Hilst, 2005]. For this are quantified by minimizing the L2 misfit of the cross-
study, a total of 12 direct-S events with depths >100 km are convolved traces,
initially selected for analysis. To cull these direct-S arrivals, n o
the waveforms are source normalized (see below) with our min ½s2ðmÞ*S1 S2*s1ðmÞ2 : ð2Þ
best one layer anisotropy model found using the SKS-only
data set. Then, event stack waveforms with large misfits are
assumed to have large source-side anisotropy contamination [12] The synthetic seismograms are calculated using ray
and discarded. This culling produces a final data set of five theory for a transverse anisotropic medium [Frederiksen
direct-S stack waveforms with source depths >228 km and Bostock, 2000]. The 3-D traveltime equation of Diebold
except one event at 101 km depth (Table 1). [1987] is adopted in this method, and the phase amplitudes
are computed using the refection and transmission matrices
2.2. Source Normalization and Synthetic Seismograms for planar interfaces separating homogeneous anisotropic
[10] Our 15 event waveform pairs all have high SNR and media. This ray-based approach allows us to selectively
small error estimates; therefore we first normalize the model the first-order free-surface Moho reverberations and
amplitudes of the waveform pairs by the sum of squared to assess the influence of these small secondary arrivals
amplitudes of S1 and S2 components. This ensures that upon our model fits. The anisotropy is parameterized as the
every event waveform pair will have equal weight in the percent variation in P- and S-velocity perturbation (dVp and
minimization of the cost function. Given that each event has dVs) and h, in the manner described by Farra et al. [1991].
a different source time function, we source normalize our In this study, dVp is set equal to dVs and h is set to 1.03
waveform pairs using the cross-convolution method [Menke where h is the anisotropic parameter that controls the
and Levin, 2003; Levin et al., 2006]. The method is velocity variation for directions that are not parallel or
predicated on the reasonable assumption that the source perpendicular to the symmetry axis (it creates small dimples
wavelet is well described as a linearly polarized wave before on the prolate spheroid velocity representation). Thus only
it splits into two approximately orthogonal waves during one anisotropic parameter, the percent variation of velocity
transit through an anisotropic medium. An approximate
cross-convolved relation between the observed S1 and S2
polarized S-waves (Figure 5) and synthetic s1 and s2
Table 2. Synthetic Test Models
waveforms is
Layer 1 Layer 2
S1*s2ðmÞ S2*s1ðmÞ ð1Þ Anis., FVA FVA Anis., FVA FVA
Models (%)a Strike, Dip, (%)a Strike, Dip,
where m is the synthetic anisotropic velocity model and * SD1 4 65 - - - -
SD3 4 5 - 4 70 -
indicates convolution. This equation simply states that when SD5 4 20 60 6 60 10
m is near the true model, the convolution of the observed S1 aFor 100-km thick layer. The FVA strike is rendered positive clockwise
trace with its orthogonal synthetic is approximately equal to with respect to North. The dip is positive downwards with respect to
the convolution of the observed S2 trace with its orthogonal horizontal along the FVA strike direction. Negative dip means it dips
synthetic. towards 180 of the strike.
5 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
(referred to as anisotropy magnitude hereinafter) is mod- analysis shows that one dof requires 0.8 s of S-wave
eled. The S1 and S2 polarization responses are calculated signal. Thus for 15 events with an average signal length of
via rotation of the SV and SH impulse responses. This 10 s, our data set contains 120 dof.
approximation assumes that the SV and SH raypaths are [15] In the second stage of the NA modeling, quantitative
the same which is a good weak-anisotropy approximation model information is extracted using the Bayesian approach
[e.g., Crampin, 1977; Rumpker and Silver, 2000]. [Sambridge, 1999b]. In this approach, the information
contained in the model ensemble is represented by the
2.3. Neighborhood Algorithm Sampling and Appraisal posterior probability density functions (PPD). The PPD
[13] The Neighborhood Algorithm (NA) [Sambridge, functions are then used to calculate quantities such as the
1999a, 1999b] is used to map the posterior model proba- model expectation (mean) and the 1- and 2-d posterior
bility density volume (PPD) associated with our five aniso- marginal probability density functions. Assuming a multi-
tropic velocity model parameterizations. NA is a direct variate Gaussian probability function with a uniform prior
Monte Carlo search method which has been applied to probability distribution r(m), the PPD is given by
many non-linear inversion problems. Readers are referred
to Sambridge [1999a, 1999b] for the details of the method. v
PðmÞ ¼ krðmÞexp c2ðmÞ ð5Þ
The NA methodology consists of two stages. First, a search 2 v
stage is guided by a random ‘‘walk’’ through the nearest
neighbor Voronoi cells that discretize the model parameter where k is a constant that normalizes the total probability to
2
volume. Two controlling parameters, ns, the number of unity, v is dof, and cv is the reduced chi-square value
models generated in each iteration, and nr, the number of (equation (3)). N-dimensional Bayesian integrals are used to
lowest misfit models in which the random walks are calculate the 1- and 2-D PPD marginals (note: we use the
performed, control the spatially uniformity of the model short-hand term ‘‘marginals’’ to refer to the set of 1- and
search. These search parameters have been both set to 5, 10, 2-D marginal PPD herein). To improve the accuracy of the
and 100, and we find they all can provide an extensive N-dimensional integrations, a Gibbsian resampling of the
search of the model space so that the global minimum is not probability space is used [Sambridge, 1999b].
missed [Sambridge, 2001]. The total number of models
generated by the NA is between 103–106 depending on the 2.4. Five Model Parameterizations
dimensionality of the model space which varies between 2 [16] The main goal of this paper is to test our 15 SKS-
and 6. and direct S-wave waveform pairs against five of the
[14] The reduced chi-square misfit between the cross simplest models of upper mantle anisotropy. Given that
convolved waveforms is defined as: our error estimates are well-characterized by the standard
deviation statistic, the F test is used to rank the significance
X N i ð Þi i ð Þi 2 of the relative variance reductions between the five different2ð Þ ¼ 1 S1 *s2 m S2 *s1c m
m ; ð3Þ model parameterizations. For our anisotropic velocity mod-v dof ¼ si 2 2i 1 þ2si si i el parameterizations, we assume a hexagonal anisotropyx x y þ sy
system with a fast velocity axis, which is a widely accept-
able representation for olivine lattice preferred orientation in
where i indexes the events number, N is number of events,
the upper mantle [e.g., Park and Levin, 2002]. The five
dof is the number of degrees of freedom, the star operator is
anisotropic model parameterizations are denoted P1–P5
convolution, S1, S2 and s1, s2 are the stacked waveforms
(Figure 6) and can be described as follows: P1 and P2 are
and synthetics, respectively, and a sum over the time series
single layer models with either a flat or potentially dipping
is implicitly assumed. The standard deviation of the two
FVA; P3 is a double layer model with a flat FVA in both
cross-convolved traces is defined as
layers; P4 has an upper flat FVA layer and lower potentially
¼ ð Þ ¼ ð Þ ð Þ dipping FVA layer; P5 has two layers with potentiallysix s1 m *si i is1 and sy s2 m *ss2; 4 dipping FVA.
[17] The thickness of each layer is fixed to 100 km. The
where siS1 and s
i
S2 are the standard errors of the S1 and S2 layer thickness is fixed because our data set has very limited
components. The errors for the two cross-convolved terms sensitivity to the trade-off between layer thickness and
are considered perfectly correlated and hence the error cross anisotropy magnitude. Therefore the anisotropy magnitude
term in the denominator of equation (3) is present. The dof found for each layer is only meaningful if one believes that
for each event waveform pair is calculated using the spectral the 100 km thick layers are close to their true values. We
bandwidth of the data [Silver and Chan, 1991]. This
Figure 7. Synthetic model PPD marginals constructed using input parameters in Table 2. (a) SD1 data set. (b) SD3 data
set. (c) SD5 data set. Each subplot of the PPD marginals is labeled by the model parameterization (P1–P5). The synthetic
model parameter input mimics the real data inversion results (shown in section 3.2). The diagonal of each lower-half
‘‘matrix’’ subplot shows the 1-D PPD marginals for the model parameters. The 1-D PPD marginals are found by integrating
over the other model parameters and hence properly accounts for correlated model errors. The lower ‘‘matrix’’ shows the
2-D PPD marginals, which show the correlation of any pair of the model parameters. For the 2-D marginals, the darker
shading corresponds to higher probabilities and a contour line outlines the 80% probability contour. The model parameter
names and values are labeled along the bottom and left side of the matrices. The strike of the FVA is rendered in degrees
positive clockwise with respect to North.
6 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
7 of 20
Figure 7
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
8 of 20
Figure 7. (continued)
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
9 of 20
Figure 7. (continued)
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Table 3. Synthetic Test: Results of SD1
Layer 1 Layer 2
FVA Strike,
Anis., (%)a () FVA Dip, () Anis., (%)a FVA Strike, () FVA Dip, () DOF c 2 RPbv
P1 4.0 [3.9, 4.2] 66 [65, 67] - - - - 78 1.50 -
P2 4.3 [3.6, 4.4] 65 [62, 68] 1 [10, 3] - - - 77 1.51 52.2%
P3 3.3 [3.1, 3.7] 67 [64, 70] - 0.5 [0.2, 0.8] 65 [49, 68] - 76 1.53 53.8%
P4 4.1 [3.6, 4.5] 67 [62, 72] - 0.8 [0.4, 1.6] 24 [23, 38] 68 [81, 37] 75 1.55 56.6%
P5 4.3 [3.8, 5.5] 65 [55, 76] 2 [13, 11] 1.2 [0.0, 1.6] 59 [68, 52] 64 [52, 75] 74 1.57 58.4%
aFor 100-km thick layer.
bRejection probability. Brackets show the 80% probability region. P 1 is the correctly parameterized model.
note that mantle xenoliths from eastern Montana suggest a and SD3 correspond to the one- and two-layer models with
120–140 km thick lithosphere [Carlson et al., 2004]. flat FVA anisotropy and SD5 corresponds to the two-layer
Therefore given the 48–54 km thickness of the crust model with dipping FVA anisotropy. We choose to test these
beneath the Billings array [Yuan et al., 2006] at least 70– three models because SD1 and SD3 are most commonly
90 km of lithospheric mantle is present based on the tested anisotropic model parameterizations, while SD5 is
xenolith studies. If the anisotropic layers beneath the array our most complex model parameterization. The anisotropy
have sharp anisotropic velocity contrasts, then modeling of parameters of the synthetic models (Table 2) are chosen
the S-P precursors could provide useful information to to simulate the real-data modeling results (shown in
constrain the trade-off between layer-thickness and anisot- section 3.2). A 4% and 6% anisotropy magnitude is used
ropy magnitude. for the upper and lower anisotropic layers. The synthetic
[18] The crustal thickness and shear velocity has been waveforms are generated using Frederiksen and Bostock
constrained in our model parameterization using a combi- [2000] ray theoretical code. Random white noise filtered to
nation of teleseismic P- and S-wave traveltime measure- the bandwidth of the S-wave amplitude spectrums is scaled
ments [Yuan and Dueker, 2005], P wave receiver function to match the observed data signal-to-noise ratio and added
analysis [Yuan et al., 2006], ballistic Rayleigh wave imag- to the synthetic waveforms. The sum of the number of
ing [Schutt et al., 2008] and combined diffusive and ballistic degrees of freedom for each synthetic data set is about 80.
Rayleigh wave imaging [Stachnik et al., in review]. These [20] To show how these three synthetic data sets are fit by
new seismic constraints show that a relatively uniform our five model parameterizations, the PPD marginals for
crustal thickness of 48-54 km is present beneath the array each of the synthetic data sets are calculated. To evaluate
with an average shear wave velocity of 3.74 km/s and a Vp/ these synthetic results, three relations between the synthetic
Vs value of 1.78 [Yuan et al., 2006]. These observations are model from which a data set is calculated and the model
consistent with the Deep Probe refraction line that sampled parameterization used, are defined: 1) a correct-parameter-
beneath the Billings array [Henstock et al., 1998; Gorman et ized model where the synthetic data model and the model
al., 2002]. This well-constrained crustal thickness and parameterization are the same; 2) an under-parameterized
velocity information means that the relatively long-period model where the synthetic data model has higher dimen-
Spms free surface reverberation present in our signal time sionality that the model parameterization; 3) an over-
window is modeled. parameterized model where the synthetic data model has
lower dimensionality than the model parameterization.
3. Results To assess the significance of variations in the data misfit
between each of the synthetic data sets, rejection probabilities
3.1. Synthetic Data Tests
are calculated via the F-test which uses the reduced chi-
[19] Synthetic tests using the same event parameters square values of each synthetic model. The rejection proba-
(back-azimuth, incidence angle, and polarization) as the bilities are calculated with respect to the lowest misfit model
true data set (Table 1) are conducted to investigate the and a 50% rejection probability means that a particular
ability of our data set coverage to constrain anisotropic model is no more probable than the lowest misfit model.
models. Three synthetic waveform data sets (SD1, SD3, and [21] Inspection of the synthetic test results (Figure 7 and
SD5) are created using the parameters in Table 2. Note SD1 Tables 2, 3, 4, 5) reveals the following conclusions. As
Table 4. Synthetic Test: Results of SD3
Layer 1 Layer 2
FVA Strike,
Anis., (%)a FVA Strike, () FVA Dip, () Anis., (%)a () FVA Dip, () DOF c 2 bv RP
P1 2.1 [2.0, 2.2] 35 [33, 38] - - - - 78 2.0 87%
P2 2.0 [1.7, 2.2] 35 [32, 43] 16 [41, 9] - - - 77 2.0 89%
P3 3.9 [3.0, 4.0] 2 [12, 7] - 3.3 [2.9, 4.4] 69 [61, 75] - 76 1.54 -
P4 3.9 [3.2, 5.3] 16 [27, 2] - 3.9 [3.4, 4.7] 72 [63, 83] 42 [52, 22] 75 1.56 52%
P5 3.9 [3.0, 5.6] 20 [36, 6] 31 [39, 43] 4.5 [3.8, 5.6] 62 [50, 87] 13 [35, 5] 74 1.57 54%
aFor 100-km thick layer.
bRejection probability. Brackets show the 80% probability region. P3 is the correctly parameterized model.
10 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Table 5. Synthetic Test: Results of SD5
Layer 1 Layer 2
Anis., (%)a FVA Strike, () FVA Dip, () Anis., (%)a FVA Strike, () FVA Dip, () DOF c 2 RPbv
P1 3.7 [3.6, 4.0] 65 [62, 67] - - - - 78 2.3 94%
P2 4.1 [3.6, 4.7] 68 [63, 71] 42 [44, 36] - - - 77 1.7 61%
P3 0.0 [0.0, 0.7] 54 [31, 65] - 3.6 [3.1, 4.3] 65 [59, 68] - 76 2.3 94%
P4 3.9 [2.5, 4.9] 65 [56, 76] - 0.0 [0.0, 0.6] 31 [42, 90] 75 [66, 89] 75 2.1 88%
P5 4.1 [2.7, 5.2] 6 [27, 8] 27 [42, 19] 8.5 [7.5, 10] 61 [46, 77] 34 [54, 17] 84 1.6 -
Anis.: Anistropy.
aFor 100-km thick layer.
bRejection probability. Brackets show the 80% probability region. P5 is the correctly parameterized model.
expected, the SD1 data set (single layer with flat FVA) can marginals are generally compact and uni-modal, indicating
be equally well fit by all five model parameterizations good resolution of the model parameters with little model
(Figure 7a and Table 3). This equally good model fitting parameter correlation (Figure 8). Most notable is that the 1-
is quantified by the F-test which shows that the rejection and 2-D marginals widen for the more complex models
probabilities for the over-parameterized models are negligi- consistent with the increasing non-linearity and non-unique-
ble (52% to 58%). Remarkable is that the over-parameter- ness associated with the more complex anisotropic models.
ized models (P2–P5) remained stable and do not produce [24] Ranking of the models is done using the F test whose
spurious results such as large anisotropy magnitudes or utility depends on the accuracy of our waveform error
large FVA dip. The SD3 data set (two-layer model with estimates (Figure 5). Inspection of the chi-square values
both flat FVA, Figure 7b and Table 4) is best fit by the and associated F test based model rejection probability
correctly parameterized model P3. As expected, the under- levels (Figure 9) shows that the P5 model provides the best
parameterized models have high rejection probabilities L2 fit to the cross-convolved waveforms (Figure 10). The
(87% and 89%) and the over-parameterized models have reduced chi-square value of this model is 1.9, indicating a
low rejection probabilities (52% and 54%). The SD5 data slight under-fitting of the waveforms (Table 6). This under-
set (two layers with dipping FVA) is best fit by the correctly fitting of the data is expected given the approximations in
parameterized model P5, and parameterizations P1–P4 can our forward problem, for instance source-side contamina-
be rejected at 94%, 61%, 94%, and 88% probability tion of the direct S-wave data set (addressed below). The
(Figure 7c and Table 5). P1–P3 model parameterizations have larger reduced chi-
[22] The synthetic tests reveal that the correctly parame- square values and these models can be rejected at 98%,
terized models can be recovered by the full data set. In all 89%, and 97% probability levels. Given that the P4 model
cases, the under-parameterized models can be rejected with can only be rejected at 80% probability, we consider this
high probability while the over-parameterized models are model to be acceptable. Note that the P4 and P5 model
always acceptable (Tables 2–5). The PPD marginals for the parameter uncertainty bounds do overlap.
synthetic tests (Figure 7) demonstrate that the double layer [25] Model parameter uncertainties are estimated by cal-
model parameterizations have 1- and 2-D PPD marginals culating the 80% probability bounds from the 1-D PPD
that are less compact and in some cases multimodal. This marginals (Figure 8). Noteworthy is that these 1-D PPD
increased uncertainty of the model parameters is due to the marginals fully account for correlation between the model
fact that two layer models with dipping FVA produce parameters. Inspection of Table 6 shows that commonalties
greater waveform complexity with respect to single layer exist between the different models. The most consistent
models [e.g., Ozalaybey and Savage, 1994]. In practice, this result is the existence of a layer with its FVA directed
data kernel complexity produces non-uniqueness of the toward the ENE (i.e., N64–N75). This direction is 19
model space: i.e., there are a range of models that produce larger than the mean FVA of 51 found by our SKS
nearly identical waveforms within the error bars of our eigenvector analysis (Figure 2), but is consistent with the
observed waveforms. The conclusions from our synthetic absolute plate motion direction [Gripp and Gordon, 2002].
tests are that: 1) the input synthetic models for the noisy We suggest this FVA difference results from the biasing
synthetics can be recovered without model parameter insta- effects of the more complex anisotropy models found by our
bilities, and 2) the under-parameterized models can be waveform analysis upon the eigenvector SKS analysis. For
rejected in most cases using the F test. the models that permit a dipping FVA, the bottom layer
FVA consistently dips at 12 to 28 (down to the SW)
3.2. Real-Data PPD Results with respect to horizontal. The anisotropy magnitude range
[23] NA mapping of the PPD marginals with the com- for the upper 100 km thick layer is 2.9–3.6% and the
bined direct-S and SKS data set is performed for all five 100 km thick bottom layer has a larger anisotropic magni-
model parameterizations (P1–P5). Inspection of the PPD tude range of 5.9–7.5%.
marginals for the different parameterizations shows that the
Figure 8. Real-data 1- and 2-D PPD marginals. See Figure 7 for figure layout description. Each subplot of the PPD
marginals is labeled by the model parameterization (P1–P5). The model values and uncertainty bounds for the P1–P5
model parameterizations are listed in Table 6.
11 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
12 of 20
Figure 8
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 9. F test results showing the model rejection probabilities with respect to the lowest misfit
model P5. The thin dashed line plotted is the F-inverse curve using the average number of degrees of
freedom of the five parameterizations. The two flat FVA model (P1 and P3) can be rejected at >97%
probability. Because the F test is a comparison of two models misfit variances, a rejection probability of
50% mean that a model is not a statistical improvement with respect to the lowest misfit model.
[26] For our lowest misfit model parameterization (P5),
the upper layer has an anisotropic velocity magnitude of
3.1% –4.1% with an FVA dip of 36 to 59 (down to
the SSE) and a FVA strike of N29W-N10W. The lower
layer has an anisotropic velocity magnitude of 6.4%–8.0%
with a FVA strike of 52–74 and a FVA dip between 20
and 1 (down to the SW). This best model is very similar
to the P4 model which is also an acceptable model. Our
most robust conclusion is that model parameterizations
without a dipping FVA do not provide acceptable fits to
our waveform data set.
3.3. Direct S-Wave Resolution Improvement
[27] An important question with respect to our analysis is
how much the five direct S-events contributed to the
resolution of our model parameters. Of most concern is
whether source-side anisotropy has significantly contami-
nated our direct-S waveforms. Theoretical analysis demon-
strates that the greater incidence angles of direct S arrivals
with respect to SKS arrivals (Figure 3) provide greater
sensitivity to dipping FVA anisotropy [Chevrot and van
der Hilst, 2003]. In addition, the direct S-waves provide a
range of SH-like polarizations not provided by the SV-
polarized SKS waves. As mentioned previously, our anal-
ysis finds that the P1 and P3 model parameterizations
without a dipping FVA can be rejected at high probability Figure 10. Maximum likelihood model (Table 6) pre-
(Table 6). dicted cross-convolved waveform pairs. The dotted lines
[28] Figure 11 shows the one and two dimensional PPD denote the five direct-S waveform pairs and solid lines are
marginals for three subsets of our waveforms: the five direct the SKS waveform pairs. The direct S-waves are fit as well
as the SKS waves.
13 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Table 6. Maximum Posterior Probability Density (PPD) Models
Layer 1 Layer 2
Anis., (%)a FVA Strike, () FVA Dip, () Anis., (%)a FVA Strike, () FVA Dip, () DOF c2 RPbv
P1 3.8 [3.6, 4.0] 64 [63, 65] - - - - 118 2.8 98%
P2 3.7 [3.4, 3.8] 65 [63, 67] 28 [32, 25] - - - 117 2.4 89%
P3 2.7 [2.1, 3.1] 2 [5, 8] - 5.9 [5.3, 6.2] 75 [69, 80] - 116 2.7 97%
P4 2.9 [2.6, 3.4] 10 [18, 3] - 5.5 [5.3, 6.1] 71 [66, 77] 28 [30, 9] 115 2.2 80%
P5 3.6 [3.1, 4.1] 20 [29, 10] 47 [59, 36] 7.5 [6.4, 8.0] 65 [52, 74] 12 [20, 1] 114 1.9 -
aFor 100-km thick layer.
bRejection probability. Brackets show the 80% probability region. The FVA strike is rendered positive clockwise with respect to North. The dip is
positive downwards with respect to horizontal along the FVA strike direction. Negative dip means it dips towards 180 of the strike.
S-waveforms, the ten SKS waveforms, and the combination NA modeled. Comparison of the PPD associated with these
of both waveforms sets. These three data sets are fitted to two sub-arrays shows that the 1- and 2-D PPD marginals
our five different model parameterizations. For the simplest strongly overlap (Figure 12). This suggests that the aniso-
model parameterization (P1), all three data sets produce 1- tropic structure beneath the Billings array is reasonably
and 2-D PPD marginals that overlap. This demonstrates that spatially uniform within the sampling volume.
the direct-S and SKS data sets individually possess coherent [31] The synthetic waveform pairs associated with the
anisotropic signal with which to constrain the single aniso- best model for each of the five anisotropic model parameter-
tropic layer. For the single layer dipping FVA parameteri- izations are shown in Figure 13. This plot shows that the S1
zation (P2), the FVA strike found is the same as the P1 components are only marginally affected by the different
model (Table 6), but the anisotropy magnitude found by the anisotropy models while the S2 components are greatly
direct-S-waves is 2% greater than the SKS or combined data affected by the different anisotropy models. Obviously,
sets. For the double layer model with flat FVA parameter- modeling these S2 response variations requires an accurate
ization (P3), the strike of the bottom layer FVA is well set of waveforms with well characterized error estimates.
defined by all three data sets, but the direct S-wave data set [32] To assess how the misfit between the data and five
shows variable marginals. best models is distributed with respect to the different
[29] The two most complicated model parameterizations events, the cross-convolved waveforms and their misfits
(P4 and P5) reveal the following general characteristics with are presented (Figure 14). Because all the events are initially
respect to their PPD marginals. First, the marginals for the equal-weighted in our data fitting, the chi-square values
combined direct-S and SKS data sets are always more have been normalized by the chi-square value associated
peaked and compact than the marginals from either of the with the P1 model parameterization. This permits the trend
individual data sets. Second, the 80% probability bounds in the chi-square values with respect to the P1–P5 models
overlap between the individual SKS and direct-S wave- to be easily assessed. Inspection of the normalized chi-
forms subsets for ten out of the fifteen 2-D marginals. Third, square values for each event and model shows that the more
the five direct-S waveform data set do have the least complicated anisotropy models generally have better fits to
compact PPD marginals. We suggest that the marginals the data (Table 6). Noteworthy is that some event wave-
for the direct S-wave data set are less compact because of forms are more poorly fit by the more complicated anisot-
the existence of the possible source-side anisotropy and the ropy models. We suggest that signal exists in our data set
limited resolving power provided by this small data subset that is not predicted by our forward model. For example, the
(without the SKS/SKKS waveforms). Yet, the five direct-S true anisotropic variations could be more complicated than
events do add significant information as demonstrated by two layers of transverse symmetry anisotropy, perhaps due
the compactification of the 1- and 2-D PPD marginals using to modest lateral heterogeneity in anisotropy. Also, source
the full data set with respect to the SKS-only data set. side anisotropy for the five direct-S events could contam-
inate these waveforms to some degree. Yet, we note that the
3.4. Geographic Variations and Model Misfit cross-convolved direct-S waveforms are not systematically
[30] As a test of robustness of our models with respect the misfit with respect to the SKS arrivals (Figure 14).
geographic size of our array, the Billings array is divided
into a northern and southern sub-array of stations. This
4. Discussion and Conclusions
division results in a northern and southern sub-array with 11
and 19 stations, respectively. From these two geographical [33] The orientation of the two anisotropic layers’ strike
sub-arrays, a set of event stack waveforms are produced and and dip for the best anisotropy model (P5) is superimposed
Figure 11. PPD marginals for the SKS-only, direct-S-only, and full data set (Figure 8). See Figure 7 for figure layout
description. Each subplot of the PPD marginals is labeled by the model parameterization (P1–P5). The marginals for the
SKS, direct-S, and combined data sets are color coded as blue, green and red. Note that both the 1- and 2-D marginals are
made more compact by the addition of the five direct-S waveform pairs. The 2-D marginals show that the direct-S wave
data set alone is insufficient to constrain the P3–P5 model parameters well.
14 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
15 of 20
Figure 11
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
16 of 20
Figure 12
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 13. Maximum likelihood model synthetic waveform pairs for the five model parameterizations
(Table 6). (a) S1 waveforms. (b) S2 waveforms. For each given event polarization values (Table 1), the
waveforms from left to right represent the waveform response of model parameterizations P1–P5. The
‘‘S’’ symbol denotes the direct-S waveform.
upon a P wave tomogram that was produced using the [34] The most common element between our five sim-
Yellowstone and Billings teleseismic P wave traveltimes plest anisotropic velocity models is that one layer (the
measurements [Yuan and Dueker, 2005]. Receiver function bottom layer for the two layer models) has a FVA strike
studies that use local surface wave models to map time to between 64–75 which is parallel to the North American
depth [Stachnik et al., in review] show that the crust beneath absolute plate motion direction. This NE-directed FVA
the Billings array is 48–54 km thick [Yuan et al., 2006]. strike is consistent with the nearby Yellowstone array single
The P wave tomogram shows that the crust is underlain by a layer shear wave splitting measurements [Waite et al.,
relatively high velocity mantle lithosphere that extends to 2005]. An additional result for our dipping anisotropy
120–140 km depth. This lithospheric thickness is consistent models is that the lower layer FVA dips at -28 to -12
with eastern Montana mantle xenolith studies that find a down to the SW with respect to horizontal. Unique inter-
120–140 km thick lithosphere from Eocene kimberlite pretation of the SW dipping FVA in the lower layer is not
eruptions [Carlson et al., 2004]. Thus the lithospheric possible until potential dipping anisotropy is characterized
mantle layer beneath the Billings array is estimated to be around the Yellowstone region. At present, only one other
at least 75 km thick which is sufficiently thick to create a analysis requires a dipping FVA axis for the National
significant anisotropic signal. Below 120–140 km depth, seismic network station HLID near the Yellowstone hot
the mantle is normal within the volume sampled by our S- spot track [Walker et al., 2004]. The most noteworthy
wave raypaths, although the westernmost raypaths come statement to be made with respect to the SW dip of the
within 40 km of the edge of the Yellowstone plume FVA is that this dip is the opposite of what a passive plate
(Figure 15). shear model of LPO evolution predicts [Bokelmann, 2002].
One could also speculate that plume entrained mantle
Figure 12. PPD marginals for the geographically divided north and south array sub-data set. See Figure 7 for figure layout
description. Each subplot of the PPD marginals is labeled by the model parameterization (P1–P5). For the 1-D PPD
marginals, the combined data set PPD marginals are drawn as unfilled lines and the north and south sub-data set PPD
marginals are filled with light gray (north data set) and dark gray (south data set). For the 2-D PPD marginals, the 80%
probability contour for the north and south data sets are filled with light/dark gray and the combined data set probability is
an unfilled contour.
17 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 14. Maximum likelihood model cross-convolved traces and residual traces. The columns labeled
P1–P5 are the five model parameterizations. The cross-convolved traces are the thin lines and the gray
filled traces are the difference between each cross convolved trace pairs (residual trace). The Chi-square
error between each of the cross-convolved waveform-pair, normalized by its P1 Chi-square value, is
given to the left of each trace. The contribution of each cross convolved trace pair to the sum of the P1
Chi-square value is shaded in gray scale.
associated with the Yellowstone plume may be perturbing 2005]. Deschamps et al. [2006] speculates that the north
the FVA dip in the asthenospheric layer beneath the Billings trending lithospheric LPO was frozen into the lithosphere
array. during northward drift of the North American plate during
[35] The N10–20W FVA strike of the upper layer for Mesozoic times.
the P4 and P5 models is broadly consistent with the FVA [36] Given that the P4 model can only be rejected at 80%
strike of the upper layer found by several other anisotropic probability, we consider the -47 south directed dip of the
studies. Love and Rayleigh wave imaging finds that most of upper lithospheric layer for model P5 to not be a robust
the mid-continental region of the United States and Canada feature. With regard to the geologic perspective, the
has a northerly trending FVA in the upper 100 km [Marone Madison Mylonite zone (Figure 15) has been interpreted
and Romanowicz, 2007]. More specifically, beneath the as either an Archean age suture zone [Hoffman, 1989] or an
Billings array the Marone and Romanowicz upper aniso- intracratonic shear zone [Erslev and Sutter, 1990]. Given
tropic layer FVA points at N10W within the error of our the uncertainties in whether the upper layer requires a
upper layer FVA orientation. Rayleigh wave imaging be- dipping FVA, we simply note that the deformation along
neath the east-central North American continent and the the Madison Mylonite zone could have produced a dipping
Slave craton also finds an upper lithospheric layer with a FVA.
north trending FVA [Deschamps et al., 2006; Chen et al., [37] In conclusion, we have developed a new method to
2007; Snyder and Bruneton, 2007]. In addition, two-layer discriminate between the five simplest models of upper
anisotropy modeling in the Slave craton, southern Alberta mantle anisotropy using source normalized and stacked
and southern Wyoming find an anisotropic lithospheric teleseismic SKS and direct-S waveforms. The addition of
layer with a northerly striking FVA underlain by a lower direct-S events has been shown to significantly improve
layer that is approximately parallel to absolute plate motion resolution of the anisotropic model parameters. The direct
[Snyder et al., 2003; Currie et al., 2004; Fox and Sheehan, Monte-Carlo NA search has proven to be an efficient
18 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Figure 15. Comparison of best anisotropic velocity model (P5) and a teleseismic P wave tomographic
model [Yuan and Dueker, 2005]. (a) Map view at 100 km depth. The thick and thin black contours are the
2% and 1% P wave velocity perturbation. The Yellowstone Caldera is outlined in white. The gray dash
line denotes the Madison Mylonite shear zone. The two-headed white arrow shows the strike of the FVA
strike of the upper anisotropic layer. (b) Map view at 200 km with lower layer FVA strike denoted by
two-headed white arrow. (c) Profile along A-A0. The white arrow shows the FVA dip in the top layer. The
two white lines show the volume sampled by the S-wave raypaths (Figure 1). (d) Profile along B-B0. The
white arrow shows the FVA dip in the bottom layer.
algorithm to map the model parameter PPD. This method- References
ology has permitted full assessment of model parameter Anglin, D. K., and M. J. Fouch (2005), Seismic anisotropy in the Izu-Bonin
uncertainties using the 1- and 2-D PPD marginals. Finally, subduction system, Geophys. Res. Lett., 32, 1–4, L09307, doi:10.1029/
2005GL022714.
the robustness of our results is largely dependent upon the Babuska, V., and J. Plomerova (1993), Lithospheric thickness and velocity
well-resolved anisotropic signals that the stacking of anisotropy; Seismological and geothermal aspects, Tectonophysics, 225,
30 broadband stations from the Billings array has permitted. 79–89.
Bank, C., and M. G. Bostock (2003), Linearized inverse scattering of tele-
seismic waves for anisotropic crust and mantle structure: 2. Numerical
[38] Acknowledgments. We thank the Billings array deployment examples and application to data from Canadian stations, J. Geophys.
group, the IRIS PASSCAL data center and the PASSCAL instrument Res., 108(B5), 2259, doi:10.1029/2002JB001950.
center. Martha Savage, the Associate Editor and an anonymous reviewer Behn, M. D., C. P. Conrad, and P. G. Silver (2004), Detection of upper
provide detailed comments that improved the manuscript. We are also mantle flow associated with the African Superplume, Earth Planet. Sci.
grateful to M. Sambridge and A. Frederiksen for providing codes in original Lett., 224, 259–274.
format, and S. Wulff and J. Pierce for useful discussions. This project was Bokelmann, G. H. R. (2002), What forces drive North America?, Geology,
funded by NSF Geophysics program grant EAR-0440432. 30, 1027–1030.
19 of 20
B03304 YUAN ET AL.: FIVE UPPER MANTLE ANISOTROPY MODELS B03304
Carlson, R. W., A. J. Irving, D. J. Schulze, and B. C. Hearn (2004), Timing Levin, V., W. Menke, and J. Park (1999), Shear wave splitting in the
of precambrian melt depletion and phanerozoic refertilization events in Appalachians and the Urals; A case for multilayered anisotropy, J. Geo-
the lithospheric mantle of the Wyoming Craton and adjacent Central phys. Res., 104, 17,975–917,994.
Plains Orogen, Lithos, 77, 453–472. Levin, V., J. Park, M. Brandon, J. Lees, V. Peyton, E. Gordeev, and A.
Chen, C.-W., S. Rondenay, D. S. Weeraratne, and D. B. Snyder (2007), Ozerov (2002), Crust and upper mantle of Kamchatka from teleseismic
New constraints on the upper mantle structure of the slave craton from receiver functions, in Structure of the continental lithosphere and upper
Rayleigh wave inversion, Geophys. Res. Lett., 34, L10301, doi:10.1029/ mantle, edited by I. M. Artemieva et al., Elsevier, Amsterdam.
2007GL029535. Levin, V., A. Henza, and J. Park (2006), Texture of mantle lithosphere
Chevrot, S., and R. D. van der Hilst (2003), On the effects of a dipping axis along the Dead Sea Rift: Recently imposed or inherited?, Phys. Earth
of symmetry on shear wave splitting measurements in a transversely Planet. Int., 158, 174–189.
isotropic medium, Geophys. J. Int., 152, 497–505. Long, M. D., and R. D. van der Hilst (2005), Upper mantle anisotropy
Crampin, S. (1977), A review of the effects of anisotropic layering on the beneath Japan from shear wave splitting, Phys. Earth Planet. Inter., 151,
propagation of seismic waves, Geophys. J. R. Astron. Soc., 49, 9–27. 206–222.
Currie, C. A., J. F. Cassidy, R. D. Hyndman, and M. G. Bostock (2004), Marone, F., and B. Romanowicz (2007), The depth distribution of azimuthal
Shear wave anisotropy beneath the Cascadia subduction zone and wes- anisotropy in the continental upper mantle, Nature, 447, 198–201.
tern North American Craton, Geophys. J. Int., 157, 341–353. Menke, W., and V. Levin (2003), The cross-convolution method for inter-
Davis, P. M. (2003), Azimuthal variation in seismic anisotropy of the south- preting SKS splitting observations, with application to one and two-layer
ern California uppermost mantle, J. Geophys., 109. anisotropic earth models, Geophys. J. Int., 154, 379–392.
Deschamps, F., S. Lebedev, T. Meier, and J. Trampert (2006), Stratification Ozalaybey, S., and M. K. Savage (1994), Double-layer anisotropy resolved
of seismic anisotropy beneath the east-central United States, EOS Trans. from S phases, Geophys. J. Int., 117, 653–664.
AGU, 87(52), Fall Meet. Suppl. Abstr. T53C-1630. Park, J., and V. Levin (2002), Seismic anisotropy; Tracing plate dynamics
Diebold, J. B. (1987), Three-dimensional traveltime equation for dipping in the mantle, Science, 296, 485–489.
layers, Geophysics, 52, 1492–1500. Rumpker, G., and P. G. Silver (2000), Calculating splitting parameters for
Efron, B., and R. Tibshitani (1986), Bootstrap methods for standard errors, plume-type anisotropic structures of the upper mantle, Geophys. J. Int.,
confidence intervals, and other measures of statistical accuracy, Stat. Sci., 143, 507–520.
1, 54–77. Sambridge, M. (1999a), Geophysical Inversion with a Neighbourhood Algo-
Erslev, E. A., and J. F. Sutter (1990), Evidence for Proterozoic mylonitiza- rithm -I. Searching a parameter space, Geophys. J. Int., 138, 479–494.
tion in the northwestern Wyoming province, Geol. Soc. Am. Bull., 102, Sambridge, M. (1999b), Geophysical Inversion with a Neighbourhood Algo-
1681–1694. rithm -II. Appraising the ensemble, Geophys. J. Int., 138, 727–746.
Farra, V., L. Vinnik, B. Romanowicz, G. Kosarev, and R. Kind (1991), Sambridge, M. (2001), Finding acceptable models in nonlinear inverse
Inversion of teleseismic S particle motion for azimuthal anisotropy in the problems using a neighbourhood algorithm, Inverse Problems, 17,
upper mantle: A feasibility study, Geophys. J. Int., 106, 421–431. 387–403.
Fischer, K. M., and X. Yang (1994), Anisotropy in Kuril-Kamchatka sub- Savage, M. K. (1999), Seismic anisotropy and mantle deformation; What
ducion zone structure, Geophys. Res. Lett., 21(1), 5–8. have we learned from shear wave splitting?, Rev. Geophys., 37, 65–106.
Fischer, K. M., E. M. Parmentier, A. R. Stine, and E. R. Wolf (2000), Savage, M. K., P. G. Silver, and R. P. Meyer (1990), Observations of
Modeling anisotropy and plate-driven flow in the Tonga subduction zone teleseismic shear-wave splitting in the Basin and Range from portable
back arc, J. Geophys. Res., 105, 16,181–116,191. and permanent stations, Geophys. Res. Lett., 17(1), 21–24.
Fischer, K. M., A. Li, D. W. Forsyth, and S. H. Hung (2005), Imaging Schutt, D. L., K. Dueker, and H. Yuan (2008), Crust and upper mantle
three-dimensional anisotropy with broadband seismometer arrays, in velocity structure of the Yellowstone hotspot and surroundings, J. Geo-
Seismic Data Analysis and Imaging With Global and Local Arrays, edi- phys. Res., doi:10.1029/2007JB005109, in press.
ted by A. R. Levander and G. Nolet, American Geophysical Union, Sherrington, H. F., G. Zandt, and A. W. Frederiksen (2004), Crustal fabric
Washington, DC. in the Tibetan Plateau based on waveform inversions for seismic aniso-
Fouch, M., and S. Rondenay (2006), Seismic anisotropy beneath stable tropy parameters, J. Geophys. Res., 109, B02312, doi:10.1029/
continental interiors, Phys. Earth Planet. Inter., 158, 292–320. 2002JB002345.
Fox, O., and A. F. Sheehan (2005), Shear wave splitting beneath the Silver, P. G., and W. W. Chan (1988), Implications for continental structure
CDROM transects, in The Rocky Mountain region–an evolving litho- and evolution from seismic anisotropy, Nature, 335, 34–39.
sphere: tectonics, geochemistry, and geophysics, edited by G. Randy Silver, P. G., and W. W. Chan (1991), Shear wave splitting and subconti-
and K. E. Karlstrom, American Geophysical Union, Washington D. C. nental mantle deformation, J. Geophys. Res., 96, 16,429–416,454.
Frederiksen, A. W., and M. G. Bostock (2000), Modelling teleseismic Silver, P. G., and M. K. Savage (1994), The interpretation of shear-wave
waves in dipping anisotropic structures, Geophys. J. Int., 141, 401–412. splitting parameters in the presence of two anisotropic layers, Geophys. J.
Frederiksen, A. W., H. Folsom, and G. Zandt (2003), Neighbourhood in- Int., 119, 949–963.
version of teleseismic Ps conversions for anisotropy and layer dip, Geo- Snyder, D. B., and M. Bruneton (2007), Seismic anisotropy of the Slave
phy. J. Int., 155, 200–212. craton, NW Canada, from joint interpretation of SKS and Rayleigh
Gledhill, K., and D. Gubbins (1996), SKS splitting and the seismic aniso- waves, Geophys J. Int., 169, 170–188.
tropy of the mantle beneath the Hikurangi subduction zone, New Zeal- Snyder, D. B., M. G. Bostock, and G. D. Lockhart (2003), Two anisotropic
and, Phys. Earth Planet. Inter., 95, 227–236. layers in the Slave craton, Lithos, 71.
Gorman, A. R., et al. (2002), Deep probe: Imaging the roots of western van Decar, J. C., and R. S. Crosson (1990), Determination of teleseismic
North America, Can. J. Earth Sci., 39, 375–398. relative phase velocity arrival times using multi-channel cross-correlation
Gripp, A. E., and R. G. Gordon (2002), Young tracks of hotspots and and least squares, Bull. Seismol. Soc. Am., 80, 150–169.
current plate velocities, Geophys. J. Int., 150(2), 321. Vinnik, L. P., V. Farra, and B. Romanowicz (1989), Azimuthal anisotropy
Hartog, R., and S. Y. Schwartz (2000), Subduction-induced strain in the in the Earth from observations of SKS at GEOSCOPE and NARS broad-
upper mantle east of the Mendocino triple junction, California, J. Geo- band stations, Bull. Seismol. Soc. Am., 79, 1542–1558.
phys. Res., 105, 7909–7930. Waite, G., D. L. Schutt, and R. B. Smith (2005), Models of lithosphere and
Henstock, T. J., et al. (1998), Probing the Archean and Proterozoic litho- asthenosphere anisotropic structure of the Yellowstone hotspot from shear
sphere of western North America, GSA Today, 8, 1–5. wave splitting, J. Geophys. Res., 110, B11304, doi:10.1029/2004JB003501.
Hoffman, P. F. (1989), Precambrian geology and tectonic history of North Walker, K. T., G. H. R. Bokelmann, and S. L. Klemperer (2004), Shear-
America, in The Geology of North America - An Overview, edited by wave splitting beneath the Snake River Plain suggests a mantle upwelling
A. W. Bally and A. R. Palmer, pp. 447–512, Geology Society of Amer- beneath eastern Nevada, USA, Earth Planet. Sci. Lett., 222, 529–542.
ica, Boulder, CO. Yang, X., and K. M. Fischer (1994), Constraints on North Atlantic upper
Holtzman, B. K., T. Hiraga, J. Hustoft, D. L. Kohlstedt, M. E. Zimmerman, mantle anisotropy from S and SS phases, Geophys. Res. Lett., 21(4),
and F. Heidelbach (2003), Melt segregation and strain partitioning: Impli- 309–312.
cations for seismic anisotropy and mantle flow, Science, 301, 1227–1230. Yuan, H., and K. Dueker (2005), P-wave Tomogram of the Yellowstone
Jung, H., and S. Karato (2001), Water-induced fabric transitions in olivine, plume, Geophys. Res. Lett., 32, L07304, doi:10.1029/2004GL022056.
Science, 293, 1460–1462. Yuan, H., K. Dueker, and D. Schutt (2006), Synoptic scale crustal thickness
Kaminski, E., and N. M. Ribe (2002), Timescales for the evolution of and velocity maps along the Yellowstone hotspot track, Eos Trans. AGU,
seismic anisotropy in mantle flow, American Geophysical Union and 87(52), Fall Meet. Suppl., Abstract S43A-1376.
The Geochemical Society, United States. Pages: 17.
Katayame, I., H. Jung, and S. Karato (2004), New type of olivine fabric
from deformation experiments at modest water content and low stress, K. Dueker, D. L. Schutt, and H. Yuan, Department of Geology and
Geology, 32, 1045–1048. Geophysics, University of Wyoming, USA. (huaiyu.yuan@berkeley.edu)
20 of 20