Crust–mantle coupling at the northern edge of the Tibetan plateau: Evidence from focal mechanisms and observations of seismic anisotropy

Tectonophysics 584 (2013) 221–229

Contents lists available at SciVerse ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

Crust–mantle coupling at the northern edge of the Tibetan plateau: Evidence from focal mechanisms and observations of seismic anisotropy Vadim Levin a,⁎, Guo-chin Dino Huang a, Steven Roecker b a b

Department of Earth and Planetary Sciences, Rutgers University, Piscataway, NJ, USA Department of Earth and Environmental Sciences, Rensellaer Polytechnic Institute, Troy, NY, USA

a r t i c l e

i n f o

Article history: Received 5 September 2011 Received in revised form 22 April 2012 Accepted 15 May 2012 Available online 25 May 2012 Keywords: Continental lithosphere strength Earthquake source mechanisms Seismic anisotropy Mantle deformation

a b s t r a c t The Tarim basin is distinct from the Tibetan plateau by the apparently low degree of internal deformation. Relative motion between the lithospheres of Tibet and Tarim likely results in coherent mantle deformation, with a diagnostic signature of anisotropic seismic wave speed. Birefringence (splitting) in core-refracted shear phases (SKS, PKS) observed along the Tibet–Tarim border is indicative of seismic anisotropy along their path. The mantle deformation direction inferred from shear wave splitting under the assumption of a single source of anisotropy is approximately E–W. The size of the splitting delays inferred in our analysis suggests that most of the anisotropic signature is coming from the mantle. Furthermore, our data suggest vertical stratification of rock fabric. Focal mechanisms of regional earthquakes, both published and derived from new data, provide evidence for ~NNE–SSW compression within the Tarim basin north of the border with Tibet. This type of internal deformation in the lithosphere is consistent with its underthrusting (subduction?) beneath Tibet. Within the crust of the Tibetan plateau we find a more complex pattern suggestive of the ~E–W extension, and generally consistent with the sense of motion on the Altyn Tagh fault that extends along the Tibet–Tarim boundary. Therefore, at the Tibet–Tarim border we find close similarity between deformation directions within the crust of the plateau and within the upper mantle on both sides of the boundary. A plausible explanation of such similarity would be the coupling of the crust and the upper mantle, with no weak zone being present in the lower crust. One scenario for such coupling would involve an extension of a deformation zone associated with the Altyn Tagh fault into the uppermost mantle, making this fault zone very similar to a plate boundary. © 2012 Elsevier B.V. All rights reserved.

1. Introduction A key constraint on the distribution of strength in the lithosphere is the presence (or absence) of earthquakes in the lower crust and upper mantle. Arguments in favor of locating continental lithosphere's strength in both the upper crust and the uppermost mantle rely on observations of crustal deformation patterns (e.g., Royden et al., 1997), the long-term response of lithosphere to loads (glaciers, mountain belts) (e.g., Burov and Diament, 1996), and the rare but significant occurrence of subcrustal earthquakes in the continental lithosphere (e.g., Monsalve et al., 2006). This conceptual model envisages two “strong” layers (brittle upper crust, brittle uppermost mantle) separated by a “weak” lower crustal zone. Whether earthquakes do occur beneath the crust of the continental lithosphere is not universally accepted (e.g., Jackson, 2002; Maggi et al., 2000a, 2000b). The alternate model places the entire strength of the continental lithosphere in the brittle crust, and envisages a weak mantle lithosphere (e.g., Jackson, 2002; Maggi et al., 2000b). A key prediction of this model is the extreme paucity (possibly a total

⁎ Corresponding author. Fax: + 1 732 445 5312. E-mail address: [email protected] (V. Levin). 0040-1951/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2012.05.013

absence) of seismic activity in the upper mantle of the continents. Thus the presence or absence of earthquake activity below the crust– mantle transition in the continental lithosphere is a likely discriminant between the two end-member models. For earthquakes of M > 5.5 in the continental lithosphere the determination of focal depth may be based on teleseismic data (e.g., Chen and Molnar, 1983; Maggi et al., 2000a). These earthquakes are not very frequent, and the short interval of the available instrumental seismic records (from 1960s to the present) has a potential to bias the overall distribution. Smaller events (Mb 5.5) are at least an order of magnitude more frequent, but it is difficult to obtain reliable locations and focal mechanisms for them using only teleseismic observations. Linked with arguments regarding vertical distribution of strength in the lithosphere is the consideration of crust–mantle coupling in actively deforming regions. If the crust is strong throughout, it is likely to transfer the stress across the Moho, and we can expect coherent deformation of the crust and the uppermost mantle (e.g., England and McKenzie, 1982). In the weak lower crust scenario ductile flow is envisaged within it (e.g., Bendick et al., 2000; Molnar and Lyon-Caen, 1989; Royden et al., 2008), and consequently the deformation patterns of the crust and the upper mantle do not have to be linked. Interestingly, arguments in favor of both scenarios often depend on observations of directional dependence

222

V. Levin et al. / Tectonophysics 584 (2013) 221–229

(anisotropy) of seismic wave speed. Thus, Shapiro et al., 2004 interpret the discrepancy in speeds of vertically and horizontally polarized surface waves as evidence for “flattening” (and hence horizontal flow) of the middle crust of Tibet. On the other hand, Holt, 2000 and Flesch et al., 2001 use similarity between deformation directions predicted by crustal deformation simulations and the direction of fast shear wave speed in the upper mantle as evidence of coherent deformation throughout the lithosphere. In this paper we focus on the northwestern part of the Tibetan plateau and the border with the Tarim basin (Fig. 1). This area is characterized by significant seismicity, including a recent Mw = 7.2 earthquake. Gravity modeling suggests Tarim lithosphere underthrusting the Tibetan plateau (Lyon-Caen and Molnar, 1984), while teleseismic imaging by Wittlinger et al. (2004) suggests the presence of the Tarim lithosphere in the upper mantle beneath the Tibetan plateau, a configuration consistent with subduction. The sinistral Altyn–Tagh Fault (ATF) is a major tectonic structure in this region. Evidence for an abrupt change in crustal thickness across it was reported by Wittlinger et al. (2004), while the compilation of active-source studies by Zhang et al. (2011) shows a gradual change from Tarim to Tibet. Numerous estimates of the motion rate on the ATF converge on ~10 mm/y (see review in Gold et al., 2011). How this motion is accommodated at depth is significant for the overall arguments on the lithospheric strength. Further east Zhang et al. (2007) combined GPS and geomorphological observations to argue for distributed “strike-slip” deformation in the crust around (primarily — south of) the ATF, while Hilley et al. (2009) used more recent GPS data in conjunction with earthquake cycle models to infer relatively high viscosity in the lower crust and the upper mantle adjacent to the ATF. In a previously published study (Huang et al., 2011) we identified a number of moderate sized earthquakes beneath this region that have lower-crustal hypocenters. Their focal mechanisms show a mix of normal and strike-slip motion, and suggest an overall NNE–SSW compression in the deep part of the crust. Among nine earthquakes deeper than 40 km found in the region by us and others (Chen and Yang, 2004; Fan and Ni, 1989) only one, at the far eastern edge of the region, shows an oblique thrusting mechanism that would be consistent with subduction of the Tarim lithosphere.

75˚

78˚

81˚

87˚

84˚

42˚

42˚

39˚

39˚

36˚

36˚

33˚

33˚

75˚

78˚

81˚

84˚

87˚

Fig. 1. A map of the study area and seismic stations used by this study. White inversetriangles show the seismic stations that contributed data used for moment tensor inversion and shear wave splitting analysis. Blue dashed rectangle shows the study area, green circles denote all earthquakes with focal mechanisms, both published and determined in this study. Gray lines show faults: ATF — Altyn–Tagh, KKF — Karakorum. Seismic activity from 1997 to 2008 from the China Earthquake Administration catalog is shown by red crosses.

In this paper we present a compilation of published and newly determined focal mechanisms for northwestern Tibet and the adjacent Tarim basin that, taken together, characterize the stress within the crust and, possibly, the uppermost mantle. We also explore the deformation of the upper mantle using published and newly analyzed measurements of shear wave splitting in teleseismic core-refracted waves. A comparison of these independent measures of deformation in the crust and the mantle provides a means for assessing the degree of coupling between them. 2. Data We use seismic data from various sources including temporary field deployments (Kao et al., 2001; Roecker et al., 2001; Wittlinger et al., 2004) and permanent seismic networks (the Kyrgyzstan telemetry seismic network— KN, and the Global Seismic Network — GSN). Selection of local earthquakes for focal mechanism analysis is based on the catalog archived by the China Earthquake Administration (CEA). In addition to the network of Kao et al. (2001) in the Tibet– Tarim border region, the seismic network in the Tien-Shan (Roecker et al., 2001) and the GSN station NIL in Pakistan provides good azimuthal constraints for seismic events that are located in our study area but were not recorded by Kao et al. (2001) (Fig. 1). Sites of the network operated by Kao et al. (2001) are listed in the Supplementary Table S1, all others are documented in the IRIS Data Management System archive (www.iris.edu). Events we analyze took place from 1992 to 2008. To study teleseismic core-refracted shear waves (PKS and SKS phases) we use data recorded by the temporary network in the Tarim–Tibet border region (Kao et al., 2001). In the discussion we also incorporate previously published results from Levin et al. (2008) that used data from a Sino-French temporary deployment (Wittlinger et al. (2004)). 3. Methods We perform regional Centroid Moment Tensor Inversion (rCMTI) analysis that uses full waveforms of the vertical, radial and transverse components to invert for earthquake source parameters and determine the focal depth (Huang et al., 2011; Kao et al., 1998). The algorithm involves testing a range of trial CMTs, and comparing resulting synthetic seismograms to selected observations. Azimuthal coverage and a choice of the velocity model influence the solution. While the synthetic seismogram computation technique we use is limited to 1D distributions of velocity, our algorithm allows the use of path-specific velocity models. We can thus accommodate, to a degree, the lateral changes in seismic wave speed within the region of study. Each source–receiver pair in the computation can have a path-specific velocity model. Given what we know about regional variations in wave speed, however, we incorporate only a small number of velocity profiles that reflect major changes in geologic structure. Errors introduced into rCMTs by potentially mislocated epicenters can be significant (Huang, 2007). As a first step in determining new focal mechanisms we test available published epicenters for each earthquake and compare the inversion results. We use catalogs of CEA, PDE/USGS, the International Seismological Center (ISC), EHB (Engdahl et al., 1998), and the global CMT (Dziewonski et al., 1981; http://www.globalcmt.org). For those rCMTs with published epicenters resulting in a larger waveform misfit, we perform epicentral relocation with additional constraints of wave azimuth in grid search (Uhrhammer et al., 2001). Additional details on this procedure are described in Huang et al. (2011). We use a basic cross-correlation technique (Levin et al., 1999), and also a group inversion method of Menke and Levin (2003) to perform shear wave splitting analysis. The former method assumes that a shear wave traversing an anisotropic region has been split exactly once, and

V. Levin et al. / Tectonophysics 584 (2013) 221–229

223

a 2-13 41

Date of Event = 1997 / 5 / 30 Origin Time = 17 : 54 : 54.6 (UT) Epicenter = 37.28 N 78.31 E Centroid depth = 16 km Mw = 4.9 Nodal Plane 1 = 101 / 25 / 60 Nodal Plane 2 = 312 / 68 / 103

2

3

Avg waveform misfit

Event ID:

0.385 0.380 0.375 0.370 0.365 0.360 0.355 0.350 0.345 0.340 0.335

Depth (km) P axis T axis 10 15 20 25 30 Azi= 10˚ Dist= 0442km Avg_misfit= 0.108 Filter_band= 0.020-0.050 Hz Station 1: WUS Z (misfit= 0.030 weight= 1.00) R (misfit= 0.136 weight= 1.00) T (misfit= 0.158 weight= 1.00)

Azi= 45˚ Dist= 1076km Avg_misfit= 0.258 Filter_band= 0.010-0.030 Hz Station 2: WMQ Z (misfit= 0.223 weight= 1.00) R (misfit= 0.265 weight= 1.00) T (misfit= 0.285 weight= 1.00) Azi= 230˚ Dist= 0609km Avg_misfit= 0.255 Filter_band= 0.010-0.040 Hz Station 3: NIL Z (misfit= 0.121 weight= 1.00) R (misfit= 0.346 weight= 1.00) T (misfit= 0.297 weight= 1.00)

Azi= 352˚ Dist= 0667km Avg_misfit= 0.721 Filter_band= 0.015-0.040 Hz Station 4: TLG Z (misfit= 0.735 weight= 1.00) R (misfit= 0.701 weight= 1.00) T (misfit= 0.729 weight= 1.00)

Time scale (sec.)

Average misfit for 4 stations = 0.335

0

50

Event ID:

2-32

Date of Event = 1999 / 7 / 18 Origin Time = 13 : 11 : 0.5 (UT) Epicenter = 37.11 N 77.35 E Centroid depth = 33 km Mw = 4.2 Nodal Plane 1 = 304 / 65 / 166 Nodal Plane 2 = 40 / 78 / 25

3 1

2 4

P axis T axis

Avg waveform misfit

b

0.60 0.58 0.56 0.54 0.52 0.50 0.48

20 25 30 35 40 Depth (km) Azi= 93˚ Dist= 0137km Avg_misfit= 0.344 Filter_band= 0.050-0.100 Hz Station 1: DUWA Z (misfit= 0.297 weight= 1.00) R (misfit= 0.213 weight= 1.00) T (misfit= 0.521 weight= 1.00)

Azi= 90˚ Dist= 0229km Avg_misfit= 0.564 Filter_band= 0.040-0.090 Hz Station 2: HTN1 Z (misfit= 0.446 weight= 1.00) R (misfit= 0.714 weight= 1.00) T (misfit= 0.533 weight= 1.00)

Azi= 328˚ Dist= 0343km Avg_misfit= 0.560 Filter_band= 0.035-0.080 Hz Station 3: WQIA Z (misfit= 0.455 weight= 1.00) R (misfit= 0.644 weight= 1.00) T (misfit= 0.582 weight= 1.00)

Azi= 105˚ Dist= 0318km Avg_misfit= 0.488 Filter_band= 0.040-0.090 Hz Station 4: WRKS Z (misfit= 0.267 weight= 1.00) R (misfit= 0.843 weight= 1.00) T (misfit= 0.354 weight= 1.00)

Average misfit for 4 stations = 0.489

Time scale (sec.) 0

50

Fig. 2. Examples of the rCMT analysis for two events, one using only data from remote stations (a) and another only data from the network of Kao et al. (2001) (b). Data are shown in black, synthetic seismograms in red. A misfit-depth dependence curves display clear global minima. We evaluate the uncertainty of the hypocentral depth determination using a misfit value of 5% over the global minimum (values shown in Table S3).

determines two parameters (a fast shear wave polarization and a delay between fast and slow components) that make horizontal components most similar. Choices made in identifying analysis windows and

band-pass filters can affect the outcome. We use this technique primarily as an assessment of consistency of measurements between sites, and also for different observation directions at individual

224

V. Levin et al. / Tectonophysics 584 (2013) 221–229

sites (or groups of nearby sites). We expect to see directional changes in splitting measurements if anisotropy of seismic wave speed varies with depth (Levin et al., 1999), and we expect site-tosite changes if it varies laterally. To determine parameters of anisotropy, we employ a more sophisticated technique that seeks a model that satisfies a group of observations simultaneously by computing synthetic seismograms in simple layered structures, and matching them with observed waveforms. As shown in Menke and Levin (2003), the values of parameters thus determined (orientations of the fast shear wave polarization and delays accumulated in vertically propagating S waves) represent true properties of the subsurface while determinations of fast polarizations and delays obtained from individual records are, necessarily, effective measures that depend on the mutual orientation of the anisotropic tensor and the ray path. Successful application of this method requires observations of (potentially) split shear waves from different backazimuths. A detailed discussion of this methodology and its application may be found in Levin et al. (2008).

75

76

77

78

79

4. Earthquake focal mechanisms 4.1. Focal depth distribution The seismic stations we used are located along the margins of the Tarim basin where seismic rays propagate through two main geological regimes, the Tarim basin where a thick sedimentary layer is present (for stations in the Tien-Shan), and Tibet where a thick crust is observed (for stations on the Tibetan plateau). Consequently, we use two 1D wave speed structures to account for these ray paths in Green's function computation (Supplementary Table S2). The selection of wave speed structure depends on the ray path to the seismic stations. We were able to determine new focal mechanisms for 26 events. Fig. 2 shows two examples of our inversions. One is based on data from permanent stations outside the region, and the other is based entirely on data from the temporary deployment in the southern Tarim basin (Kao et al., 2001). Similar images for all newly analyzed events are presented in the supplementary materials. Fig. 2 contains curves of total misfit

80

81

82

83

40 50 km 39

4 1

38

2 3

37

36

35

34

8000 4000 0

Topo(m)

-30

Depth (km)

Altyn-Tagh fault

-60 Line 1 Line 2 Line 3 Line 4

-90 -120

-150

-100

-50 0 50 Lateral distance (km)

100

150

Fig. 3. Map view and stacked vertical projections for all events with focal mechanisms. Red lines on the map represent two major strike-slip faults in northwestern Tibet (ATF, KKF, see Fig. 1), gray lines locate the profiles of hypocenter projection. The center of each projection plane is aligned with the ATF. All projections are in turn stacked onto one profile to exhibit the relative position of earthquakes with respect to the ATF. Hypocenters within each projection plane are color-coded. Dashed lines at 40 and 80 km delineate deepest earthquakes beneath the Tarim and the likely extent of the crust beneath Tibet (see text for details).

V. Levin et al. / Tectonophysics 584 (2013) 221–229

value as a function of hypocentral depth. If a number of epicentral locations were available for the same event (e.g., from the CMT catalog and the CEA catalog) we compared depth vs. misfit curves for them, and chose the epicenter yielding the smallest misfit. In many cases we found that relocating the epicenter leads to significant improvement of waveform fit. We used finite-difference location algorithm (Roecker et al., 2006) and 1-D models shown in Table S2 to determine new epicenter position. Table S3 notes which epicenter (published or newly determined) was chosen as the preferred result. We evaluate the depth uncertainty of the hypocenter using the width of the minimum in the misfit vs. depth curve. We draw a line at a misfit value 5% greater than the global minimum, and use resulting range of values as an estimate of hypocentral depth uncertainties. These values are given in Table S3. We assembled a total of 50 focal mechanisms and focal depths for our study region, including 26 earthquakes determined in this study, and 24 seismic events from the global CMT catalog (www.globalcmt. org) and others (e.g., Chen and Yang, 2004). Focal depths of these 50 earthquakes range from near-surface to about 100 km, while their moment magnitudes range from 2.6 to 7.1. To display the vertical distribution of these focal depths and their correspondence with the topography we project them onto four cross sections normal to the trace of the ATF, and plot resulting distance– depth distributions of hypocenters relative to the position of the ATF trace (Fig. 3). This projection insures that every hypocenter in the cumulative projection is placed correctly with respect to the curving fault trace. We can see that north of the ATF and within the Tarim basin, seismicity is restricted to the upper 40 km within the crust.

225

South of the ATF earthquakes are present from the surface to about 60 km, and there is a cluster of events at ~90 km depth.

4.2. Focal mechanisms Focal mechanisms form two unequally sized groups with distinct strain patterns, geographically separated by the Tibet–Tarim boundary. North of it, four available focal mechanisms (epicenters denoted by red squares) show thrust or oblique-thrust faulting. Their P axes are nearly horizontal, and cluster around the N30°E direction, while the T axes are near vertical. These mechanisms reveal ongoing internal deformation within the Tarim lithosphere, and reflect a dominant lateral compression regime in the crust just north of the boundary between Tibet and Tarim (Kao et al., 2001). South of the boundary, three distinct areas of seismicity are seen in the Tibetan part of our study area. One of them (the eastern-most purple oval in Fig. 4) is associated with the largest earthquake (Mw=7.1 on March 30, 2008, according to the global CMT) instrumentally recorded in this region. Clearly, normal faulting mechanisms dominate in this seismic zone, with only one rare case (event #36) of thrust faulting. Most P axes are nearly vertical, and T axes are close to E–W direction. Farther to the west, earthquakes form two broader zones where most focal mechanisms reflect strike-slip faulting. The P axes are more horizontal, and trend NNE– SSW; most T axes are sub-horizontal as well, with most aligned in the SSE–NNW direction. Overall, P-axes of the earthquakes in the northwestern part of the Tibetan plateau agree with the crustal deformation direction suggested by GPS observations (Wang et al., 2001). Extensional

Fig. 4. Focal mechanisms of earthquakes in northwestern Tibet. The beach-ball diagrams in blue stand for the solutions from this study, and the black ones are from the global CMT (www.globalcmt.org) and others (e.g., Chen and Yang, 2004) (see Table S3 for details). The text next to each diagram shows event number in Table S3 and its focal depth. Focal sphere size is proportional to the moment magnitude.

226

V. Levin et al. / Tectonophysics 584 (2013) 221–229

axes scatter over nearly 90° of range, but if the earthquake size is taken into account the dominant direction in this area is close to E–W (Fig. 5). To summarize, in the Tibet–Tarim border region earthquakes take place from near-surface to the bottom of the crust of the Tibetan plateau (~60 km), but reach only ~30 km depths beneath the Tarim basin. A small number of events likely take place in the upper mantle. Earthquakes ranging from Mw = 2.6 to Mw = 7.1 have either normal-fault or strike slip mechanisms, broadly consistent with the east–west extension of the crust. A notable exception is found in a group of events within the Tarim basin that show thrust-fault mechanisms consistent with dominant ~north–south compression within the crust.

5. Shear wave splitting We used data from the year-long deployment of Kao et al. (2001) to recover new observations for shear wave splitting in SKS and PKS phases. In addition to SKS phases from the South–Western Pacific that are typically analyzed by similar studies in Tibet, we identified a number of clear records of PKS phases from events in South and Central American subduction zones. Our data set thus contains two clusters of sources, with backazimuths in ranges of 95°–110° (SKS) and 325°– 345° (PKS). The difference in incoming directions of ~130° is close to optimal for resolving directional variation in splitting. It is approximately half-way between opposite and orthogonal directions, both of which may yield similar observations even if anisotropic properties vary with depth (see Levin et al., 1999). Individual measurements performed using a cross-correlation algorithm (Levin et al., 1999) are shown in Fig. 6 together with measurements previously obtained in the area by Levin et al. (2008). That earlier study used the same analytical strategy; thus the results are compatible. Values and errors of 28 new splitting observations, as well as parameters of earthquake sources used to obtain them are listed in Table S4. Most measurements obtained in the region yield delay values Δt on the order of 1 s

a

N

Scale

N

PTMw

or less. Relatively low values may reflect a technique-specific bias as cross-correlation method tends to return smaller delays in cases of noisy data (Menke and Levin, 2003). Combining two data sets on a map we find overall consistency, with many measurements yielding NNW– SSE fast polarization directions, sub-parallel to the strike of the Altyn– Tagh fault. However, the scatter in observations points to likely complexity of the anisotropic structure. Applying the Menke and Levin (2003) group inversion approach to develop constraints on anisotropic properties of seismic velocity, we use the strategy discussed in detail in Levin et al. (2008), with a few modifications. First, we explore only two classes of anisotropic models, those with 1 layer of anisotropic material and those with 2 layers. In both cases the symmetry axis of anisotropic velocity is horizontal, and may take any orientation. Second, capitalizing on the close spacing of sites in the Kao et al. (2001) network we group together observations from nearby sites and solve for common anisotropic models that fit them best. The groups overlap in their sampling of the subsurface, forming different realizations of data sampling the same volume, and each group contains more data, making inversions less dependent on individual records. Fig. 7 presents data distribution, maps of individual measurements, and inversion results for three groupings of new data, and also results from Levin et al. (2008). The similarity of inversion results for three subsets of the new data suggests that our assumption of a common structure for the entire network is reasonable. We obtain fast directions of 118°, 98° and 93°, and cumulative delays of 0.7, 0.6 and 0.7 s for the eastern, middle and western subsets, respectively. These results are quite similar to a finding of a fast direction of 108° and a cumulative delay of 1 s for a site further west (Levin et al., 2008). The good agreement between results obtained with independent data sets is particularly notable for the fact that the area sampled by Levin et al. (2008) is within the Tibetan plateau, while new results presented here sample on the other side of the Altyn–Tagh fault, mostly beneath the Tarim basin. For all subsets of our data we see a preference for 2-layer anisotropic models, as represented by the existence of 2-layer models that yield data fit measures smaller than 1-layer models. In the inversion for the 2-layer

7 6

W

E

5 4 S

3

S

b

38˚

36˚ 1 second

76˚ Fig. 5. (a) Strain axes projected on the focal sphere (separately, for P and T axes, respectively) and (b) classification of the focal mechanisms using a ternary diagram of Frohlich and Apperson (1992) for all events used in this study. A clear grouping of P and T axes is seen in focal sphere projections, with P axes forming a band consistent with NE–SW shortening, and T axes mostly indicative of near-horizontal extension, in general E–W direction with fairly wide variation. Note almost total absence of axes showing vertical shortening. The ternary diagram of focal mechanism types shows near absence of thrusts, a mechanism expected to dominate in compressive regime. Red symbols in all diagrams show four thrust events within the Tarim basin that display a consistently different stress pattern.

78˚

80˚

82˚

Fig. 6. The map of shear wave splitting measurements. Circular inset in the upper right shows teleseismic sources used in shear wave splitting analysis, where inverted triangle in the center shows our network. Seismic rays are coming from northwest and southeast directions. The splitting measurements (red — this study, blue — Levin et al. (2008)) are presented at 250 km deep piercing points of their respective rays. Observations with delays less than 0.05 s are shown by crossed black lines, the longer line being aligned with the backazimuth of a respective observation. These observations are treated as “nulls”, i.e. measurements that show no evidence of anisotropy. All other observations are shown as bars scaled with delay Δt and aligned with fast direction φ.

V. Levin et al. / Tectonophysics 584 (2013) 221–229

227

Fast directions in the lower layers of all chosen solutions fall in the range 114°–138°, with an average of 126.6°SE. Fast directions in the upper layer fall in the range from 50° to 100°, with an average of 73.8°NE. Delay values are approximately equal in both layers, and range from 0.24 to 0.72 s, so that the cumulative delay is between 0.8 and 1.3 s in all chosen solutions. This is in a good agreement with observed splitting delays in individual observations. Our 2-layer solution for the Tibet–Tarim border region differs from the 2-layer solution obtained by Levin et al. (2008) for the site on the Tibetan plateau. There, chosen models had a lower layer with a fast direction in a range 36°–86°, and an upper layer with a fast direction of 112°–120°. The lower layer made a significantly smaller contribution to the delay than the upper layer. In summary, our results of shear wave splitting analysis at the Tibet–Tarim border suggest the presence of anisotropy throughout the area. To first order (i.e., assuming there is only 1 layer of fabric at depth), attributes of anisotropy at all sites along the Altyn Tagh fault appear to be similar: we find ~1 s of delay in split waves, the orientation of the fast polarization is approximately east–west. Vertical changes are favored by data, and vertically-varying anisotropic parameters on the opposite (Tibet and Tarim) sides of the fault do differ in detail.

model we examine a total of 8100 combinations of fast directions, from 0 to 180 with a step of 2° for both the top and the bottom layers of the model. For each combination a range of delays from 0 to 3 s is searched. Depending on the data subset, between 730 and 910 combinations (or 8.5–11%) yield data fit values smaller than the 1-layer solution. These combinations of fast directions are shown in Fig. 7d by shading, and the combination of two fast axes yielding the absolute smallest data fit value is marked. We note that the 2-layer models with best “formal” measure of fit tend to have near-orthogonal fast directions, a condition that leads to a complete trade-off between layer thickness and the predicted splitting delay (Levin et al., 2008; Menke and Levin, 2003). As a result, the cumulative delay in such 2-layer models tends to be very large, a condition unlikely to be representative of the real Earth. To mitigate this complication, from the range of 2-layer models we select a subset of “realistic” models as follows: first we select 50 solutions with smallest values of data fit, and subsequently select from these 10 solutions that have the smallest cumulative delay. The logic of this choice, based on Occam's Razor, reflects the fact that in our parameterization the delay quantifies anisotropy strength in a layer, thus solutions with smallest delays correspond to the least amount of anisotropy required. Two-layer anisotropy solutions thus chosen are shown by open circles within the color field in lower row of diagrams in Fig. 7. They clearly fall within the field of solutions that are “better”, in the formal sense, then those with 1 layer of anisotropy, and they do not suffer from the “crossed layers” tradeoff. For all three subsets of the new data the outcomes of two-layer inversions are similar — shaded areas in plots of Fig. 7d occupy similar ranges of values, and “realistic” solutions within these areas are close.

6. Discussion Taken together, focal mechanisms from previous studies and those determined in this paper provide a consistent description of the deformation within the crust, and plausibly the uppermost mantle. In

a 76.5˚E 77˚E 77.5˚E 37˚N

b

79˚

81˚

79˚

81˚

79˚

81˚

116

38˚

1s

36.5˚N 37˚ 36˚N

from Levin et al. 2008

36˚

fast direction, deg fast direction, deg fast direction, deg fast direction, deg 0

45

90 135 1800

45

90 135 180 0

45

90 135 180 0

45

90 135 180 3

2

2

1

1

0

0

180 135

upper layer

d

3

delay, s

c

90 45 0

0

45

90 135 1800

45

90 135 1800

45

90 135 180 0

45

90 135 180

lower layer fast direction, deg from N Fig. 7. Results of the group inversion of shear wave splitting for 3 subsets of the network in the Tibet–Tarim border region. A result for a site farther south (from Levin et al., 2008) is shown for comparison.(a) Source distributions for each set.(b) Maps of individual measurements of shear wave splitting illustrating lateral coverage for each subset of data. Symbols as in Fig. 6. Thick purple bars show single-layer group solution values from (c).(c) Results of inversions for a single layer of anisotropy. Data fit measure for each data set graphed as a function of choices of two model parameters (fast polarization φ and delay time Δt). Orange shading shows the smaller values, and the global minima are marked by stars.(d) Results of inversions for two layers of anisotropy. Shading corresponds to all combinations of fast polarizations in two layers that yield data fit measures smaller than the best one-layer solution (darkest shading is for smallest fit measures). Crossing lines in the diagram show fast direction of the corresponding one-layer solution. A global minimum of the data fit measure is marked by a star on each plot. For each data set open dots show ten solutions with the smallest cumulative delay (a sum of Δt in both layers) drawn from 50 solutions with the smallest misfits. See text for details.

228

V. Levin et al. / Tectonophysics 584 (2013) 221–229

northwestern Tibet individual focal mechanisms (Figs. 4, 5) populate the entire range between purely normal and purely strike-slip solutions. Significantly, we find consistency of the inferred stress regime throughout the crust beneath the northwestern edge of the Tibetan plateau. Even earthquakes likely located in the uppermost mantle (events with ~ 90 km depth) are part of the same pattern. All these earthquakes are consistent with horizontal shearing of the crust, in general agreement with both the geological and the short-term motion inferred for the ATF (Cowgill et al. 2009; Gold et al. 2011; Hilley et al., 2009). In contrast, four thrust mechanisms within the Tarim basin reflect a different regime (Figs. 4, 5). Only these events are consistent with the proposed subduction of the Tarim lithosphere beneath that of the Tibet plateau (Wittlinger et al., 2004). Comparing SKS splitting and receiver function results, Levin et al. (2008) showed that most of the splitting signal in this area has to originate in the mantle. We believe this conclusion applies to new data as well. Beneath a continent we expect fast polarizations of split shear waves to represent deformation-induced systematic orientation of olivine crystals in the upper mantle (Nicolas and Christensen, 1987; Ribe, 1992), and can use them as proxies for the orientation of mantle deformation. Making an assumption of a single anisotropic layer (Fig. 7b,c), we find very close agreement between the orientation of the inferred mantle deformation and the strike of the ATF. Allowing vertical changes in anisotropic properties (by introducing 2 layers) we obtain a better fit to data. In these models one layer has the fast shear wave polarization aligned close to the strike of the fault. This layer is constrained better than the other, and also contributes most of the delay. The primary difference in models is that beneath Tibet this layer is the “top” one, while beneath the Tarim it is the “bottom”. Our inversions for layered anisotropy offer no constraints for the specific location of anisotropic layers at depth, only for their order from top to bottom. Thus we can speculate that the layer of fabric aligned to the Altyn–Tagh fault may be “correlated” across it, and that the fabric within it reflects the long-term motion of the fault. Then the differences in layered models of anisotropy would reflect the change in the nature of the lithosphere across the fault. For example, anisotropy within the upper layer beneath the Tarim basin may reflect pre-existing fabric of the Asian lithosphere. The lithosphere of the Tibetan plateau, on the other hand, appears to have E–W fabric aligned with the ATF, as shown by Levin et al. (2008) (left-most column in Fig. 7). This orientation of the fast shear wave speed is found throughout much of the plateau as well (e.g. Chen et al., 2010). In Fig. 8 we summarize indicators of deformation in the crust and the upper mantle beneath the Tibet–Tarim border. For the crust, we emphasize the clustering of T axes in a general WNW–ESE direction. For the mantle, we show orientations of fast polarizations in single and double-layer models for the new data and for the site previously studied by Levin et al. (2008). On the basis of our findings the ATF in our study area appears to be a lateral boundary separating distinct lithospheric blocks. To the south of it crustal earthquakes extend to depths ~60 km, and some earthquakes take place in the upper mantle at ~90 km, while to the north seismicity is restricted to the upper 40 km of the lithosphere (Fig. 3). Best-fitting models of seismic anisotropy distribution also differ for regions north and south of the fault. This view of the ATF aligns well with findings by Hilley et al. (2009), but differs somewhat from that of Zhang et al. (2007), who argue for continuous deformation across the ATF further to the east. We note that ample seismicity consistent with lateral shearing of the Tibetan plateau crust is in good agreement with Zhang et al. (2007) notion that the distributed deformation accommodating India–Eurasia collision is restricted to the plateau. Summarizing our results, at the Tibet–Tarim border we find close similarity between deformation directions within the crust of the plateau and within the upper mantle on both sides of the boundary. We also observe

a 38˚

37˚ N

36˚ E

35˚ 34˚ 75˚

76˚

77˚

78˚

79˚

80˚

81˚

82˚

83˚

84˚

b lower

38˚

upper

37˚ N

36˚ E

35˚ 34˚ 75˚

76˚

77˚

78˚

79˚

80˚

81˚

82˚

83˚

84˚

Fig. 8. Lithospheric deformation directions at the border of Tibet and Tarim. Crustal extension inferred from focal mechanisms is illustrated by the stereo-plot of T axes (from Fig. 5). Shear wave splitting fast directions (a proxy for upper mantle deformation) are shown by colored bars. (a) Purple bar shows fast directions obtained if one layer of anisotropy is assumed (from Fig. 7c). (b) Two-layer models for the Tarim basin are illustrated by black and magenta bars for lower and upper layers, respectively. Crossing bars show ranges of directions within each layer. See Fig. 7d and text for details. Results for a site close to Altyn Tagh fault are from Levin et al. (2008).

seismic activity throughout the crust of the plateau, attesting to the brittle regime within it. A plausible explanation of such similarity would be the coupling of the crust and the upper mantle, with no weak zone being present in the lower crust. One scenario for such coupling would involve an extension of a deformation zone associated with the ATF into the uppermost mantle, making this fault zone very similar to a plate boundary. Acknowledgments Authors are grateful to Dr. Kao for donating data from the temporary network, and to Dr. Engdahl for the use of the earthquake catalog he maintains. Online data bases maintained by IRIS DMC (waveforms), China Earthquake Administration, Global CMT project, and NEIC (earthquake hypocenters and mechanisms) were invaluable in this research. Funding for the project came from NSF grant EAR 0440062. Appendix A. Supplementary data Supplementary data to this article can be found online at http:// dx.doi.org/10.1016/j.tecto.2012.05.013. References Bendick, R., Bilham, R., Freymueller, J.T., Larson, K., Yin, G., 2000. Geodetic evidence for a low slip rate in the Altyn Tagh fault system. Nature 404, 69–72, http://dx.doi.org/ 10.1038/35003555. Burov, E., Diament, M., 1996. Isostasy, equivalent elastic thickness, and inelastic rheology of continents and oceans. Geology 24, 419–422.

V. Levin et al. / Tectonophysics 584 (2013) 221–229 Chen, W.-P., Molnar, P., 1983. Focal depths of intra-continental and intraplate earthquakes and their implications for the thermal and mechanical properties of the lithosphere. Journal of Geophysical Research 88, 4183–4214. Chen, W.-P., Yang, Z.-H., 2004. Earthquakes beneath the Himalayas and Tibet: evidence for strong lithospheric mantle. Science 304, 1949–1952. Chen, W.P., Martin, M., Tseng, T.L., Nowack, R.L., Hung, S.H., Huang, B.S., 2010. Shear-wave birefringence and current configuration of converging lithosphere under Tibet. Earth and Planetary Science Letters 295, 297–304. Cowgill, E., Gold, R.D., Chen, X., Wang, X.-F., Arrowsmith, J.R., Southon, J.R., 2009. Low Quaternary slip rate reconciles geodetic and geologic rates along the Altyn Tagh fault, northwestern Tibet. Geology 37 (7), 647–650, http://dx.doi.org/10.1130/ G25623A.1. Dziewonski, A.M., Chou, T.-A., Woodhouse, J.H., 1981. Determination of earthquake source parameters from waveform data for studies of global and regional seismicity. Journal of Geophysical Research 86, 2825–2852. Engdahl, E.R., van der Hilst, R., Buland, R., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bulletin of the Seismological Society of America 88, 722–743. England, P., McKenzie, D., 1982. A thin and viscous sheet model for continental deformation. Geophysical Journal of the Royal Astronomical Society 70 (2), 295–321. Fan, G., Ni, J.F., 1989. Source parameters of the 13 February 1980 Karakorum earthquake. Bulletin of the Seismological Society of America 79, 945–954. Flesch, L.M., Haines, A.J., Holt, W.E., 2001. Dynamics of the India–Eurasia collision zone. Journal of Geophysical Research 106, 16435–16460. Frohlich, C., Apperson, K.D., 1992. Earthquake focal mechanisms, moment tensors, and the consistency of seismic activity near plate boundaries. Tectonics 11, 279–296. Gold, R.D., Cowgill, E., Arrowsmith, J.R., Chen, X., Sharp, W.D., Cooper, K.M., Wang, X.-F., 2011. Faulted terrace risers place new constraints on the late Quaternary slip rate for the central Altyn Tagh fault, northwest Tibet. GSA Bulletin 123 (5/6), 958–978, http://dx.doi.org/10.1130/B30207.1. Hilley, G.E., Johnson, K.M., Wang, M., Shen, Z.-K., Bürgmann, R., 2009. Earthquake-cycle deformation and fault slip rates in northern Tibet. Geology 37 (1), 31–34, http:// dx.doi.org/10.1130/G25157A.1. Holt, W.E., 2000. Correlated crust and mantle strain fields in Tibet. Geology 28, 67–70. Huang, G.-C., 2007. 3-D lithospheric structure and seismotectonics of the central Himalayan region, Ph.D. thesis, 233 pp., State Univ. of N. Y. at Binghamton, Binghamton. Huang, G.-C.D., Roecker, S.W., Levin, V., 2011. Lower-crustal earthquakes in the West Kunlun range. Geophysical Research Letters 38, L01314, http://dx.doi.org/ 10.1029/2010GL045893. Jackson, J., 2002. Strength of the continental lithosphere: time to abandon the jelly sandwich? GSA Today 12, 4–10. Kao, H., Jian, P.-R., Ma, K.-F., Huang, B.-S., Liu, C.-C., 1998. Moment tensor inversion for offshore earthquakes east of Taiwan and their implications to regional collision. Geophysical Research Letters 25, 3619–3622. Kao, H., Gao, R., Rau, R.-J., Shi, D.-N., Chen, R.-Y., Guan, Y., Wu, F., 2001. Seismic image of the Tarim basin and its collision with Tibet. Geology 29, 575–578. Levin, V., Menke, W., Park, J., 1999. Shear-wave splitting in Appalachians and Urals: a case for multilayered anisotropy. Journal of Geophysical Research 104, 17975–17994. Levin, V., Roecker, S.W., Graham, P., Hosseini, A., 2008. Seismic anisotropy indicators in Western Tibet: shear wave splitting and receiver function analysis. Tectonophysics 462, 99–108. Lyon-Caen, H., Molnar, P., 1984. Gravity anomalies and the structure of western Tibet and the southern Tarim basin. Geophysical Research Letters 11, 1251–1254.

229

Maggi, A., Jackson, J.A., Mckenzie, D., Priestley, K., 2000a. Earthquakes focal depth and effective elastic thickness, and the strength of the continental lithosphere. Geology 28, 495–498. Maggi, A., Jackson, J.A., Priestley, K., Baker, C., 2000b. A re-assessment of focal depth distributions in southern Iran, the Tien Shan and northern India: do earthquakes really occur in the continental mantle? Geophysical Journal International 143, 629–661. Menke, W., Levin, V., 2003. The cross-convolution method for interpreting SKS splitting observations, with application to one and two-layer anisotropic earth models. Geophysical Journal International 154, 379–392. Molnar, P., Lyon-Caen, H., 1989. Fault plane solutions of earthquakes and active tectonics of the Tibetan plateau and its margins. Geophysical Journal International 99, 123–153. Monsalve, G., Sheehan, A., Schulte-Pelkum, V., Rajaure, S., Pandey, M., Wu, F., 2006. Seismicity and 1-D velocity structure of the Himalayan collision zone: earthquakes in the crust and upper mantle. Journal of Geophysical Research 111, B10301, http://dx.doi.org/10.1029/2005JB004062. Nicolas, A., Christensen, N.I., 1987. Formation of anisotropy in upper mantle peridotites — a review. Composition, Structure and Dynamics of the Lithosphere–Asthenosphere System: In: Fuchs, K., Froidevaux, C. (Eds.), Geodynamics Series, 16. American Geophysical Union, Washington DC, pp. 111–123. Ribe, N.M., 1992. On the relationship between seismic anisotropy and finite strain. Journal of Geophysical Research 97, 8737–8748. Roecker, S.W., Xu, Y., Martynov, V., 2001. Local and regional earthquake depths in the Tien Shan: evidence for a strong lower crust in an active orogen. EOS. Transactions of the American Geophysical Union 82, 47. Roecker, S., Thurber, C., Roberts, K., Powell, L., 2006. Refining the image of the San Andreas Fault near Parkfield, California using a finite difference travel time computation technique. Tectonophysics 426, 189–205. Royden, L.H., Burchfiel, B.C., King, R.W., Wang, E., Chen, Z., Shen, F., Liu, Y., 1997. Surface deformation and lower crustal flow in eastern Tibet. Science 276, 788–790. Royden, L.H., Burchfiel, B.C., van der Hilst, R.D., 2008. The geological evolution of the Tibetan plateau. Science 321, 1054–1058. Shapiro, N.M., Ritzwoller, M.H., Molnar, P., Levin, V., 2004. Thinning and flow of Tibetan crust constrained by seismic anisotropy. Science 305, 233–326. Uhrhammer, R.A., Dreger, D., Romanowicz, B., 2001. Best practice in earthquake location using broadband three-component seismic waveform data. Pure and Applied Geophysics 158, 259–276. Wang, Q., et al., 2001. Present-day crustal deformation in China constrained by Global Positioning System measurements. Science 294, 574–577. Wittlinger, G., Vergne, J., Tapponnier, P., Farra, V., Poupinet, G., Jiang, M., Su, H., Herquel, G., Paul, A., 2004. Teleseismic imaging of subducting lithosphere and Moho offsets beneath western Tibet. Earth and Planetary Science Letters 221, 117–130. Zhang, P.-Z., Molnar, P., Xu, X., 2007. Late Quaternary and present-day rates of slip along the Altyn Tagh fault, northern margin of the Tibetan Plateau. Tectonics 26, TC5010, http://dx.doi.org/10.1029/2006TC002014. Zhang, Z., Deng, Y., Teng, J., Wang, C., Gao, R., Chen, Y., Fan, W., 2011. An overview of the crustal structure of the Tibetan plateau after 35 years of deep seismic soundings. Journal of Asian Earth Sciences 40, 977–989, http://dx.doi.org/10.1016/j.jseaes. 2010.03.010.