Abstract
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].
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.
(a)
(b)
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.
(a)
(b)
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 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.
Acknowledgments
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.