Highresolution imaging of the deep anisotropic structure of the San ...

Geophysical Journal International Geophys. J. Int. (2011) 186, 418–446

doi: 10.1111/j.1365-246X.2011.05082.x

High-resolution imaging of the deep anisotropic structure of the San Andreas Fault system beneath southern California GJI Geodynamics and tectonics

Vadim Monteiller and Sébastien Chevrot Dynamique Terrestre et Planétaire, Observatoire Midi-Pyrénées, Université Paul Sabatier, Toulouse 14 Avenue Edouard Belin, 31400 Toulouse, France. E-mail: [email protected]

Accepted 2011 May 16. Received 2011 May 16; in original form 2011 March 10

SUMMARY Geological and geodetic records show that motion is distributed over a broad and complex system of faults in southern California, with many large earthquakes occurring off the main San Andreas fault. This unusually complex plate boundary geometry is imposed by the big bend of the San Andreas Fault and the dynamics of the underlying lithosphere, but lithosphere rheology and distribution of applied forces are still poorly constrained. Geodynamic modelling is thus very uncertain, and distribution of lithospheric deformations highly controversial. Here, exploiting the property that deformation orients olivine crystals, we show that deep lithospheric deformations can be mapped directly with a new imaging approach of seismic anisotropy, with a resolution fine enough for detailed geodynamic interpretation. This method relies entirely on finite-frequency effects in the splitting of SKS waves. We build an extensive data set of more than 3400 SKS splitting measurements from the analysis of SKS waves recorded by all the broad-band stations in southern California. Resolution tests demonstrate that, using this data set, we are able to resolve anisotropic structures smaller than the size of the first Fresnel zone of SKS waves. The good vertical resolution in our tomographic images results from the short inter-station spacing in southern California. In our 3-D anisotropy model, we find no evidence for localized lithospheric shear deformation beneath the San Andreas Fault. Instead, the lithospheric plate boundary is localized along a broad shear zone beneath the East California Shear Zone. Therefore, surface and deep deformation patterns are poorly correlated and most likely decoupled. Active N–S convergence in the Transverse Ranges results in strongly coherent E–W alignment of olivine fast axes in the shallow lithosphere. At the top of the asthenosphere, the observed low level of anisotropy suggests weak lithospheric basal tractions and thus a strong decoupling between the lithosphere and the asthenosphere. Keywords: Seismic anisotropy; Seismic tomography; Transform faults; Dynamics of lithosphere and mantle; Rheology: crust and lithosphere.

1 I N T RO D U C T I O N In southern California, about 70 per cent of the total 50 mm yr−1 relative motion between the Pacific and North America plates is accommodated by strike-slip motion on the San Andreas Fault (SAF) in Central California (deMets & Dixon 1999; Becker et al. 2005; Meade & Hager 2005). North of the San Bernardino Mountains, the Mojave segment of the SAF accommodates 9–16 mm yr−1 whereas the East California Shear Zone (Dokka & Travis 1990) accommodates 15–18 mm yr−1 (Becker et al. 2005). Further south, motion is distributed between the Indio segment of the SAF (∼23 mm yr−1 ) and the San Jacinto Fault (∼15 mm yr−1 ) (Fig. 1). In southern California, it thus appears that the different segments of the SAF have very different slip rates, and that deformation is distributed over a broad system of faults subparallel to the San Andreas Fault. The East California Shear Zone and the Big Bend, a left step in the

418

SAF, initiated with the opening of the Gulf of California, around 5–12 Ma (Liu et al. 2010). The Big Bend results in shortening and uplift of the Transverse Ranges. This complex geometry should promote slip on offshore faults and on the East California Shear Zone (Li & Liu 2006). Downwelling of the lithosphere beneath the Transverse Ranges, evidenced by seismic tomography (Humphreys & Clayton 1990; Kohler et al. 2003; Yang & Forsyth 2006; Tian et al. 2007), is thought to have maintained the SAF as the dominant fault of the system (Humphreys & Hager 1990; Li & Liu 2006; Fay et al. 2008). Understanding crustal tectonics and fault mechanics thus requires characterizing the dynamics of the lithosphere and its interaction with surface processes. In principle, lithospheric deformations can be modelled if both applied forces and rheology of the lithosphere are known. However, in practice, the strong spatial variability of lithospheric strength makes this modelling extremely difficult, especially in the vicinity of plate boundaries (Conrad et al. C 2011 The Authors C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California

419

Figure 1. Topographic relief and main active faults in southern California. ECSZ, East California Shear Zone; WTR, Western Transverse Ranges; CTR, Central Transverse Ranges; ETR, Eastern Transverse Ranges; ST, Salton Trough.

2007). Consequently, two simplified approaches are generally applied to model lithospheric deformations (Thatcher 1995). The first emphasizes the role of faults and discontinuous deformation localized along the boundaries of well-defined lithospheric blocks. In the other, strain is continuously and broadly distributed over the deforming region. Which of these two end-member models more accurately describes deformation of the continental lithosphere remains an open problem. Above ∼300-km depth, deformation occurs as dislocation creep, and aligns olivine crystals to produce macroscopic anisotropy. In principle, 3-D imaging of seismic anisotropy could thus provide crucial constraints on the distribution and localization of lithospheric deformation and mantle dynamics. This has not been done yet, owing to insufficient spatial resolution of classical tomographic methods. However, theoretical and methodological developments have recently opened new perspectives for improved resolution (Chevrot 2006; Chevrot & Monteiller 2009). We apply here for the first time this new imaging approach on southern California data. With its particularly dense station coverage, especially in the Los Angeles region, southern California provides one of the most favourable settings for a first attempt to image seismic anisotropy with a fine resolution. After a short overview of the theory, which relies entirely on finite-frequency effects in SKS splitting, we describe the method to measure splitting intensity from radial and transverse SKS records. We then perform different synthetic experiments to evaluate the resolution in our anisotropy models. Finally, we explore in detail all the geological and geodynamical implications of our final 3-D anisotropy model. 2 O U T L I N E O F T H E T H E O RY Our tomographic approach relies on shear wave splitting, which is usually quantified by measuring apparent fast directions and splitting delay times on individual SKS or SKKS records. The difficulty to relate these parameters to anisotropic perturbations of the elasticity tensor has hampered for a long time the formulation of a tomographic problem for imaging upper mantle anisotropy. This obstacle has been crossed with the introduction of a new seismological observable, the so-called splitting intensity (SI). It is obtained by C

2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

projecting the derivative of the radial component on the transverse component (Chevrot 2000) T (t)R (t) dt . (1) S = −2 R 2(t) dt An equivalent definition of SI has been formulated recently by Sieminski et al. (2008). As first noted by Chevrot (2006), the similarity of eq. (1) with the definition of a traveltime delay given by Dahlen et al. (2000) suggests that one can derive 3-D kernels for splitting intensity by plugging into this equation the Born expression of the perturbation of displacement produced by anisotropic perturbations of the elasticity tensor. Because a general anisotropic elasticity tensor is described by 21 coefficients, it is necessary to introduce a simplified (but accurate) description of upper mantle anisotropy. Such a description has been proposed by Becker et al. (2006a), to which we refer the reader for more detail. The main conclusions of this study were that (1) upper-mantle anisotropy is to first-order hexagonal and that (2) the hexagonal parameters (describing P-wave anisotropy), γ (describing S-wave anisotropy) and δ (an extra parameter describing the shape of P- and S-wave slowness surfaces), are strongly correlated. This allows us to considerably simplify the description of upper-mantle anisotropy which then requires only three parameters: one parameter, say γ , to quantify the amplitude of shear anisotropy, and two angles that define the orientation of the symmetry axis. Following Montagner & Nataf (1988), a tomographic problem can thus be formulated, which consists in retrieving the spatial variations of an anisotropy vector, defined by its length (representing the anisotropy parameter γ ) and its orientation. The description of anisotropy can be further simplified by assuming that the symmetry axis is horizontal. This is a safe assumption because for subvertically propagating SKS or SKKS waves, the signature of the dip of the symmetry axis on splitting measurements is very small (Chevrot & van der Hilst 2003). In this case, as already shown in Favier & Chevrot (2003), splitting intensity is related to only two parameters: γ c = γ cos 2φ and γ s = γ sin 2φ, where φ is the azimuth of the projection of the symmetry axis on the horizontal plane. Note

420

V. Monteiller and S. Chevrot

that these two parameters are identical to the parameters G c and G s , introduced by Montagner & Griot-Pommera (2000), following a completely different line of thought, to relate body wave and surface wave anisotropy. At this point, the computation of Fréchet kernels for splitting intensity is a well-posed problem. It can be solved with different approaches, involving different levels of approximations. In their pioneering study, Favier & Chevrot (2003) computed SI kernels for vertically propagating SKS waves in a homogeneous background

medium. They later found that it is necessary to account for the near field contribution of Green’s functions to get correct sensitivity values for depths smaller than two wavelengths (Favier et al. 2004). These results were then generalized to account for the nonvertical incidence of SKS waves in Chevrot (2006). A comparison with the results of full wave modelling with the spectral element method has shown that this approach gives very accurate results. We thus adopt it to compute all the 3-D finite-frequency kernels that are used in this study. In our computations, because we

Figure 2. Sensitivity kernels for γ c (left panel) and γ s (right panel) for an SKS wave coming from an azimuth of 22.5◦ . The wavelet is a 15 s dominant period second derivative of a Gaussian. Top panel: cross sections along the E–W vertical plane. Bottom panel: horizontal cross sections at 150-km depth. C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California assume an homogeneous background medium, we thus neglect the curvature of SKS rays between 30- and 300-km depth. This introduces a modelling error in the inversion, but this error is actually small and negligible compared to the errors on SI measurements (discussed in the next section). Exact sensitivity kernels in model PREM have been obtained with a spectral element method by Sieminski et al. (2007), but this required high-end computational fa-

421

cilities to get a single kernel. In contrast, our approach requires only very modest computational resources and can be actually performed on a laptop. Therefore, while possible, we think that the quality of the data does not warrant spending heavy computational resources to get more accurate sensitivity kernels. In addition, large-scale problems such as the one we considered in this study is probably still beyond the reach of modern supercomputers. More efficient

Figure 3. Examples of splitting intensity measurements on three SKS waves recorded by station PFO. Wiener-filtered radial and transverse components are traced with solid lines, and radial components derivatives with dashed lines. C


422


methods to compute exact and complete kernels in spherically symmetric models are now being developed (Zhao & Chevrot 2011a,b), opening new perspectives for their utilization in high-resolution imaging of the mantle in the near future. In classical regional ray tomography, resolution results from ray crossing. In regional finite-frequency tomography, resolution results from the partial overlap of Fresnel zones at adjacent stations. Fig. 2 shows γ c and γ s kernels for an SKS wave with a dominant period of 15 s coming from a backazimuth of 22.5◦ . Because the width of the first Fresnel zone is about 100 km at 100-km depth, Fresnel zones at adjacent stations will have a significant overlap if they are separated by less than about 50 km. Therefore, contrarily to common belief, it is possible to obtain images of upper mantle anisotropy with a fine lateral and vertical resolution, even if we consider vertically propagating SKS waves. The only requirement is a small station spacing. This approach, which relies entirely on finite-frequency effects, can thus be only applied in regions densely covered by broad-band seismological sensors. With its particularly dense station coverage, especially in the Los Angeles region, southern California provides one of the most favourable settings to test this new imaging approach.

3 D ATA A N A LY S I S

Amplitude of averaged splitting intensity

We measured splitting intensity of SKS and SKKS waves from large (magnitude ≥5.8) teleseismic earthquakes in the distance range 85◦ –120◦ for SKS and 120◦ –145◦ for SKKS, recorded by all the broad-band stations in southern California during the time period 2000–2009. Individual measurements of SI have large uncertainties because the measurement process involves comparing radial and transverse component waveforms which are both contaminated by noise. From a detailed statistical analysis of SI measurements obtained on a set of permanent broad-band stations distributed worldwide, various improvements were recently proposed by Monteiller & Chevrot (2010) to obtain more robust measurements. We thus carefully followed their approach, which is described in the following. A Wiener filter is first applied to each radial and transverse component records so that the radial components are standardized to a

Station PFO, Φ = 78°

Ricker (second derivative of a Gaussian) wavelet with a prescribed dominant period (here 15 or 20 s). We thus exploit SKS signals in the period range between approximately 10 and 30 s. We found that at shorter period, measurements become less reliable, owing to microseismic noise, and that at longer period the energy on SKS transverse components becomes very weak, and the records difficult to exploit. Wiener filters allow us to fine tune the frequency content of seismic waves, and thus probe anisotropy at different scales. In addition, choosing a reference Ricker wavelet leads to semi-analytic formula for the sensitivity kernels (Favier & Chevrot 2003), which makes their computation very efficient. After filtering, all the records are standardized, and SI measurements can be averaged over 10◦ azimuthal windows, which has the effect of reducing measurement errors. Because these errors scale as N −1/2 with N the number of measurements that are averaged, it is essential to use the most comprehensive data set possible for each station. After analysing 107 earthquakes recorded by more than 200 permanent broad-band stations installed in southern California, we obtained 1734 bin-averaged splitting intensities inside each frequency band, or a total of 3468 bin-averaged SI measurements. Fig. 3 illustrates our procedure to measure splitting intensities on three SKS waves recorded by station PFO. The radial components are basically identical after their standardization to a 15 s dominant period Ricker wavelet. Note also that the temporal derivatives of the radial components (dashed lines) are very similar to the transverse components (solid lines). The splitting intensity simply quantifies the amplitude ratio between the two. Fig. 4 summarizes all our SI measurements as a function of backazimuth at the same station PFO. The azimuthal coverage at this station is very good, and representative of the data set that can be obtained on permanent stations that have been operating for a long time. For stations installed during the last few years, the azimuthal coverage is usually degraded. Apparent splitting delays and fast directions can be determined by fitting a 180◦ periodic sinusoid through splitting intensity measurements at each individual station. The amplitude of the sinusoid gives the splitting delay δt and its phase the fast direction φ (Chevrot 2000). These splitting parameters represent the average or apparent anisotropy beneath a particular seismological station. Table A1 and Fig. B1 compile the SI measurements and apparent splitting parameters that we determined at each station. For station PFO, we obtain

δ

1.5 1 0.5 0

0

45

90

135 180 225 270 Azimuth (in degrees)

315

360

Figure 4. Variations of splitting intensity with backazimuth at station PFO. Vertical error bars indicate the 95 per cent confidence level (2 SD) of azimuthally averaged splitting intensity measurements. C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


423

Figure 5. (Top panel) Map of apparent splitting delays and fast directions plotted on station locations. The orientation and length of black bars respectively indicate fast direction and splitting delay. The colour background shows the interpolated splitting delays, which range from 0.3 to 2.0 s. The values of splitting parameters at each station are given in Table A1. SAF, San Andreas fault; GF, Garlock fault; ECSZ, East California Shear Zone; WLB, Walker Lane Belt; WTR, western Transverse Ranges; CTR, Central Transverse Ranges; ETR, eastern Transverse Ranges; ST, Salton Trough. (Bottom panel) Map of splitting parameters predicted by our 3-D anisotropy model.

C


424


Figure 6. Results of the inversion on a single-layer checkerboard model between 30- and 150-km depth, obtained after solving a linear inverse problem to find γ c and γ s . The checkerboard pattern has alternating blocks with φ = 45◦ and γ = −0.04 (4 per cent anisotropy), and φ = 90◦ and γ = −0.02 (2 per cent anisotropy). Note that the fast directions are correctly retrieved but that the amplitudes of shear anisotropy are strongly biased.

a fast direction of 78 ± 2◦ and a splitting delay of 0.98 ± 0.10 s. We can compare these values to those obtained in older studies on much smaller data sets. For example, from the analysis of 18 events, Liu et al. (1995) obtained φ = 92 ± 0.2◦ and δt = 1.2 ± 0.1 ¨ s, whereas Ozalaybey & Savage (1995) obtained φ = 89 ± 2◦ and δt = 1.45 ± 0.1 s. Two other studies (Helffrich et al. 1994; Polet & Kanamori 2002) report inconsistent measurements, with apparent splitting parameters varying with the polarization of the incoming wave. The fast direction that we determined at station PFO is thus in good agreement with previous studies, but the splitting delay is significantly smaller. This discrepancy most likely results from the utilization of the Silver & Chan (1991) method in the older studies, which tends to overestimate splitting delays on noisy records, as demonstrated in Monteiller & Chevrot (2010). The solution to this problem is known and consists in stacking error maps obtained for different events (and not averaging the individual measurements, which results in important biases) as proposed by Wolfe & Silver (1998), but for obscure reasons this method does not seem to be very popular and is seldom used. Overall, we are confident that our

splitting measurements currently represent the most complete and reliable data set for studying anisotropy in southern California. The compilation of the splitting parameters for all the stations in southern California is given in Table A1 and a map of these splitting parameters is shown in Fig. 5. On this map, we have represented the apparent splitting parameters at each individual station by arrows oriented along the fast directions, with their lengths proportional to the splitting delays. The colour background shows the interpolated splitting delays, which vary from 0.3 to 2.0 s. The dominant feature in this map is the preferential E–W alignment of fast directions, which has been attributed to asthenospheric flow in the slab window left behind the subducted Farallon plate (Liu et al. ¨ 1995; Ozalaybey & Savage 1995). Geodynamical modelling later confirmed this view and showed that a deep eastward flow is indeed able to predict observations of shear wave splitting (Becker et al. 2006b). However, owing to sparse data coverage and lack of vertical resolution, splitting measurements alone could not discriminate between models with a strong plate-mantle coupling (Becker et al. 2006b) from those with a well-developed asthenospheric C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California decoupling (Silver & Holt 2002). The absence of rotation of the fast directions in the vicinity of the San Andreas Fault suggests that the plate boundary has no signature and thus little influence in shaping the deformations in the upper mantle. Overall, our map is in excellent agreement with previously published studies of shear wave ¨ splitting in southern California (e.g. Liu et al. 1995; Ozalaybey & Savage 1995; Polet & Kanamori 2002; Becker et al. 2006b). However, a more careful analysis reveals the presence of small-scale variations in both fast directions and splitting delays, which are clearly seen here owing to the improved spatial coverage of our data set. These variations, of a few tens of degrees for the fast directions and of more than 1 s for splitting delays, are well above our measurement uncertainties. We thus believe that they are real and that they require complex and heterogeneous anisotropic structures beneath southern California. We will now try to image them in detail, using the high-resolution finite-frequency tomographic approach described above.

425

4 INVERSION METHOD The tomographic problem that we want to solve consists in finding the 3-D distribution of an anisotropy vector, defined by its projections γ c and γ s on the axes of a geographical reference frame. Note that it is very similar to the classical problem of azimuthal anisotropy imaging from surface waves (e.g. Montagner & Tanimoto 1991). It is linear, but relies on the implicit assumption that γ c and γ s follow a Gaussian distribution which may introduce biases in the results of the inversion (Chevrot & Monteiller 2009). An alternative is to parameterize the anisotropy vector by its length and its orientation angle, which leads to a non-linear inverse problem (Chevrot 2006). In principle, this problem can be solved with a Gauss–Newton algorithm, but it was found that, owing to the strong non-linearity of the problem, this approach does not work and always gets trapped inside a local minimum (Chevrot & Monteiller 2009). To overcome these difficulties, a new algorithm was proposed, which consists in separating the inversion in two steps. In the first step, the fast

Figure 7. Models of φ at 66-km depth for the same single layer checkerboard test as in Fig. 6 for different values of damping and smoothing parameters. C


426


directions are determined with a Gauss–Newton algorithm by imposing a constant but small value of γ everywhere in the model. The solution provides a small variance reduction, but synthetic tests have shown that fast directions are reasonably well recovered after a few iterations. Then, using these fast directions, the amplitudes of shear wave anisotropy are determined by formulating and solving a new linear inverse problem for parameter γ (Chevrot & Monteiller 2009). Although the second step does not raise any particular problem, different approaches can be followed to invert for the fast directions in the first step. We found that it is possible and actually much more efficient to obtain fast directions by solving the linear inverse problem: ⎛ ⎜ ⎜ ⎝

−1/2

Cd

G

μ1 μ2 S

⎞

⎛

⎟ ⎜ ⎟·m=⎜ ⎠ ⎝

−1/2

Cd

0

d

⎞ ⎟ ⎟, ⎠

0

(2)

where G is the sensitivity matrix whose lines contain the sensitivity kernels for γ c and γ s , S is a horizontal Laplacian operator which constrains the final model to be smooth, μ1 and μ2 are, respectively, damping and smoothing parameters, and m the model vector containing the spatial variations of γ c and γ s . The model is discretized inside a regular 3-D cartesian grid with 12 km cubic blocks inside which parameters γ c and γ s are assumed constant. From the solution of (2), the fast directions are given by φ = 0.5 tan−1 (γ s /γ c ). The critical point is to choose properly the values of the regularization parameters. This will be described in detail in the next section and it will suffice here to say that we use regularization parameters that are sufficiently large to suppress the rapid oscillations of φ, but sufficiently small to still get a good data fit. The resulting model, owing to its unrealistically strong fluctuations of γ , would usually be rejected, but nevertheless after performing extensive synthetic tests we have found that the fast directions so retrieved are usually correct. This linear approach is much more efficient because it avoids the time-consuming recomputation of all the sensitivity

Figure 8. Same as Fig. 7 but at 114-km depth. C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California kernels for the orientation of the fast directions at each iteration in the Gauss–Newton algorithm. Because resolving the tomographic problem (2) involves comparing the results obtained with a large number of different damping and smoothing parameters, the practical implementation of the new imaging algorithm would not have been possible otherwise. The linear inverse problem that needs to be solved in the second step is similar to (2) but with a vector model that contains the spatial variations of γ only. We distinguish two distinct depth domains while imposing the regularization constraints in the inversion for γ . The first interval lies between depths of 30 km, the average Moho depth (e.g. Zhu & Kanamori 2000) and 90 km, the average thickness of the lithosphere given by previous tomographic studies (e.g. Humphreys & Hager 1990; Yang & Forsyth 2006) in southern California. Mapping of the depth of the lithosphere–asthenosphere boundary from S-to-P conversions gives a shallower value than the one deduced from tomographic images, around 70 km (Li et al. 2007). However, we note that the interpretation of these detections is still controversial and

Figure 9. Same as Fig. 7 but at 174-km depth. C


427

that in many places, in particular in stable continental interiors, the depth of the discontinuity is clearly located inside the high-velocity lid. In addition, we note that Li et al. (2007) did not observe any deepening of the lithosphere–asthenosphere boundary beneath the Transverse Ranges, where a high-velocity anomaly is seen in virtually all tomographic images of southern California upper mantle. This also suggests that the detected mantle discontinuity may not be simply and unequivocally related to the lithosphere–asthenosphere boundary in this region. This lead us to favour the lithospheric thickness deduced from regional tomographic studies. In any case, we have checked that using a thinner lithospheric domain did not change our final anisotropy model significantly. In the following, we will show maps corresponding to a depth of 66 km in our models as representative of anisotropy in the lithosphere, but we could have similarly shown maps at say 50-km depth, which are almost indistinguishable. The second depth interval lies between 90- and 294-km depth, from the top to the bottom of the asthenosphere. The asthenospheric thickness is fixed rather arbitrarily, but we found

428


that increasing or decreasing it by 50 km did not change the results significantly, except by a constant scaling ratio for the amplitude of anisotropy. We introduce a priori values of γ in these two depth intervals of 4 and 2 per cent, respectively. This is simply implemented by replacing the zero coefficients in the right-hand side of eq. (2) by γ 0 = − 0.04 or γ 0 = − 0.02. At the end of the inversion, we get a 3-D model of anisotropy described in terms of φ, the orientation of the fast axis, and γ , the amplitude of shear wave anisotropy. 5 R E S O LU T I O N T E S T S We first performed a classical checkerboard test with a pattern of anisotropic alternating blocks with φ = 90◦ and γ = 2 per cent, and φ = 45◦ and γ = 4 per cent. The square blocks have a lateral size of 120 km, from 30- to 150-km depth. Theoretical splitting intensities are computed for each SKS or SKKS path in our complete data set by integrating the products of 3-D sensitivity

kernels with the distribution of γ c and γ s in the 3-D anisotropic models. Gaussian random noise is added to the theoretical values of splitting intensities, with the same standard deviation (0.3 s) as the one determined on real measurements (see Monteiller & Chevrot 2010, for more detail). Fig. 6 shows the best model that we were able to obtain after a single inversion step. The φ checkerboard pattern is rather well recovered, but the amplitudes of anisotropy are strongly biased. Small values of γ are observed where fast directions vary spatially over Fresnel volumes, that is around the edges of the blocks. To compensate for these anomalously small γ , larger values are required in the centre of the blocks. These artefacts in the γ model are a direct consequence of model parameterization and regularization (Monteiller & Chevrot 2010). Indeed, in this approach, as explained earlier, we make the implicit assumption that γ c and γ s follow a Gaussian distribution, which favours models in which these parameters are small when fast directions vary over the scale of the first Fresnel zone. Because the checkerboard has a pattern comparable in scale as the size of the first Fresnel zone, we do

Figure 10. Same as Fig. 7 but at 234-km depth. C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


429

Figure 11. Models of γ at 66-km depth for the same single layer checkerboard test as in Fig. 6 for different values of damping and smoothing parameters.

expect a strong signature of finite-frequency effects in the inversion results. Also noticeable is the vertical smearing of anisotropy below 150-km depth. As a consequence, the thickness of the anisotropic domain is overestimated, and the amplitude of the checkerboard pattern underestimated. All these artefacts result from the limited resolution potential of our data set. Although spatial coverage is excellent, with broad overlaps of Fresnel volumes in most parts of the model, the strong uncertainties of SI measurements require rather strong damping and smoothing to regularize the inversion, which strongly degrades the resolution. We believe that this is the main limitation of this approach. We next demonstrate how the two-step inversion scheme allows us to improve the results of the inversion. Figs 7–10 show the fast directions models obtained using various damping and smoothing parameters at respectively 66-, 114-, 174-, and 234-km depth. In general, fast directions are poorly retrieved with a large damping because strong fluctuations of γ are required to match the splitting C


intensities, as can be seen in Fig. 7 for example. When damping is strong, the inversion only allows small values of γ , and the fast directions are poorly retrieved. With little damping and smoothing, the models show strong and sharp fluctuations of φ. The best results are thus obtained with intermediate damping and smoothing parameters around respectively 10 and 50. This model would usually be discarded, owing to its unrealistically strong fluctuations of γ , much stronger than the ones observed in Fig. 6. However, the fast directions in this model are the closest from the input model and can be used in a new inversion to find a model for γ . Using these fast directions, we compute sensitivity kernels for γ and solve a new linear inverse problem. The resulting models obtained with different regularization parameters are shown in Figs 11–14, at the same depths as before. Using damping and smoothing parameters of respectively 10 and 50, we retrieve the checkerboard pattern in γ with a reasonable accuracy. The final model obtained with the two-step inversion algorithm is shown in Fig. 15. The γ

430


Figure 12. Same as Fig. 11 but at 114-km depth.

checkerboard pattern is now better recovered at shallow depth, but its amplitude is still underestimated beneath 100-km depth. Many artefacts present in the model obtained with the single-step algorithm (Fig. 6) have been suppressed or strongly attenuated. The checkerboard pattern is still visible below 150-km depth, but its amplitude is small. There is still some vertical smearing, but its effect is rather limited. Note that this resolution test is rather extreme, because the checkerboard pattern we attempt to recover is comparable in size to the SKS Fresnel zones. We then considered a 4 per cent anisotropy two-layer checkerboard model with 120 km blocks with φ = 90◦ and φ = 45◦ fast directions, between 30- and 294-km depth. The two layers have equal thickness. The final model obtained with the two-step inversion algorithm is shown in Fig. 16. The φ checkerboard pattern is rather well retrieved at the top and bottom of the model, which suggests that our finite-frequency approach does provide depth resolution even though it relies entirely on SKS waves with a subvertical

incidence. However, owing to limitation in resolution and imposed smoothing and damping constraints, the fast directions can deviate locally from the input directions, which introduces some biases in the reconstruction of γ . Note also that the amplitudes of anisotropy anomalies are underestimated, a classical feature also observed in checkerboard tests in most regional tomographic studies of isotropic seismic velocities. We then tested a model with a 36-km-thick shear zone with 4 per cent anisotropy and a fast direction at 135◦ going from 30 down to 150-km depth surrounded by a background medium with homogeneous E–W 2 per cent anisotropy. Fig. 17 shows the output model at four different depths. The shear zone is clearly seen in the φ model but its width and thickness are both overestimated. As a consequence, amplitude of anisotropy in the shear zone is underestimated, with values similar to the surrounding medium. Again, this suggests that the underestimated values of shear anisotropy simply result from the overestimation of the size of the shear zone with C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


431


anomalous orientations, owing to the imposed smoothing and to the limited resolution potential of the data set. Stronger anisotropy in small, anomalously oriented regions, will make their detection easier, but increase their apparent size. Thus, although the resolution of our data set is sufficient to detect a 36-km-wide shear zone, which is much smaller than the size of the Fresnel zone, it is not sufficient to retrieve its anisotropic strength correctly. In general, our resolution tests show that there is some trade-off between fast directions, thickness of anisotropic domains, resolution, and amplitudes of anisotropy. Similar trade-offs between size and strength of isotropic seismic velocity anomalies are also present in classical seismic tomography. The lower the resolution, the stronger the trade-offs. In any case, it is important to note that our data set is indeed able to image anisotropic structures smaller than the size of the Fresnel zone, which is about 150 km at 150-km depth.

C


6 THE MODEL OF SEISMIC A N I S O T R O P Y B E N E AT H S O U T H E R N CALIFORNIA We inverted the complete splitting intensity data set to determine the 3-D anisotropic structures beneath southern California, using the two-step inversion algorithm described above. We selected values of damping and smoothing parameters of respectively 10 and 50, which correspond to the values that gave the best results on all the synthetic tests we ran. The final model predicts splitting intensities with an error around 0.3 s, a value close to our estimate of measurement error, and provides a variance reduction of about 80 per cent. The predicted splitting parameters (Fig. 5 bottom) at each station are also in excellent agreement with the observed ones, but their pattern looks smoother (see also Fig. B1 for a comparison between observed and predicted bin-averaged SI at each station). The stations that show significant deviations from this map have very

432



few measurements (Table A1), and their splitting parameters are actually poorly constrained. Our model thus contains the simplest and smoothest anisotropic structures that explain the data within their error bars. The anisotropy model (Fig. 18) shows many features that can be related to structures found previously by isotropic tomographic imaging. At 66-km depth, a level representative of the rather thin southern California lithosphere, we find strong lateral variations of seismic anisotropy. The Garlock Fault separates a lithospheric block characterized by weak and spatially variable anisotropy to the north from another block with strong, coherently oriented E–W anisotropy to the south. The Garlock Fault has thus a deep signature, probably located along an inherited zone of lithospheric weakness. Strong E–W anisotropy is also observed in a lithospheric block located south of the eastern Transverse Ranges and north of the Salton Trough. It is tempting to relate these strong and penetrative fabrics to the transpression resulting from the ‘Big Bend’ in the SAF, which is expected

to generate a N–S compression that would align the olivine a axes along the E–W direction (Fig. 19). Active orogenic processes in the Transverse Ranges are thus currently creating a strong and coherent anisotropic fabric in the shallow lithospheric mantle. East of these strongly anisotropic lithospheric domains, we find evidence of a rather broad and deep lithospheric shear zone in which the fast directions rotate toward a NW–SE direction. Interestingly, this dextral shear zone, which is located beneath the trace of the dextral East California Shear Zone, coincides with the region where the Big Bend of the SAF favours release of plastic strain energy (Fay et al. 2008), and could correspond to the actual lithospheric plate boundary. All these findings suggest a rather complicated geometry of the plate boundary in the lithosphere, rather poorly correlated with its surface expression, the SAF. Thus, although block models have been successful in explaining plate kinematics, they seem to poorly represent the behaviour of the deep continental lithosphere, at least in the vicinity of major plate boundaries. North of the C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


433

Figure 15. Final anisotropy model obtained with the two-step inversion scheme on the single-layer checkerboard model.

Garlock Fault and around the trace of the SAF, anisotropy shows a very different pattern, with fast directions coherently oriented along the strike of the SAF. The simple structure of the SAF in central California thus seems to be associated to simpler and more coherent lithospheric deformations. This observation is in good agreement with previous investigations of shear wave splitting around the SAF in central and northern California (e.g. Savage & Silver 1993; Bonnin et al. 2010; Hartog & Schwartz 2001). In these older studies, the azimuthal dependence of splitting parameters measured with the Silver & Chan (1991) technique was interpreted in terms of a twolayer model of anisotropy, the upper layer being related to a coherent lithospheric shearing of the lithosphere around the SAF decoupled with the lower layer, the underlying asthenosphere. Unfortunately, owing to sparser distribution of seismological stations, this part of the model is not well resolved in our study. Low anisotropy around the California-Nevada border can be explained by extension and upwelling resulting from a combination of strong upwelling of buoyant material beneath southern Sierra Nevada and the south C


ern Walker Lane Belt, and deep downwelling beneath southeastern Nevada (Fay et al. 2008). A low anisotropy anomaly is also observed beneath the Salton Trough, from 50 km to at least 120-km depth, which coincides with the low-velocity anomaly imaged with classical tomographic approaches. Both signatures are simply explained by extension in the lithosphere, which generates a strong upwelling of anomalously hot mantle (Lachenbruch et al. 1985; Humphreys & Hager 1990). Fig. 19 summarizes our main interpretations of anisotropic structures in the lithosphere seen in the 3-D anisotropy model. At 114-km depth or in the shallow asthenosphere, a striking feature is the low anisotropy anomaly found beneath the Transverse Ranges. It coincides with the prominent high-velocity anomaly seen in all tomographic maps (Humphreys & Clayton 1990; Kohler et al. 2003; Yang & Forsyth 2006; Tian et al. 2007), which has been interpreted either as a lithospheric downwelling (Bird & Rosenstock 1984; Humphreys & Hager 1990) or as a gravitational instability of a cold lithospheric mantle (Houseman et al. 2000). In any case, these

434


Figure 16. Anisotropy model obtained with the two-step inversion algorithm on the double-layer checkerboard model. The checkerboard has 120 km blocks, with an alternating pattern of fast directions φ = 45◦ and φ = 90◦ , between 30- and 294-km depth. The two layers have equal thickness with 4 per cent anisotropy.

dynamical models would all predict predominantly vertical mantle flow beneath the Transverse Ranges, and little, if any, shear wave splitting, in very good agreement with our results. Another striking feature is the generally low level of anisotropy, much smaller than in the lithosphere or in the deep asthenosphere. This is a robust feature that we consistently observed in all the inversions that we performed. This is an important result because weak anisotropy in the shallow asthenosphere implies small tractions at the base of the lithosphere. This is in good agreement with the fact that plate deformations can be explained remarkably well by shear coupling across the plate boundary alone, with very little, if any, basal tractions (Humphreys & Coblentz 2007; Liu et al. 2010; Yang & Liu 2010). Our results thus point to a zone of decoupling in the asthenosphere beneath North America. Silver & Holt (2002) came to a similar conclusion from a comparison of the fast directions inferred from SKS splitting with the motion of the Pacific and North America plates.

At deeper levels, the anisotropy pattern becomes simple, with very low anisotropy and variable fast directions beneath the North America Plate, and larger anisotropy and more coherently oriented E–W fast directions beneath the Pacific Plate. We found that stronger anisotropy is required in the deep asthenosphere, in particular beneath the Pacific Plate. The exact depth location of this increased level of anisotropy depends on the a priori assumptions made on the average shear wave anisotropy profile. However, the results that the amplitude of anisotropy increases with depth in the asthenosphere seems robust, even though there is a clear trade-off between the absolute value of shear wave anisotropy and the assumed thickness of the asthenosphere, here arbitrarily fixed at 204 km. The E–W fast directions in the asthenosphere have been previously related to regional mantle flow produced by the Farallon slab window (e.g. Hartog & Schwartz 2001; Silver & Holt 2002; Becker et al. 2006b). However, this simple model does not explain the C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


435

Figure 17. Anisotropy model obtained with the two-step inversion algorithm on the shear zone model. In the 36-km-thick shear zone, φ = 135◦ and γ = 4 per cent anisotropy between 30- and 150-km depth. In the surrounding medium, φ = 90◦ and γ = 2 per cent.

apparent decrease of anisotropy in the deep asthenosphere toward the eastern part of the model, which is a robust observation already visible in the surface map of apparent SKS splitting (Fig. 5). It is possible that beneath North America density-driven mantle flow becomes more vertical, as in the model of Fay et al. (2008), for example, which would tend to disorganize horizontal fabrics and thus weaken SKS splitting. Another possibility is that the eastward flow is deflected downward by the thicker North American lithosphere. Alternatively, the apparent small splitting observed on the North America Plate could result from the combined effect of the contribution of a deep asthenospheric layer with predominantly E–W flow and of a lithospheric layer with fabrics coherently oriented along the N–S direction, possibly produced when the Farallon Plate was still subducting beneath the North America Plate, resulting in a predominantly E–W compression. Such a simple model would explain the small apparent splitting on the North America Plate with a coherent regional model of eastward mantle flow in the astheno C


sphere. These issues could probably be addressed in the future by detailed geodynamical modelling.

7 DISCUSSION To identify the main factor limiting resolution in our newly developed imaging approach, we have made a simple experiment in which we inverted the apparent splitting parameters given in Table A1, using 3-D kernels obtained by averaging all the kernels corresponding to each individual station used in the global inversion. We obtained a 3-D model of anisotropy very similar to the one shown in Fig. 18. This demonstrates that the main limitation factor for the resolution is not the accuracy with which we can compute each individual kernel but the strong microseismic noise that contaminates horizontal components records, which requires introducing rather strong damping and smoothing constraints to stabilize the inversion. Indeed, many permanent stations were installed

436


Figure 18. Horizontal cross-sections in the anisotropy model at 66-, 114-, 174-, and 234-km depth. The orientation and length of black bars respectively indicate fast direction and shear wave anisotropy. The colour background shows the amplitude of shear wave anisotropy.

during the last few years and for these stations we were able to obtain very few high quality SI measurements (Table A1). We thus expect that significant improvements of our high-resolution images will have to wait for the construction of a larger high-quality data set from a longer period of observation. However, this approach already offers a tremendous improvement in the resolution of our images of upper-mantle anisotropy compared to those classically obtained from azimuthal anisotropy of surface waves. Various attempts have been made recently to improve this latter approach by performing a joint inversion of apparent splitting parameters of SKS waves and phase velocities of surface waves (e.g. Yuan et al. 2011). We think that one should be cautious with these joint tomographic approaches because these two types of data have very different resolution potential. Indeed, the lateral resolution of surface waves being rather weak, these methods will tend to put the anisotropy signal of SKS waves that is discrepant from surface waves in the deepest layers. In other words, there will be

a tendency to map the strong, small-scale, variations of anisotropy localized in the lithosphere, evidenced in our study, inside the deep asthenosphere. A new method of surface wave tomography based upon the eikonal equation has been introduced recently to obtain a local measurement of the azimuthal dependence of the phase velocity of surface waves (Lin et al. 2009). With this approach, the 180◦ periodic variation of the Rayleigh waves speed was directly observed for the first time, a quite remarkable achievement. Such measurements were later used to determine a model of seismic anisotropy in the western United States (Lin et al. 2011). Because in eikonal tomography there is no explicit forward operator, resolution is difficult to determine. Rough estimates can be obtained by characterizing the correlation length of speed measurements, which gives values smaller than 70 km, the typical inter-station spacing of USArray stations. The resolution in eikonal tomography thus seems to be much better than the one in classical azimuthal anisotropy imaging, C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International


437

Figure 19. Schematic block diagram summarizing and crudely interpreting the lithospheric anisotropic structures observed in the 3-D anisotropy model shown in Fig. 18.

with the important limitation that the eikonal equation is only valid at high frequency, which means that it can probe only the very shallow part of the upper mantle. The distribution of anisotropy in the uppermost mantle in the model obtained by Lin et al. (2011) share many similarities with our model shown in Fig. 18. This encouraging observation suggests that combining SKS and surface wave data may improve our ability to image the fine anisotropic structures in the lithosphere and asthenosphere. For example, the result of eikonal tomography could be used to construct an a priori model for anisotropy in the shallow lithosphere. This would certainly improve the resolution of SI imaging in the deeper layers. Note that it would be straightforward to incorporate such an a priori constraint in our inversion algorithm. Compared to isotropic velocity tomography, vectorial tomography raises new theoretical and practical problems. The two-step inversion method is a practical way to solve inverse problems on rather large data sets, but more work is probably required to find a better parameterization of the anisotropy models and more robust inversion scheme. 8 C O N C LU S I O N S Using data from the dense array of broad-band stations in southern California, we have shown that it is possible to obtain highresolution images of upper-mantle anisotropy from SKS splitting measurements alone. Lateral and vertical resolution results from the partial overlap of Fresnel volumes at adjacent stations. The new anisotropy model reveals anisotropic structures beneath southern California with unprecedented resolution, many of which can be related to isotropic structures previously imaged by regional body and surface wave tomography. Low anisotropy is observed beneath the Transverse Ranges, where fast velocities have been interpreted as a downwelling of the lithosphere. Similarly, low anisotropy is found beneath the Salton Trough, where lithospheric extension is responsible for a strong mantle upwelling. We found no evidence of a deep C


anisotropic signature beneath the San Andreas Fault. Instead, we detect a broad shear zone beneath the East California Shear Zone, which could be related to the actual plate boundary. This suggests that surface and deep deformation processes are decoupled. The last two decades have seen the rapid development of SKS splitting analysis to infer apparent anisotropy beneath individual seismological stations. Our results open the perspective of direct 3-D imaging of seismic anisotropy from these measurements to reveal lithospheric fabrics and asthenospheric flow fields with a fine resolution, at least in regions densely equipped with seismological sensors. This holds the key for a better understanding of the rheological behaviour of the continental lithosphere and of the dynamic mantle. AC K N OW L E D G M E N T S This work was supported by the French ANR Jeunes Chercheurs program (grant ANR-JC-0014). We thank the Southern California Earthquake Data Center for providing the seismological records used in this study, and the constructive comments of two anonymous reviewers. REFERENCES Becker, T.W., Hardebeck, J.L. & Anderson, G., 2005. Constraints on fault slip rates of the southern California plate boundary from GPS velocity and stress inversions, Geophys. J. Int., 160, 634–650. Becker, T.W., Chevrot, S., Schulte-Pelkum, V. & Blackman, D.K., 2006a. Statistical properties of seismic anisotropy in the upper mantle explored by geodynamic models, J. geophys. Res., 111(B08309), doi:10.1029/2005JB004095. Becker, T.W., Schulte-Pelkum, V., Blackman, D.K., Kellog, J.B. & O’Connell, R.J., 2006b. Mantle flow under the western United States from shear wave splitting, Earth planet. Sci. Lett., 247, 235–251. Bird, P. & Rosenstock, R.W., 1984. Kinematics of present crust and mantle flow in southern California, Geol. Soc. Am. Bull., 95, 946–957.

438


Bonnin, M., Barruol, G. & Bokelmann, G.H.R., 2010. Upper mantle deformation beneath the North American-Pacific plate boundary in California from SKS splitting, J. geophys. Res., 115, B04306, doi:10.1029/2009JB006438. Chevrot, S., 2000. Multichannel analysis of shear wave splitting, J. geophys. Res., 105, 21 579–21 590. Chevrot, S., 2006. Finite frequency vectorial tomography: a new method for high resolution imaging of upper mantle anisotropy, Geophys. J. Int., 165, 641–657. Chevrot, S. & Monteiller, V., 2009. Principles of vectorial tomography – the effects of model parametrization and regularization in tomographic imaging of seismic anisotropy, Geophys. J. Int., 179, 1726–1736. Chevrot, S. & van der Hilst, R., 2003. On the effects of a dipping symmetry axis on shear wave splitting measurements in a transversely isotropic medium, Geophys. J. Int., 152, 497–505. Conrad, C.P., Behn, M.D. & Silver, P.G., 2007. Global mantle flow and the development of seismic anisotropy: differences between the oceanic and continental upper mantle, J. geophys. Res., 112, B07317, doi:10.1029/2006JB004608. Dahlen, F.A., Hung, S.H. & Nolet, G., 2000. Fréchet kernels for finitefrequency traveltimes—I. Theory, Geophys. J. Int., 141, 157–174. deMets, C. & Dixon, T.H., 1999. New kinematic models for Pacific-North America motion from 3 Ma to present. I: Evidence for steady motion and biases in the NUVEL-1A model, Geophys. Res. Lett., 26, 1921–1924. Dokka, R.K. & Travis, C.J., 1990. Role of the eastern California shear zone in accomodating Pacific-North American plate motion, Geophys. Res. Lett., 17, 1323–1326. Favier, N. & Chevrot, S., 2003. Sensitivity kernels for shear wave splitting in transverse isotropic media, Geophys. J. Int., 153, 213–228. Favier, N., Chevrot, S. & Komatitsch, D., 2004. Near-field influences on shear wave splitting and traveltime sensitivity kernels, Geophys. J. Int., 156, 467–482. Fay, N.P., Bennett, R.A., Spinler, J.C. & Humphreys, E.D., 2008. Small-scale upper mantle convection and crustal dynamics in southern California, Geochem. Geophys. Geosyst., 9, Q08006, doi:10.1029/2008GC001988. Hartog, R. & Schwartz, S.Y., 2001. Depth-dependent mantle anisotropy below the San Andreas fault system: apparent splitting parameters and waveforms, J. geophys. Res., 106, 4155–4167. Helffrich, G., Silver, P. & Given, H., 1994. Shear-wave splitting variation over short spatial scales on continents, Geophys. J. Int., 119, 561–573. Houseman, G.A., Neil, E.A. & Kohler, M.D., 2000. Lithospheric instability beneath the Transverse Ranges of California, J. geophys. Res., 105, 16 237–16 250. Humphreys, E.D. & Clayton, R.W., 1990. Tomographic image of the southern California mantle, J. geophys. Res., 95, 19 725–19 746. Humphreys, E.D. & Coblentz, D.D., 2007. North American dynamics and western U.S. tectonics, Rev. Geophys., 45, RG3001, doi:10.1029/2005RG000181. Humphreys, E.D. & Hager, B.H., 1990. A kinematic model for the late Cenozoic development of southern California crust and upper mantle, J. geophys. Res., 95, 19 747–19 762. Kohler, M.D., Magistrale, H. & Clayton, R.W., 2003. Mantle heterogeneities and the SCEC reference three-dimensional seismic velocity model version 3, Bull. seism. Soc. Am., 93, 757–774. Lachenbruch, A.H., Sass, J.H. & Galanis, S.P., 1985. Heat-flow in southernmost California and the origin of the Salton Trough, J. geophys. Res., 90, 6709–6736. Li, Q.S. & Liu, M., 2006. Geometrical impact of the San Andreas Fault on stress and seismicity in California, Geophys. Res. Lett., 33, L08302, doi:10.1029/2005GL025661. Li, X., Yuan, X. & Kind, R., 2007. The lithosphere-asthenosphere boundary beneath the western United States, Geophys. J. Int., 170, 700–710. Lin, F., Ritzwoller, M.H. & Snieder, R., 2009. Eikonal tomography: Surface wave tomography by phase front tracking across a regional broad-band seismic array, Geophys. J. Int., 177, 1091–1110. Lin, F.C., Ritzwoller, M.H., Yang, Y., Moschetti, M.P. & Fouch, M.J., 2011, Complex and variable crustal and uppermost mantle seismic anisotropy in the western United States, Nature Geoscience, 4, 55–61.

Liu, H., Davis, P.M. & Gao, S., 1995. SKS splitting beneath southern California, Geophys. Res. Lett., 22, 767–770. Liu, M., Wang, H. & Li, Q., 2010. Inception of the eastern California shear zone and evolution of the Pacific-North American plate boundary: from kinematics to geodynamics, J. geophys. Res., 115, B07401, doi:10.1029/2009JB007055. Meade, B.J. & Hager, B.H., 2005. Block models of crustal motion in southern California constrained by GPS measurements, J. geophys. Res., 110, B03403, doi:10.1029/2004JB003209. Montagner, J. & Tanimoto, T., 1991. Global upper mantle tomography of seismic velocities and anisotropies, J. geophys. Res., 96, 20 337–20 351. Montagner, J.P. & Griot-Pommera, D.A., 2000. How to relate body wave and surface wave anisotropy?, J. geophys. Res., 105, 19 015–19 027. Montagner, J.P. & Nataf, H.C., 1988. Vectorial tomography—I. Theory, Geophys. J., 94, 295–307. Monteiller, V. & Chevrot, S., 2010. How to make robust splitting measurements for single-station analysis and three-dimensional imaging of seismic anisotropy, Geophys. J. Int., 182, 311–328. ¨ Ozalaybey, S. & Savage, M.K., 1995. Shear-wave splitting beneath western United States in relation to plate-tectonics, J. geophys. Res., 100, 18 135–18 149. Polet, J. & Kanamori, H., 2002. Anisotropy beneath California: shear wave splitting measurements using a dense array, Geophys. J. Int., 149, 314–328. Savage, M.K. & Silver, P.G., 1993. Mantle deformation and tectonics: constraints from seismic anisotropy in the western United States, Phys. Earth planet. Inter., 78, 207–227. Sieminski, A., Liu, Q., Trampert, J. & Tromp, J., 2007. Finite-frequency sensitivity of body waves to anisotropy based upon adjoint methods, Geophys. J. Int., 171, 368–389. Sieminski, A., Paulssen, H., Trampert, J. & Tromp, J., 2008. Finite-frequency SKS splitting: measurement and sensitivity kernels, Bull. seism. Soc. Am., 98, 1797–1810. Silver, P.G. & Chan, W.W., 1991. Shear wave splitting and subcontinental mantle deformation, J. geophys. Res., 96, 16 429–16 454. Silver, P.G. & Holt, W.E., 2002. The mantle flow field beneath western North America, Science, 295, 1054–1057. Thatcher, W., 1995. Microplate versus continuum descriptions of active tectonic deformation, J. geophys. Res., 100, 3885–3894. Tian, Y., Zhao, D.P. & Teng, J.W., 2007. Deep structure of southern California, Phys. Earth planet. Inter., 165, 93–113. Wolfe, C.J. & Silver, P.G., 1998. Seismic anisotropy of oceanic upper mantle: shear wave splitting methodologies and observations, J. geophys. Res., 103, 749–771. Yang, Y. & Forsyth, D.W., 2006. Rayleigh wave phase velocities, smallscale convection, and azimuthal anisotropy beneath southern California, J. geophys. Res., 111, B07306, doi:10.1029/2005JB004180. Yang, Y. & Liu, M., 2010. What drives short- and long-term crustal deformation in the southwestern United States?, Geophys. Res. Lett., 37, L04306, doi:10.1029/2009GL041598. Yuan, H., Romanowicz, B., Fischer, K.M. & Abt, D., 2011. 3-d shear wave radially and azimuthally anisotropic velocity model of the North American upper mantle, Geophys. Res. Lett., 184, 12437–1260. Zhao, L. & Chevrot, S., 2011a. An efficient and flexible approach to the calculation of three-dimensional full-wave Fréchet kernels for seismic tomography—I. Theory, Geophys. J. Int., 185, 922–938. Zhao, L. & Chevrot, S., 2011b. An efficient and flexible approach to the calculation of three-dimensional full-wave Fréchet kernels for seismic tomography—II. Numerical results, Geophys. J. Int., 185, 939– 954. Zhu, L.P. & Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions, J. geophys. Res., 105, 2969–2980.

A P P E N D I X A : A P PA R E N T S P L I T T I N G PA R A M E T E R S F O R S TAT I O N S I N SOUTHERN CALIFORNIA C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California Table A1. Compilation of splitting parameters (fast directions φ, and splitting delays δt) in southern California determined by fitting a sinusoid through azimuthally averaged splitting intensity measurements at each station. φ and δt are the errors on respectively φ and δt. The last column gives the number of individual splitting intensity measurements that were stacked in 10◦ azimuthal windows.

C

Station

Latitude

Longitude

φ

φ

δt(s)

δt(s)

Nbr.

ADO AGA ALP ARV BAK BAR BBR BBS BC3 BCC BCW BDM BEL BFS BKR BLA BLY BOR BRE BTC BTP BZN CAP CAR CCC CGO CHF CHN CIA CLC CLT CMB CRP CRY CTC CVS CWC DAN DEC DEV DGR DJJ DLA DNR DPP DRC DRE DSC DVT EDW EDW2 EML ERR FIG FMP FRD FUR GMR GOR GR2

34.55046 33.63838 34.68708 35.12690 35.34444 32.68005 34.26230 33.92139 33.65515 33.57508 34.94009 37.95397 34.00060 34.23883 35.26930 34.06948 33.75034 33.26820 33.80776 33.01213 34.68224 33.49150 33.38822 35.30819 35.52495 36.55040 34.33341 33.99879 33.40190 35.81574 34.09284 38.03455 34.13624 33.56540 33.65516 38.34526 36.43988 34.63745 34.25353 33.93597 33.65001 34.10618 33.84822 33.56667 32.99862 32.80540 32.80532 35.14255 32.65915 34.88303 34.88110 32.89083 33.11675 34.72832 33.71264 33.49470 36.46703 34.78457 33.15370 34.11816

−117.43391 −116.40110 −118.29946 −118.83009 −119.10445 −116.67215 −116.92075 −116.98058 −115.45366 −117.26119 −119.41313 −121.86554 −115.99820 −117.65853 −116.07030 −116.38896 −114.52375 −116.41716 −117.98116 −115.21987 −118.57398 −116.66700 −117.19495 −119.84583 −117.36453 −117.80295 −118.02585 −117.68044 −118.41502 −117.59751 −117.31687 −120.38651 −118.12705 −116.73730 −115.99006 −122.45840 −118.08016 −115.38115 −118.33383 −116.57794 −117.00947 −118.45505 −118.09624 −116.63056 −116.94153 −115.44654 −115.44679 −116.10395 −116.10061 −117.99106 −117.99388 −116.84566 −115.82267 −119.98803 −118.29381 −116.60220 −116.86322 −115.65994 −117.22921 −118.30024

70 77 76 53 −88 86 81 88 −83 71 −83 −74 78 −89 49 72 71 86 70 −69 87 67 80 −77 72 −88 87 68 77 −83 65 83 −88 80 −83 −73 73 75 −87 90 83 83 66 74 −52 86 −77 88 43 85 83 62 76 −68 89 79 85 90 −85 81

1 1 2 1 1 1 3 1 1 1 5 1 1 2 11 2 3 2 2 3 2 2 2 3 4 2 3 2 1 2 1 2 2 1 2 2 3 2 2 2 1 1 2 2 2 9 6 3 11 1 1 2 1 2 2 2 5 1 2 3

0.62 1.25 0.87 0.60 0.85 1.00 0.66 1.39 0.82 1.25 0.80 0.90 0.86 0.90 0.65 0.90 0.53 0.53 1.27 0.93 0.80 0.84 0.90 0.61 0.30 0.45 0.69 1.00 0.90 0.61 1.01 0.69 1.60 1.20 0.89 0.97 0.69 0.66 0.63 0.69 0.91 1.10 1.07 0.80 0.67 0.27 0.55 0.47 0.39 1.66 1.15 0.47 1.16 0.88 0.41 0.68 0.37 0.72 0.82 1.72

0.03 0.05 0.04 0.07 0.06 0.03 0.03 0.04 0.04 0.04 0.19 0.03 0.04 0.05 0.20 0.08 0.06 0.05 0.07 0.11 0.06 0.06 0.06 0.05 0.04 0.04 0.04 0.06 0.05 0.04 0.05 0.04 0.06 0.03 0.07 0.05 0.07 0.05 0.05 0.05 0.03 0.04 0.07 0.04 0.03 0.22 0.19 0.07 0.02 0.19 0.03 0.04 0.01 0.05 0.05 0.05 0.05 0.02 0.05 0.37

31 18 51 33 28 47 38 36 35 15 8 26 39 39 9 40 25 60 26 17 35 23 33 18 33 31 45 22 37 52 22 49 17 45 43 36 24 43 38 40 48 45 27 43 43 6 11 39 20 11 39 43 24 36 31 45 32 21 45 7


439

440

V. Monteiller and S. Chevrot Table A1. (Continued.) Station

Latitude

Longitude

φ

φ

δt(s)

δt(s)

Nbr.

GRA GSC HEC HLL HOP IRM ISA JCC JCS JEM JRC JRC2 JRS JVA KNW LAF LCG LCP LDF LDR LFP LGB LGU LJR LKL LLS LMR LRL LTP LUG LVA MCT MGE MHC MLA MLS MON MOP MPI MPM MPP MSJ MTP MUR MWC NBS NEE NEE2 NJQ NSS OLI OLP ORV OSI PAS PASC PDE PDM PDR PDU PER PFO

36.99608 35.30177 34.82940 34.17643 38.99349 34.15738 35.66278 40.81745 33.08590 33.08098 35.98230 35.98249 37.40373 34.36622 33.71410 33.86889 34.00033 34.73551 35.13066 34.99060 34.30529 33.97530 34.10819 34.80762 34.61594 33.68447 34.93439 35.47954 33.88110 34.36560 33.35160 34.22645 33.81841 37.34164 37.63019 34.00460 32.89270 34.28085 34.81258 36.05799 34.88848 33.80801 35.48434 33.59999 34.22362 34.78035 34.82490 34.76759 34.53412 33.55553 33.94539 32.60783 39.55451 34.61450 34.14844 34.17141 34.44199 34.30336 33.96273 34.12070 33.86164 33.61170

−117.36621 −116.80574 −116.33500 −118.35967 −123.07234 −115.14513 −118.47403 −124.02955 −116.59590 −116.59755 −117.80760 −117.80885 −122.23868 −116.61266 −116.71190 −118.33143 −118.37794 −120.27996 −115.18416 −118.34156 −118.48805 −118.14918 −119.06587 −118.86775 −117.82493 −117.94304 −117.69629 −117.68212 −118.17568 −117.36683 −116.56150 −116.04073 −116.36874 −121.64257 −118.83605 −117.56162 −116.42250 −118.90490 −119.14530 −117.48901 −119.81362 −116.96789 −115.55320 −117.19543 −118.05832 −116.55798 −114.59941 −114.61881 −120.17737 −115.94586 −117.92372 −116.93036 −121.50036 −118.72350 −118.17117 −118.18523 −118.58215 −114.14152 −118.43702 −117.63808 −117.20529 −116.45940

84 79 86 72 −55 86 71 56 −75 73 −88 −83 −71 75 77 73 81 −75 53 84 48 −60 79 −71 86 82 88 65 77 −84 78 86 82 −77 58 86 −86 73 −84 88 −75 82 68 77 83 81 64 32 −64 −77 71 80 −86 −79 84 77 88 58 79 −88 −88 79

1 2 4 3 1 2 2 5 1 4 1 3 2 2 2 3 2 5 1 1 1 1 2 1 1 1 1 2 2 2 2 1 1 3 1 1 2 1 1 2 1 1 2 3 1 1 2 2 2 8 2 2 2 1 2 2 1 5 2 1 2 1

1.07 0.67 0.43 1.10 1.10 0.46 0.73 0.29 0.43 0.67 1.53 0.84 0.95 0.69 0.85 0.35 1.35 0.49 0.90 0.94 0.98 0.97 0.96 0.98 1.21 1.99 0.76 0.89 1.30 1.35 0.68 0.76 1.02 0.81 1.20 1.17 0.44 0.97 0.39 0.68 1.02 1.40 0.92 0.71 1.16 0.52 0.66 0.73 0.60 0.67 1.23 0.71 0.66 1.49 1.20 1.16 1.01 0.29 1.18 1.50 0.89 1.04

0.06 0.04 0.06 0.15 0.07 0.03 0.05 0.05 0.02 0.09 0.17 0.08 0.06 0.05 0.05 0.03 0.08 0.09 0.02 0.03 0.08 0.04 0.06 0.04 0.25 0.07 0.02 0.05 0.05 0.10 0.06 0.02 0.06 0.07 0.04 0.09 0.04 0.05 0.01 0.03 0.05 0.02 0.06 0.05 0.03 0.03 0.03 0.03 0.03 0.16 0.07 0.04 0.02 0.03 0.10 0.07 0.10 0.04 0.10 0.05 0.03 0.04

38 49 38 14 32 35 25 34 39 7 7 27 23 44 34 26 32 19 57 46 21 35 40 46 8 26 18 36 26 20 43 43 49 28 50 21 27 25 23 55 30 34 32 27 47 30 22 22 30 13 35 34 38 51 39 22 25 47 17 53 32 43

C 2011 The Authors, GJI, 186, 418–446 C 2011 RAS Geophysical Journal International

Seismic anisotropy beneath southern California Table A1. (Continued.)

C

Station

Latitude

Longitude

φ

φ

δt(s)

δt(s)

Nbr.

PHL PKD PLC PLM PLS PMD PSD PSR RCT RDM RIN RIO RPV RRX RUS RVR RXH SAL SAO SBB SBC SBI SBP SCI SCZ SDD SDG SDP SDR SES SHO SLA SLR SMB SMM SMR SMS SNC SND SOT SPF SPG SRN SSW STC STG STS SVD SWS SYP TA2 TEH TFT THX TIN TOR TOV TRO TUQ USC VCS VES VTV

35.40773 35.94517 33.82436 33.35361 33.79530 33.64785 33.82393 34.09181 36.30523 33.63000 34.28196 34.10473 33.74346 34.87533 34.05073 33.99351 33.18313 33.28010 36.76403 34.68850 34.44076 33.48046 34.23240 32.97990 33.99543 33.55259 32.78400 34.56547 32.73561 34.43692 35.89953 35.89095 33.83359 34.90218 35.31420 35.37699 34.01438 33.24787 33.55190 34.41600 34.05933 36.13550 33.82854 33.17747 34.30302 33.66400 33.79033 34.10647 32.94080 34.52775 34.38203 35.29130 35.14592 33.63495 37.05422 33.57526 34.15607 33.52340 35.43584 34.01919 34.48364 32.75904 34.56065

−120.54556 −120.54160 −116.51195 −116.86265 −117.60906 −116.37769 −116.55026 −117.80709 −119.24384 −116.84780 −118.47925 −117.97956 −118.40412 −116.99684 −118.08085 −117.37545 −115.62257 −115.98585 −121.44722 −117.82448 −119.71492 −119.02986 −117.23484 −118.54697 −119.63510 −117.66171 −117.13805 −120.50137 −116.94241 −119.13750 −116.27530 −117.28332 −116.79737 −120.44709 −119.99581 −120.61247 −118.45617 −119.52437 −116.61290 −118.44900 −118.64614 −118.81099 −117.78938 −115.61564 −119.18676 −117.76856 −118.19878 −117.09822 −115.79580 −119.97834 −117.67822 −118.42079 −119.41946 −116.16402 −118.23009 −116.22584 −118.82039 −116.42570 −115.92389 −118.28631 −118.11783 −115.73161 −117.32960

−76 −77 72 76 85 82 85 −85 80 69 63 76 83 84 88 87 −79 89 −89 −86 −75 82 −80 82 −86 81 87 −72 81 82 73 68 76 −80 −68 −80 −89 83 67 −86 −87 81 80 −78 −68 71 72 73 −76 −79 82 −88 −78 −83 60 65 81 85 70 66 87 −74 77

6 8 1 1 1 1 1 1 7 1 10 2 2 2 2 1 1 6 2 1 1 3 4 2 4 1 2 4 3 2 3 3 2 4 1 1 2 1 2 8 2 1 2 1 4 1 1 2 1 3 3 1 4 1 2 5 1 2 2 2 1 4 1

0.64 0.53 1.47 0.85 1.90 0.85 1.17 1.21 0.31 0.88 0.32 1.16 0.87 0.71 1.18 0.96 1.18 0.58 0.48 1.30 1.13 0.90 0.48 1.40 0.95 1.03 1.10 0.96 0.42 0.70 0.48 0.60 1.09 1.34 0.79 0.56 0.94 0.89 0.71 0.98 0.74 0.41 1.03 1.18 0.40 0.95 1.21 0.85 0.83 0.75 1.15 1.42 0.74 1.57 0.84 0.56 1.22 1.21 0.51 0.95 1.02 0.45 0.36

0.11 0.10 0.02 0.04 0.03 0.04 0.03 0.07 0.07 0.04 0.09 0.09 0.05 0.05 0.11 0.04 0.02 0.06 0.06 0.05 0.05 0.15 0.04 0.10 0.14 0.06 0.03 0.14 0.05 0.06 0.05 0.06 0.05 0.30 0.04 0.03 0.03 0.02 0.05 0.46 0.03 0.01 0.05 0.03 0.05 0.04 0.07 0.05 0.04 0.08 0.16 0.03 0.11 0.04 0.05 0.09 0.01 0.04 0.05 0.07 0.02 0.07 0.01

35 18 14 56 13 23 24 26 43 42 17 23 47 44 22 43 26 23 24 29 38 12 38 22 20 40 24 11 38 28 41 34 27 6 47 12 17 21 32 5 38 20 41 9 28 25 30 41 29 31 14 45 16 16 37 17 29 16 39 29 42 25 25


441

442

V. Monteiller and S. Chevrot Table A1. (Continued.) Station

Latitude

Longitude

φ

φ

δt(s)

δt(s)

Nbr.

YBH WBS WDC WEN WER WES WGR WLT WMC WSS WTT NE70 NE71 214 S08 U05 V05

41.73204 35.53664 40.57988 37.62211 35.06053 32.75904 34.51085 34.00948 33.57360 34.17170 33.94869 32.42095 31.68973 31.95590 37.49930 36.33560 35.86670

−122.71039 −118.14035 −122.54113 −121.75697 −119.02711 −115.73161 −119.27407 −117.95077 −116.67470 −118.64971 −118.25547 −115.26078 −115.90526 −112.81150 −118.17110 −120.12050 −119.90280

64 87 56 80 84 −78 −89 76 80 72 78 −69 82 25 58 84 71

1 1 1 2 2 3 1 3 1 1 1 10 2 20 1 9 11

0.96 1.10 0.62 1.10 1.31 0.66 0.96 0.98 0.80 1.22 1.15 0.88 0.89 0.19 1.41 0.61 1.13

0.04 0.01 0.03 0.09 0.38 0.09 0.05 0.10 0.04 0.08 0.04 0.30 0.09 0.14 0.04 0.18 0.48

35 33 46 22 5 23 36 22 35 32 20 6 13 9 9 8 4

APPENDIX B: FIT OF SPLITTING I N T E N S I T I E S AT I N D I V I D UA L S TAT I O N S

Figure B1. Bin averaged splitting intensities (open circles) and their best sinusoidal fits (black lines) at each individual stations. The phase and amplitude of these optimal sinusoids give, for each station, respectively the apparent fast direction and splitting delay. The sinusoidal fits predicted by the 3-D anisotropy model are also displayed with blue lines.



Figure B1. (Continued.)


C


443

444








C


445

446


