Seismic Imaging of Microblocks and Weak Zones in the Crust Beneath the Southeastern Margin of the Tibetan Plateau

The studies of Earth's history and of the physical and chemical properties of the substances that make up our planet, are of great significance to our understanding both of its past and its future. The geological and other environmental processes on Earth and the composition of the planet are of vital importance in locating and harnessing its resources. This book is primarily written for research scholars, geologists, civil engineers, mining engineers, and environmentalists. Hopefully the text will be used by students, and it will continue to be of value to them throughout their subsequent professional and research careers. This does not mean to infer that the book was written solely or mainly with the student in mind. Indeed from the point of view of the researcher in Earth and Environmental Science it could be argued that this text contains more detail than he will require in his initial studies or research.


Introduction
The southeast margin of the Tibetan Plateau lies between the heartland of the plateau to the west and the stable south China block to the east, spanning from western Sichuan to central Yunnan in southwest China. Based on low-gradient topographic slope and lack of largescale young crustal shortening at the southeast plateau margin, Royden et al. (1997) and Clark and Royden (2000) proposed a channel-flow model in which a weak (low-viscosity) zone exists in the mid-to lower crust. Gravitational potential drives crustal materials from the Tibetan Plateau outward through the channel, creating a broad and topographically gentle margin and also accumulating stress near the strong crust of the Sichuan Basin. Using GPS data collected from the Crustal Motion Observation Network of China between 1998and 2004, Shen et al. (2005 showed that the crust is fragmented into tectonic blocks of various sizes, separated by strike-slip and transtensional faults ( Figure 1). They proposed a model for Tibetan Plateau deformation in which a mechanically weak lower crust experiences distributed deformation underlying a stronger, highly fragmented upper crust. On May 12, 2008, a destructive Ms 8.0 earthquake occurred along the Longmen Shan Fault, located between the eastern margin of the Tibetan Plateau and the Sichuan Basin (Burchfiel et al., 2008). It ruptured mainly toward the northeast over a length of ~270 km along the northeast-trending fault, with coseismic slip mainly consisting of thrust-and right lateral strike-slip components (Wang et al., 2008b). No noticeable precursors were observed before the main shock, which was anticipated because GPS modeling showed very low right-slip (~1 mm/yr) and convergence (<~3 mm/yr) rates along the Longmen Shan boundary (Meade, 2007). A deep process involving channel flow is hypothesized to be responsible for the 2008 Wenchuan Ms 8.0 earthquake (Burchfiel, et al., 2008;Teng et al., 2008;Zhang et al., 2008). Other models than the channel flow model such as the block model were also proposed for causing this earthquake (e.g. Hubbard and Shaw, 2009). www.intechopen.com Regional seismic tomography studies using body waves (Huang et al., 2002;Wang et al., 2003;Wang et al., 2007;Huang et al., 2009;Xu and Song, 2010) and surface waves Huang et al., 2010;Li et al., 2010) found widespread low velocity zones in the mid-and lower crust, supporting the channel-flow model proposed by Clark and Royden (2000). Receiver function analysis on stations in southwest China also identified low velocity zones (LVZs) in the mid-and lower crust and high average Poisson's ratio in the crust (e.g. Xu et al., 2007;Wang et al., 2008a;Liu et al., 2009;Zhang et al., 2009c). In addition, magnetotelluric (MT) sounding detected low resistivity layers in the middle and lower crust (e.g. Sun et al., 2003;Zhao et al., 2008;Bai et al., 2010). These low velocity and low resistivity zones were interpreted to be caused by partial melt. In this article, we present the results of a joint inversion for Vp, Vs, and Vp/Vs models, applying a modified double-difference seismic tomography method to the catalog picks collected by the Seismological Bureau of Sichuan Province for the period [2001][2002][2003][2004]. The joint interpretation of three models permits a more complete characterization of the mechanical properties and geological identity of crustal materials and therefore is helpful for better understanding the cause of the low velocity and low resistivity layers. Compared to the previous regional tomography studies in the Sichuan region, this is the first time that a Vp/Vs model is directly inverted from S and P arrival times instead of from dividing Vp by Vs. The three-dimensional (3D) shear-wave velocity model of Yao et al. (2008) indicated that the LVZs vary considerably in strength and depth range and faults may mark lateral boundaries of the LVZs. Our high-resolution 3D Vp, Vs, and Vp/Vs models are utilized to examine the spatial distribution of and interconnectivity between LVZs, which is important for understanding the tectonic block motions (Shen et al., 2005). For accurately calculating ray paths and travel times between events and stations in the case of strong velocity heterogeneity, a spherical-earth finite-difference (SEFD) travel time calculation method is developed and tested.

Spherical-Earth Finite-Difference (SEFD) travel time calculation
Since their introduction to seismology by Vidale (1988), finite difference solutions to the eikonal equation have enjoyed widespread application as a robust and efficient technique for computing travel times in heterogeneous media. To the extent that one can easily access the travel time tables produced by such techniques, they can be readily incorporated into earthquake location and tomographic imaging algorithms (e.g. Nelson and Vidale, 1990;Hole, 1992). With few exceptions (Fowler, 1994;Schneider, 1995) where s is the local slowness. To the extent that there is no significant spatial regularity in the heterogeneity that we are attempting to parameterize, the bias that we introduce by a particular choice of grid system, Cartesian or otherwise, will not be significant or in the worst case will increase the level of model noise.
As a simple consequence of gravity and temperature, wavespeeds in the earth are primarily a function of depth; lateral variations in wavespeed often tend to be only a few percent. Over regional distances on the order of ~200 km or less, such depth variations should for most purposes be modeled adequately by a Cartesian grid. However, there is a potential for introducing a model induced signal into an inversion when at greater distances the radial variations in wavespeed do not correlate well with the Cartesian grid. One strategy for coping with sphericity is to employ earth flattening (e.g. Abers and Roecker, 1991) but the transformations for flattening are not appropriate for a laterally heterogeneous medium, and moreover there are issues with computing distance properly in the flattened frame (in particular they should always be computed along great circles). Another strategy is to simply put a round earth in a rectangular box, known as the sphere-in-a-box method (Flanagan et al., 2007), but this can artificially introduce anisotropy into the model because radial gradients are not represented the same way in all directions. Of course, such artifacts can be reduced by decreasing the grid spacing but resulting increase in the number of grid points could make the computations intractable.
As an alternative, one might consider solving the eikonal equation in a spherical coordinate system, so that radial gradients are parameterized equally throughout the model with a reduced number of grid points. The eikonal equation in spherical coordinates is: where r is the radius from center of the earth, dr is positive away from the center, and |dr| = h;  is the co-latitude (0° at north pole, 90° at equator), d is positive to the south, and |d| = ;  is longitude, d is positive to the east, and |d| = ; and s is slowness. To solve this system, we must be account for the differences in r, , and  for each node in the mesh. For each node i we assign r i ,  i ,  i , and also signs for directional purposes (Table  1). We derive expressions for each of the finite difference (FD) "stencils" used in the algorithm. For example, when applying Scheme A of Vidale (1990), we compute the time at one point given the times at 7 adjacent points in the 8-point cell. Referring to Figure 2 and Table 1, the FD derivatives are:   Given the values for t 0 through t 6 , this expression can be rewritten in the form at 7 2 + bt 7 + c = 0 to solve for t 7 , with coefficients a, b and c defined as follows:   Comparable equations, which are included in the Appendix, can be derived for the "edge" and "face" stencils of Vidale (1990). One of the problems encountered with these finite difference techniques is that the travel times at the grid points in the immediate neighborhood of the starting point need to be assigned somehow. As long as the wavespeeds are not overly heterogeneous near the starting point, integration of slowness along a straight line path provides a reasonable estimate of travel time. This may not always be the case, however, and in any event as Vidale (1988) pointed out the finite difference approach does not work well when there is significant wavefront curvature over the size of the grid volume element. One efficacious way to solve both of these problems is to use a cascading approach by defining a fine grid in www.intechopen.com the vicinity of the starting point and a coarser grid outside that region. We have adopted this approach. We tested the SEFD method by calculating travel times in an analytical velocity model V=V 0 (r 0 /r), where V 0 =4.0km/s, r 0 is the Earth's radius and r is the distance between the source and the Earth's surface. Figure 3a shows the analytical travel times for a source located at latitude 21.2° and longitude 121.75°. We discretized the model into a 3D grid with a grid interval of 0.1° in latitude and longitude and 10 km in depth. The source region is set up to be 3 grid nodes in which the travel times are calculated along a straight-line path. The differences in travel times compared to analytic times are shown in Figure 3b. The travel time error around the source is as much as 1.08 s. Outside the source region, the mean travel time error is 0.108 s, and is everywhere generally smaller than 0.3 s. Along the latitude www.intechopen.com and longitude directions and the directions between them, the travel time errors are relatively small due to the design of the stencils. To deal with the inaccuracy problem near the source region, we applied a cascading-grid strategy, in which a fine grid is used near the source region and a coarse one is used outside the source region. The grid interval inside the source region is 10 times smaller than that outside. The resulting travel time error near the source is much smaller than before, down to 0.17 s and the mean travel time error decreases to 0.087 s. The tests show that the cascading-grid strategy improves the travel time accuracy near the source region and can also decrease the travel time error away from the source region. We also calculated the travel times using the "sphere-in-abox" method, in which the travel times are calculated on a 3D Cartesian grid with a uniform grid interval using the finite-difference eikonal solver of Podvin and Lecomte (1991). The velocity values on Cartesian grid nodes are linearly interpolated from 8 surrounding spherical grid nodes. The grid interval is set to be 5 km, about 2 times smaller than that used for the SEFD travel time calculation. The travel time errors from Cartesian grid FD method are plotted in Figure 3d. It can be seen that the travel time errors around the source region are small. This is because the FD scheme used in Podvin and Lecomte (1991) adopted an initialization procedure to accurately calculate the travel times around the source. Similar to our SEFD method, the travel time errors are small along latitude, longitude and their middle intersections. However, the travel time errors outside the source region are relatively large. The overall mean travel time error is 0.312 s, much greater than 0.108 s and 0.087 s for the two SEFD cases. This is mainly due to the inaccuracy in velocity values on Cartesian grid nodes when they are interpolated from the exact spherical grid nodes. Even when the Cartesian grid interval is finer, the travel time errors are still greater compared to the case using spherical grid.

Seismic tomography method
We employed a new version of the double-difference (DD) seismic tomography method that simultaneously solves for Vp, Vs, Vp/Vs and event locations using both absolute and differential P, S, and S-P times (Zhang, 2003;Zhang et al., 2009a, b). This new code, named tomoDDPS, avoids the pitfalls of inferring Vp/Vs from Vp and Vs models via division (Eberhart-Phillips, 1990). We briefly summarize the method as follows. The P and S arrival times p T and s T from an earthquake i to a seismic station k are expressed using ray theory as path integrals where i  is the origin time of event i , p u and s u are the P-and S-wave slowness fields and dl is an element of path length. The source coordinates 123 (, ,) , xxx origin times, ray paths, and the slowness field are the unknowns. By assuming the ray paths of P and S waves are identical, which is true when Vp/Vs is constant, Vp/Vs can be determined from S-P arrival times Ts Tp  , as follows (Thurber, 1993), www.intechopen.com Note here because P and S waves from the same event share the same origin time, the unknown origin times are removed from this equation. In the simul2000 algorithm (Thurber and Eberhart-Phillips, 1999), Equations (6) and (8) are used to solve for Vp and Vp/Vs using P and S-P times and Vs is later calculated by dividing Vp by Vp/Vs. However, as noted by Wagner et al. (2005), the Vs model may be biased if calculated in this way because the anomaly in Vp may leak into Vs. In the new tomoDDPS algorithm, Vp, Vs, and Vp/Vs are determined simultaneously in a system using P, S, and S-P times based on Equations (6), (7) and (8) (Zhang, 2003;Zhang et al., 2009a, b). To meet the assumptions made for Equation (8), only S-P times from similar P and S ray paths are selected to solve for Vp/Vs. Similar to the DD tomography code tomoDD, differential P and S times are also used in tomoDDPS to better constrain seismic event locations and Vp and Vs models (Zhang and Thurber, 2003). In addition, differential S-P times are also used to determine the Vp/Vs structure based on the differential time version of Equation (8), which can be directly constructed from differential P and S times. One advantage of using differential S-P times is to remove the effect of different ray paths of P and S waves outside the source region. Near the source region, P-and S-wave ray paths are generally close to each other. Smoothing weights are applied to P-and S-wave slowness perturbations and Vp/Vs perturbations for neighboring inversion grid nodes to stabilize the tomographic inverse problem. The complete tomographic system is represented as follows (Zhang et al., 2009a):  is the origin time perturbation, u  is the P or S slowness perturbation, (/) Vp Vs  is the V p /V s perturbation, w 1 and w 2 are data weights for the absolute and differential P or S data, w 3 and w 4 are data weights for the absolute and differential S-P data, w 5 and w 6 are smoothing weights for slowness and V p /V s models, and m and n indicate neighboring inversion grid nodes. The complete system is solved using a www.intechopen.com damped least squares inversion method LSQR in which the weighted data residuals are minimized (Paige and Saunders, 1987).

Data and inversion details
For the Sichuan region, we collected ~38,600 P-and ~36,500 S-wave first arrival times from 4878 earthquakes observed on 55 stations for the period 2001 to 2004 (Figure 1). These arrival times are selected from the original catalog data based on the major trend of travel time curves (Figure 4). There are obvious 60-second clock shift errors and other reading errors in catalog picks. For each event included in the analysis, there are at least 6 P and 2 S observations, increasing the likelihood of reliable relocations. From the absolute P and S arrival times, we constructed ~269,000 P and ~261,000 S differential times. The average number of differential times (links) per event pair is 11 and the average hypocentral separation (based on catalog locations) for the linked event pairs is ~11 km. The inversion grid interval for the velocity model in latitude and longitude is 0.5°. In depth, the grid nodes were positioned at 0, 5, 10, 17.5, 25, 35, 45, 65, and 90 km. In the Sichuan region, the Moho depth varies from ~60 km in the Songpan-Ganze terrane to ~46 km in the Sichuan basin (Xu et al., 2007). Therefore our model mainly reflects the crustal structure of the southeastern Tibetan Plateau. We first derived a minimum one-dimensional (1D) velocity model for the region based on the regional 1D velocity model of Zhao et al. (1997) ( Figure 5). The travel times were calculated using the new SEFD method described above. Both damping and first-order smoothing were used to stabilize the inversion. A trade-off analysis between data variance and model variance was used to select optimum damping and smoothing parameters. The initial unweighted root-mean-square (RMS) travel time residual of 1.78 s was reduced to a final value of 0.48 s, a reduction of approximately 73%. We assess the model quality by a checkerboard resolution test. ±5% velocity anomalies were added to the final 3D Vp and Vs models with an anomaly size of one grid node (Figures 6  and 7). The velocity anomalies for Vp and Vs are made opposite in sign so that the Vp/Vs anomaly ranges from approximately -9% to 11% (Figure 8). A combination of constant noise for each station and random noise at a level comparable to the final inversion misfit is added to the absolute P and S times. The checkerboard resolution test showed that both Vp and Vs models are relatively well resolved for the depth range of 5 to 65 km except for the depth slice of 17.5 km. For the Vp/Vs model, it is also well resolved from a depth of 5 to 45 km except for the depth slice of 17.5 km. For the depth slice of 65 km, the Vp/Vs model has some resolution in the middle part of the model. All three models have poor resolution at depth 0 km.        Figure 9 shows the horizontal slices of the Vp, Vs and Vp/Vs models at depths of 0, 5, 10, 17.5, 25, 35, 45 and 65 km. Figures 10, 11 and 12 show the cross sections of the Vp perturbations, Vs perturbations and Vp/Vs model at latitudes of 28°, 29°, 30°, 31° and 32°, respectively. Vp and Vs perturbations are with respect to the 1D Vp and Vs models averaged from the 3D models ( Figure 5). Compared to previous tomography studies in the Sichuan region, to our knowledge this is the first time that the Vp/Vs model is directly determined from the S-P arrival times instead of being derived from the Vp and Vs models. Pei et al. (2010) used the same methodology of Zhang et al. (2009a, b) to solve for the Vp, Vs and Vp/Vs models around the Longmenshan Fault, but the depth extent of the model is only down to 30 km. At shallow depths (0-5 km), our tomographic velocity models are consistent with the local geology. The Sichuan Basin is clearly imaged as a low velocity region with high Vp/Vs ratios (>1.9). The basin is mainly composed of Tertiary, Quaternary to Mesozoic sediments derived from uplift resulting from the collision. Outside the Sichuan Basin, velocities are generally higher and the Vp/Vs ratio is lower than 1.7, mainly corresponding to the Songpan-Ganze Terrane. One low velocity zone (low Vp and Vs and high Vp/Vs) in the Sichuan Basin is located between latitudes 31° to 32° and longitudes 104° to 105°, extending all the way from the surface down to the depth of a depth of 25 km. Although this low velocity zone is located around the model edge, the checkerboard resolution tests showed this area has good resolution. Previous tomography studies also identified this low velocity anomaly ). At deeper depths (10-25 km), strong velocity variations are present across the region. At 10 km depth, the Longmen Shan Fault (LMSF) separates a higher velocity region to the northwest from lower velocities to the southeast. At 17.5 km depth, the velocity contrast is still clear north of latitude ~30.5°, especially for Vs, indicating the LMSF may penetrate at least down to ~18 km. The Longquan anticline separates a relatively lower velocity region to the west and a higher velocity to the east at depths 10 and 17.5 km.

Results and discussion
In the Songpan-Ganze Terrane, there are scattered low velocity regions bounded by high velocity bodies in the depth slices of 10 km and 17.5 km. Especially at 10 km, the velocity pattern resembles the deformation block model found by modeling the GPS data by Shen et al. (2005). The low velocity anomalies generally follow the derived sub-block boundaries. For example, the Songpan-Xihe deformation zone separating the Ahba block to the northwest and the Longmen Shan block to the southeast corresponds to a broad low velocity zone.       Starting from the depth slice of 25 km and down to 65 km, we see widespread low velocity zones outside the Sichuan Basin, which itself generally corresponds to a high velocity anomaly. These low velocity zones are not uniformly distributed but vary in amplitude and they are mostly connected to each other. Previous surface wave and body wave tomography studies also found crustal low velocity zones in this region and interpreted them as weak zones for possible channel flow (e.g. Yao et al., 2008Yao et al., , 2010Wang et al., 2009) Figures 10 and 11), we see low velocity anomalies below ~20 km depth underneath the Songpan-Ganze block. In comparison, the region beneath the Sichuan basin is shown as a high velocity body. For example, in the cross section of latitude 30°, there is an evident low velocity layer in both Vp and Vs around 20 km depth from longitude 101.5° to 102.5°. This low velocity layer was previously detected from a deep seismic sounding profile and by receiver function analysis  along latitude 30°. A magnetotelluric (MT) survey between longitudes 102° to 104° and slightly to the north of latitude 29° also showed a low resistivity layer around 20 km depth (Zhao et al., 2008). This low resistivity layer is associated with low Vp, low Vs and high Vp/Vs (Figure 12). Receiver function analysis at station MC09 (latitude 29° and longitude 102.8°) also showed a low Vs layer in the crust associated with a high Vp/Vs ratio (Xu et al., 2007). Because of high regional surface heat flow values (Hu et al., 2000), the low Vp and Vs anomalies in the eastern Tibetan Plateau have been suggested by many researchers to be due to elevated temperatures or partial melt (e.g. Yao et al., 2008). However, the recent receiver function analysis by Robert et al. (2010) found low Vp/Vs (=1.69) beneath the Songpan-Ganze terrane, which is lower than the mean value for continental areas (Zandt and Ammon, 1995). This observation led them to dispute the existence of a thick and extensive zone of partial melt in the crust of the Songpan-Ganze terrane. In contrast, Wang et al. (2008) showed high Vp/Vs (or Poisson ratio) perturbations and Xu and Song (2010) showed high Poisson ratios in the middle and lower crust of the Songpan-Ganze terrane. Both of their models are obtained by directly dividing Vp by Vs. In comparison, our Vp/Vs model is obtained by directly inverting absolute and differential S-P times and thus is more reliable. Above ~35 km, most of the Songpan-Ganze terrane has a Vp/Vs value of ~1.6. Below ~35 km, the Vp/Vs value starts to increase and reaches up to ~1.85. By averaging from 0 to 60 km, the Vp/Vs value from our study is close to what Robert et al. (2008) found from the receiver function analysis. Although low velocity anomalies start from ~20 km depth in the crust of the eastern Tibetan Plateau, partial melt may not exist until 35 to 40 km depth where the Vp/Vs ratio is relatively high (Christensen, 1996). In the shallower part, the low velocity anomalies could be caused by the existence of aqueous fluids (Li et al., 2003).
The distribution of aqueous fluids in spheroidal pores can result in low Vp and Vs anomalies and an average to low Vp/Vs anomaly (Takei, 2002). In this study area, the earthquakes are mostly located above 30 km (Figure 13), supporting the argument that partial melt does not occur above this depth. Most of the earthquakes are located in the region with Vp>4.4 km/s and Vs>2.6 km. Laboratory experiments show that even a small percentage of melt dramatically reduces the viscosity of rock and it essentially loses its solid nature and behaves like fluid when the melt content reaches 20% to 55% (Kohlstedt and Zimmerman, 1996). A rock containing a fluid content greater than 5% is 10 times weaker than the surrounding material with the same composition (Rosenberg and Handy, 2005). An MT survey through latitude 30° showed a strong low resistivity anomaly below ~30 km depth, which may contain 5% to 20% of fluid content (Bai et al., 2010). The low velocity layer in the shallower part of the crust (such as depth 17.5 km) may play a role of decoupling the upper crust from the mid-and lower crust and the low velocity layer in the lower crust may decouple the crust from the upper mantle. If both low velocity layers exist, they can act as upper and lower sliding planes for the crustal materials to move eastward from the Tibetan Plateau (Teng et al., 2008). Jamieson et al. (2006) found that simply adding a weakened layer in the upper crust to a channel-flow model allowed the model to reproduce several geological observations. Our tomography results support the existence of low velocity zones in the shallower part of the crust.
The 2008 Wenchuan Ms 8.0 earthquake occurred at latitude 31° and longitude 103.4°, where there is a high velocity body seen in the 17.5 km and 25 km depth slices. This high velocity body may act as a local barrier to the channel flow so that it cannot flow to the east and north. As a result, the strain was continuously built up around the corner and the high velocity body acted as an asperity for the main shock. From a local-scale seismic tomography study around the Longmen Shan Fault using aftershocks of the 2008 Wenchuan earthquake, Pei et al. (2010) found two high velocity bodies around Wenchuan and Beichuan, associated with two large slip patches there. These two high velocity bodies act as asperities for the strain to accumulate and lead to large slip during earthquakes. These two high velocity bodies can also be identified in the depth slice of 17.5 km. There exists a relatively low velocity zone between two high velocity bodies along the Longmen Shan fault, where the aftershocks are relatively sparse (Pei et al., 2010).

Conclusions
New three-dimensional velocity models including Vp, Vs and Vp/Vs for the southeastern margin of the Tibetan Plateau covering most of Sichuan, China, provide new insights into the geodynamics of the region. The tectonic subblocks found by modeling the GPS data (Shen et al., 2005) are associated with high velocity (or strong) bodies, surrounded by low velocity (or weak) regions. Widespread low velocity zones are found below ~20 km depth in the crust with a complicated spatial distribution. At some depths, the low velocity zones are clearly bounded by faults. Aqueous fluids may exist in the mid-crust above ~35

Acknowledgements
The research presented here was supported by the Chinese government's executive program for exploring the deep interior beneath the Chinese continent (SinoProbe-02). The work was also supported by the Air Force Research Laboratory under contract number FA8718-05-C-0016, and by the Department of Energy under contract number DE-FC52-06NA27325 and under grant number DE-FG3608GO18190.

Appendix
In this appendix, we give details of "edge" and "face" stencils used in the spherical finitedifference travel time calculation method.

Two dimensional stencils
For each edge there are 8 stencils: 4 parallel to the face and 4 perpendicular to the face. We number the nodes (t 0 , t 1 , t 2 , t 3 ), with t 0 being the unknown, t 1 , t 2 being the adjacent points and t 3 being at the opposite corner. The perpendicular stencils all have (0-2) as a common edge, and the parallel stencils all have point 0 in common and either (0-1) or (0-2) as a common edge.