Table of Contents Author Guidelines Submit a Manuscript
International Journal of Geophysics
Volume 2011 (2011), Article ID 802346, 7 pages
Research Article

Analysis of the Crust Deformations before and after the 2008 Wenchuan Ms8.0 Earthquake Based on GPS Measurements

1Department of Surveying and Geo-Informatics, Tongji University, Shanghai 200092, China
2Department of Land Surveying and Geo-Informatics, The Hong Kong Polytechnic University, Kowloon, Hong Kong
3Institute of Earthquake Science, China Earthquake Administration, Beijing 100036, China
4Institute of Geodesy and Geophysics, Chinese Academy of Surveying and Mapping, Beijing 100830, China

Received 20 December 2010; Revised 18 March 2011; Accepted 19 March 2011

Academic Editor: Yun-Tai Chen

Copyright © 2011 Jicang Wu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


We use the movement velocities of GPS stations in western Sichuan province, China, to determine the distribution of horizontal strain accumulation before the 2008 Wenchuan Ms8.0 earthquake with a piecewise approximation approach and the coseismic displacements to determine the fault slips of the earthquake with an inversion analysis method. The results show that the distribution of the principal strain rates is strongly related to the active faults in the region, but along Longmenshan fault where the earthquake occurred, the strain rates are much lower than the others. The fault slip distribution shows mainly a thrusting with dextral striking, and the fault slips in the upper parts of the fault plane are in general bigger. Using the current strain accumulation rate and the released energy by the earthquake, we predict such a big earthquake in Longmenshan fault zone will happen in 460 to 1380 years.

1. Introduction

On May 12, 2008, a disastrous earthquake of magnitude Ms8.0 occurred in Wenchuan, Sichuan province, China, killing more than 69,000 people and injuring over 370,000 people [1]. Figure 1 shows geological setting in this area. The earthquake occurred at the east boundary of Tibet plateau and the west front of Sichuan basin and is believed to be caused by the collision between Indian plate and Eurasian plate. The red lines show the traces of active faults [2]. The Ms8.0 Wenchuan earthquake occurred in the Longmenshan fault zone (LMSFZ). The LMSFZ mainly consists of three parallel faults running from northeast to southwest [1]: Wenchuan-Maoxian fault (WMF), Yingxiu-Baichuan-Qingchuan fault (YBQF), and Guanxian-Anxian fault (GAF). The Wenchuan Ms8.0 earthquake occurred along over 300 km long YBQF and 40 km long GAF simultaneously [3]. There are 10 aftershocks with magnitude greater than 5.5 within three months near the main earthquake fault YBQF [1].

Figure 1: Geological setting. WMF is Wencuan-Maoxian fault, YBQF is Yingxiu-Baichuan-Qingchuan fault, GAF is Guanxian- Anxian fault, XSHF is Xianshuihe fault, ANHF is Anninghe fault, and LQSF is Longmenshan fault. The red circles represent the main- and aftershocks with magnitude greater than 6.0. The red lines show the active faults based on [2]. The inset map shows relative plate motion and the location of our study area. The black arrows show direction and magnitude of plate movements relative to the fixed Eurasian plate, and the black box shows the location of our study area, the western Sichuan province, China.

After the earthquake, the area attracted a lot of attentions in scientific community. Zhang et al. [3] studied the earthquake fault slips by using the inversion analysis of seismic waves and found that the LMSFZ has a high dip angle near surface and lower dip angle at depth to show a kind of listric shape which favors strain or energy to form great earthquakes. Meng et al. [4] analyzed the GPS measurements before the earthquake and deduced the average fault strike of LMSFZ is about 1~2 mm/a, much smaller than its adjacent active faults, such as the XSHF and ANHF with fault movements of about 3~8 mm/a. Meng et al. [4] also conducted the strain rates analysis using the GPS velocities and found that the strain rate near the LMSFZ was lower but with a high strain rate gradient. Kobayashi et al. [5] derived the LOS (Line of sight) displacements and offsets by using SAR technology and concluded that the earthquake rupture is characterized by a right-lateral motion in the northeast on the Beichuan fault, while an oblique slip consisting of right-lateral and reverse motions is suggested in the southwest. However, much study is still needed to better understand the mechanism of this disastrous earthquake using the existing welfare of information and data.

However, the derived strain rates are contaminated by the errors in GPS measurements. Their significance must be carefully assessed to distinguish between those computed strain rates generated from GPS observation errors and those of crust strain accumulation in the region. Moreover, the earthquake fault slips can also be estimated by the inversion analysis of coseismic displacements with a dislocation model in addition to the analysis of seismic waves as studied in [3]. Therefore, the objectives of this study are (1) to understand better the crustal deformation pattern before the earthquake by studying significant strain accumulation in the LMSFZ and its adjacent regions and (2) to understand better the fault slip distribution of the earthquake by using the inversion analysis of the coseismic displacements.

2. GPS Surveys in the Region

Under a research project “Dynamic Process and Strong Earthquake Prediction of Block Boundaries,” 61 campaign-based and 4 continuously operating GPS stations in western Sichuan region were established and surveyed in 2005 by China Earthquake Administration (CEA). In addition, 22 GPS stations of the Crustal Movement Observation Network of China (CMONOC) were surveyed twice before 2005. All of these 83 GPS stations were resurveyed in 2005, 2006, and 2007. The surveys were carried out using Trimble 5700 GPS receivers with choke ring antennas and lasted for 22–24 hours each day for 3 consecutive days in each survey campaign. GAMIT and GLOBK software were used for the determination of the velocities of GPS stations in ITRF2005 frame based on the campaign mode and continuous GPS measurements. The rotation of the Eurasian plate was eliminated by carefully choosing Eurasia-fixed frame [4]. Some GPS networks established and surveyed in the same accuracy level by local organizations are also available, totaling 176 GPS stations in this area. Before the earthquake, the area was monitored with multiple GPS campaigns. These results can be used to derive the movement velocities at the GPS stations (see electronic supplement of [4]). Figure 2(a) shows the movement velocities of these GPS stations and their 95% confidence error ellipses before the earthquake.

Figure 2: Movement observations of GPS stations. (a) The movement velocities of GPS stations before Wenchuan earthquake (In ITRF2005 with Eurasian plate fixed). (b) Horizontal coseismic displacements at GPS stations; red star is the locations of Wenchuan Ms8.0 earthquake.

Upon the earthquake, the State Bureau of Surveying and Mapping (SBSM) and CEA conducted GPS surveys in the area. The coseismic displacements at the GPS stations are obtained from the results of the continuous and campaign GPS surveys before and after the earthquake. Dual-frequency GPS receivers with choke ring antenna were used, and each station was occupied for two observation sessions of 23 hours each. The GPS surveys were completed on June 2, 2008. Also CEA resurveyed at their GPS stations in the area. These GPS observations were processed with GAMIT and GLOBK software to estimate the coseismic displacements. Figure 2(b) shows the horizontal coseismic displacements, where the red dots are the locations of the Wenchuan Ms8.0 earthquake and the aftershocks.

3. Methodology for Analyzing the Crustal Deformations

3.1. Determination of Strain Rates in a Piecewise Approximation Approach

The approach of piece-wise approximation divides the area of interest into a number of triangles with GPS stations as their vertexes. The Delaunay method is applied to form triangle network by using the GMT software [6]. Within each triangle, a homogenous strain is assumed. The movement velocities of the GPS stations at the three vertexes of each triangle are used to calculate the principle strain rates.

Let the velocities and their standard deviations of each station of a triangle be where and are the eastern and northern velocity components, and , , and are the standard deviations and covariance, respectively. Under the assumption of homogeneous strain in each triangle, the displacement velocities and the strain rates are related by [7] where , , and are the components of strain rate tensor and is the differential rotation rate; , and are the coordinate difference in the eastern and northern directions, respectively, between stations and . For each triangle, we can form four independent equations as (2) involving four unknowns. The unknowns , , and and can then be solved. The azimuth of the principal strain rate is computable as [8] and the maximum and minimum strain rates, respectively, are Using the error propagation law with given variances and covariance of movement velocities of each GPS station, the covariance matrix of the estimated unknown can be obtained.

We then test the significance of the estimated strain tensor. The null hypothesis and alternative hypothesis are formulated as where . Let be the variance-covariance matrix of the strain tensor, which is the upper left submatrix of . The statistic is which follows a chi-square distribution with degrees of freedom 3. Selecting significance level , we can calculate the boundary value . If , we reject null hypothesis and accept alternative hypothesis, meaning the strain tensor in the triangle is significant. To better visualize the strain accumulation status in the area of study, we calculate the principal strains only for those triangles with significant estimated strain tensors.

3.2. Fault Slip Distribution by the Inversion Analysis of Coseismic Displacements

Okada [9] gives a forward model for calculating ground displacements caused by fault slips on a rectangular fault plane in a homogeneous elastic half-space. Generally, the ground displacement component is a function of fault geometry and dislocation components (fault slips) expressed as where is the horizontal coordinates of the middle point of the upper side of the fault plane, is the strike of the fault plane, , , , and are half-length, -width, and -depth of upper side and dip angle of the fault plane, respectively, and , , and are the strike, thrust, and tensile dislocation components (fault slips) on the fault plane. In (7), the ground displacement component is the observed quantities, and other quantities in the left side of the equation are the unknowns to be estimated. In practice, several quantities in the left side of equation, like the fault location parameters , , and and the fault geometry parameters , , , and , can be determined from other means, for example, geological investigations and the earthquake mechanism solutions [1, 3]. They can be treated as either known quantities or unknown parameters with prior information on their values and uncertainties in the least-squares method (also called as Bayesian method). For the former case, the number of unknown parameters will be greatly reduced, resulting in better solution for fault slip parameters; for the latter case, these unknown parameters are given prior values and larger weights in the inversion process [10]. The inversion analysis of the coseismic displacements for fault slips is done based on (7) using a Bayesian method. For the detail on the method, readers can refer to [11, 12].

4. Results and Discussion

4.1. Determination of Strain Accumulation before the Earthquake

The methodology discussed in Section 3.1 and the movement velocities of 176 GPS stations obtained by GPS surveys before the earthquake (Figure 2(a)) were used to calculate strain rate tensor for each triangle. They were statistically tested with significance level 0.1%. For those triangles with significant strain accumulation, the principal strain rates were calculated. Figure 3(a) shows the significant principal strain rates with double outward arrows representing extension and inward arrows representing compression. Figure 3(b) shows total strain rates which look a little bit mess. In this figure, red lines represent active faults, red star indicates the Wenchuan Ms8.0 earthquake, and red dots indicate aftershocks with magnitude greater than 6.0.

Figure 3: Principal strain rates distribution before Wenchuan earthquake. Double arrows show the principal strain rate labeled at the center of individual triangles. Red circles represent the main- and aftershocks with magnitude greater than 6.0. The Delaunay triangle network is formulated with GPS points as their vertexes. The red lines show the active faults based on [2]. (Only the principal strain rates passed the drawn in (a).)

We can see from Figure 3(a) that only small parts of triangles have significant principal strain rates, and most of them are in compression mode. Only two triangles have extension principal strain rates, one at the junction of XSHF and LMSFZ, near to Garze, and the other in the Chengdu basin, related to LQSF. This agrees with the compression tectonic environment of Tibet [13, 14]. The location of these significant principal strain rates is related to active fault traces, and generally the maximum principal strain rate axis is perpendicular to the related active fault strike. For fault YBQF where the Wenchuan Ms8.0 earthquake occurred, most of the principal strain rates are not significant except at its southeast end with compression strain rates, which are near to the initial rupture segment of the Ms8.0 Wenchuan earthquake. The results suggest that the areas with significant strain accumulation are hazardous and need great attention and detailed study and monitoring for earthquake precursors.

4.2. Determination of Fault Slips by the Inversion Analysis of the Coseismic Displacements

The coseismic displacements in Figure 2(b) and the methodology discussed in Section 3.2 are used to determine the earthquake fault slips. Based on Xu et al. [1] and Zhang et al. [3], we assumed that the source of the Wenchuan earthquake is a rectangular fault plane with length 350 km, width 50 km, depth of upper side 3 km, and dip angle 39°, and the azimuth of the fault strike 230°. We divided the fault plane into square grids of 10 km × 10 km each. In our inversion analysis, we fixed the fault geometry to be above values, and only the fault slips on each grid are parameters to be determined. The a priori value for U1(strike slip), U2(dip slip), and U3(tensile slip) is set to be −0.5 m ± 1 m, 0.5 m ± 1 m, and −0.5 m ± 1 m, respectively. The dislocation components (fault slips) on each grids were determined by the Bayesian inversion approach using coseismic displacements (Figure 2(b)). The estimated fault slip distributions are given in Figure 4, where the top left side shows the fault slip distribution in the horizontal geodetic coordinate system, and the lower part on the fault plane. The top right side of the figure gives the data fitting residuals. We can see that most of the residuals are less than 10 mm, indicating that the estimated fault slips fit well with the observed coseismic displacements.

Figure 4: Fault slip distribution by inversion of horizontal coseismic displacements. Top left side is fault slip distribution in horizontal geodetic coordinate system, and bottom is fault slip distribution in the fault plane, and the circles at the heads of arrows show the 95% uncertainty error ellipses. Top right side is the histogram of data residuals.

Figure 4 shows that the fault slips are mainly in the depth less than 20 km. Generally, the distribution of fault slips indicates that the earthquake fault movement is a thrust with dextral strike. In the southwest part, the thrust movement prevails, while in the northeast part, the dextral strike prevails. The maximum thrust is 5.6 m, and the maximum dextral strike is 5.2 m.

Using the obtained fault slip distribution, the seismic moment can be calculated based on [15] as where is the shear modulus of the rock in the area, is the th fault grid area, and is the magnitude of fault slip on the th fault grid with . Taking and using and the estimated fault slips, the total seismic moment released can be calculated from (8) as (Joule). The seismic moment can further be used to determine the earthquake magnitude based on [15] as Substituting into (9), we get , very close to the published earthquake magnitude of 8.0 based on the seismic waves. If we assume that the fault deficit accumulation rate is uniform at the current rate of 1 to 3 mm/a [3], the energy accumulation rate can be calculated using (8) in a similar way, resulting in to , respectively. Dividing the total seismic moment by the estimated energy accumulation rate, we predict earthquake cycle to be 460  to 1380 .

5. Conclusion

The proposed approach can provide the distribution of significant principal strain accumulation using the velocities of GPS stations, which is important information for studying crustal deformations and can help identify the areas of geological hazards. There is no significant strain accumulation in most parts of LMSFZ before Wenchuan Ms8.0 earthquake comparing to nearby active faults, such as XSHF and AZHF; this could be an indicator of earthquake potential. The fault slip distribution estimated from the coseismic displacements shows that Wenchuan Ms8.0 earthquake is mainly a thrusting in southwest part and with a significant dextral striking in the northeast part, and the fault rupture is mainly distributed at the depth less than 20 km. Using the estimated fault slip distribution, the energy released is approximately estimated, and magnitude of the earthquake is calculated, which agree well with the results derived from seismic waves. Under our assumption of uniform fault deficit accumulation rate of 1 to 3 mm/a, a large earthquake will take place in 460–1380 years.


This study was supported by the Hong Kong Research Grant Council (PolyU A/C: B-Q02B) and by China National Science Fund (no. 40674004). The figures are prepared using the GMT software.


  1. Z. Xu, S. Ji, H. Li, L. Hou, X. Fu, and Z. Cai, “Uplift of the Longmen Shan range and the Wenchuan earthquake,” Episodes, vol. 31, no. 3, pp. 291–301, 2008. View at Google Scholar · View at Scopus
  2. X. Ma, Ed., Lithospheric Dynamics Atlas of China, China Cartographic Publishing House, Beijing, China, 1989.
  3. P. Z. Zhang, X. W. Xu, X. Z. Wen, and Y. K. Ran, “Slip rates and recurrence intervals of the Longmen Shan active fault zone, and tectonic implications for the mechanism of the May 12 Wenchuan earthquake, 2008, Sichuan, China,” Chinese Journal of Geophysics, vol. 51, no. 4, pp. 1066–1073, 2008. View at Google Scholar · View at Scopus
  4. G. J. Meng, J. W. Ren, M. Wang et al., “Crustal deformation in western Sichuanregion and implications for 12 May 2008 Ms 8.0 earthquake,” Geochemistry, Geophysics, Geosystems, vol. 9, no. 11, Article ID Q11007, 2008. View at Publisher · View at Google Scholar · View at Scopus
  5. T. Kobayashi, Y. Takada, M. Furuya, and M. Murakami, “Locations and types of ruptures involved in the 2008 Sichuan earthquake inferred from SAR image matching,” Geophysical Research Letters, vol. 36, no. 9, Article ID L09304, 2009. View at Publisher · View at Google Scholar · View at Scopus
  6. P. Wessel and W. H. F. Smith, “New improved version of the generic mapping tools released,” Eos, Transactions American Geophysical Union, vol. 79, no. 47, p. 579, 1998. View at Google Scholar
  7. J. C. Wu, H. W. Tang, Y. Q. Chen, and Y. X. Li, “The current strain distribution in the North China Basin of eastern China by least-squares collocation,” Journal of Geodynamics, vol. 41, no. 5, pp. 462–470, 2006. View at Publisher · View at Google Scholar · View at Scopus
  8. A. C. Eringen, Mechanics of Continua, Robert E Krieger Publishing Company, Malabar, Fla, USA, 1980.
  9. Y. Okada, “Surface deformation due to shear and tensile faults in a half-space,” Bulletin of the Seismological Society of America, vol. 75, no. 4, pp. 1135–1154, 1985. View at Google Scholar
  10. J. C. Wu, Inverse Analysis of Crust Deformation Measurements, Ph.D. thesis, Hong Kong Polytechnic University, Hong Kong, 1999.
  11. D. D. Jackson and M. Matsu'ura, “A Bayesian approach to nonlinear inversion,” Journal of Geophysical Research, vol. 90, no. 1, pp. 581–591, 1985. View at Google Scholar · View at Scopus
  12. J. C. Wu, H. W. Tang, Y. Q. Chen, and Y. X. Li, “Inversion of GPS measurements for a layer of negative dislocation distribution in north China,” Journal of Geophysical Research, vol. 108, no. B10, 2003. View at Publisher · View at Google Scholar · View at Scopus
  13. Y. Zhang, W. P. Feng, L. S. Xu, C. H. Zhou, and Y. T. Chen, “Spatio-temporal rupture process of the 2008 great Wenchuan earthquake,” Science in China Series D, vol. 52, no. 2, pp. 145–154, 2009. View at Publisher · View at Google Scholar · View at Scopus
  14. M. Taylor and A. Yin, “Active structures of the Himalayan-Tibetan orogen and their relationships to earthquake distribution, contemporary strain field, and Cenozoic volcanism,” Geosphere, vol. 5, no. 3, pp. 199–214, 2009. View at Publisher · View at Google Scholar · View at Scopus
  15. D. Turcotte and G. Schubert, Geodynamics, Cambridge University Press, Cambridge, UK, 2nd edition, 2002.