For the 12 May 2008 Mw 7.9 Wenchuan earthquake, two imbricate faults, Beichuan fault and Pengguan fault, have ruptured simultaneously. Special attention should be paid to the point of 40 km northeast of the epicenter, in which the Xiaoyudong fault intersects the above two faults, creating a complex fault structure. Surface rupture data from field surveys and previous research of dynamics studies indicate that an important transformation may take place at the intersection. But, few studies about inversion of source rupture process have focused on this issue. We establish a multiple-segment, variable-slip, finite-fault model to reproduce the rupture process and distinguish rupture sequence. Based on the nonnegative least square method and multiple-time-window approach, the spatial and temporal distribution of slip for three rupture sequences are exhibited, using teleseismic records and coseismic displacements. The conformity between synthetic and observed teleseismic records as well as the slip value of the shallowest subfaults and the coseismic displacements is utilized to calibrate the model. The results are as follows: (1) The teleseismic records inversion alone could not distinguish different rupture sequences. However, in order to make the slip of the Hongkou and Yingxiu area coincide with the field investigation, only the Beichuan fault has a bilateral rupture on the point of intersection of Xiaoyudong fault. So the possible rupture sequence is that the earthquake started at the low dip angle part of southern Beichuan fault, and then it propagated to the Pengguan fault, which caused the rupture of Xiaoyudong fault. Then the southern part of Beichuan fault with high dip angle is triggered by the Xiaoyudong fault. (2) The coseismic displacements constraint can control the slip of subfaults near the surface and has little impact on the deeper subfaults. (3) The maximum slip on the fault is located near the Yingxiu and Beichuan area; moreover, the slip is mainly distributed at the shallow region rather than at the deep, which led to serious disasters. Meanwhile, majority of the aftershocks occur in the periphery of large slip.

1. Introduction

The May 12th 2008 Mw 7.9 Wenchuan earthquake which occurred in Wenchuan, China, caused a large number of casualties and very serious engineering damage. Field surveys showed that the earthquake ruptured the middle segment of the Longmenshan fault zone; meanwhile, two rupture spaces were approximately parallel with the northeastern Beichuan-Yingxiu fault (abbreviated as BCF) with 240 km surface rupture and the Guanxian-Jiangyou fault (abbreviated as PGF) with 70 km surface rupture [14]. The BCF was divided into two sections with different movements at the point of Gaochuan [59]. The southern part of Gaochuan (southwestern BCF), the slip on the fault, is mainly thrust, while the northern part of Gaochuan (northeastern BCF) slip type changes into strike. The research of the focal mechanism of aftershocks showed that aftershocks of the northern segment of the BCF (northern part of Gaochuan) were characterized by high dip angles [10]. The width of aftershock distribution in the northern segment of BCF was obviously narrower than that in the southern section [11], which indicated the dip angle of the northern segment to be larger than that of the southern section.

It should be noted that the Xiaoyudong fault (abbreviated as XYDF) with strike about 325° and length of 6 km [12], which was a small tear fault almost perpendicular to the PGF and BCF, located at about 40 km northeast of the epicenter. Near the intersection with the XYDF, surface rupture of PGF and BCF showed obvious dislocation and discontinuity [1, 13] (Figure 1). Surface rupture of PGF and BCF was 5 km and 7 km apart near the southern part of XYDF, while near the northern part of XYDF was 7 km and 15 km apart. Aftershocks along the XYDF demonstrated a left-lateral and strike-slip; nevertheless, a thrust slip of BCF at south of Gaochuan, which indicated that an important transformation may have taken place at the intersection. According to the result of Wang and Liu-Zeng [14], the XYDF played a positive role in linking coseismic slip transfer between the BCF and PGF, which used a quasistatic stress analysis. This raises the possibility that rupture of the XYDF caused bilateral rupture of the BCF. Inversion of the source rupture process also indicated that the BCF was initiated at the surface near the intersection of the XYDF [15, 16] although this view differed from most previous studies that have dealt with this rupture process [39, 1723]. These research studies presumed a unilateral rupture. In other words, the rupture sequence of BCF, PGF, and XYDF has been ignored, which is a valuable research issue. Furuya et al. [20] used crustal deformation data sets to develop a fault source model of Wenchuan earthquake, which considered the XYDF. But the rupture did not involve the time dimension, so the rupture sequence could not be recognized.

In this paper, we present a three-dimensional (3-D) fault model combined with an investigation of the aftershock distribution, seismic reflection profile, and surface rupture data. Based on teleseismic records as well as coseismic displacement data, we invert spatiotemporal variation of the rupture process and analyze possible rupture sequences between BCF, PGF, and XYDF. Here, two standards of the judgments are utilized to calibrate the model. One is the degree of conformity of synthetic and observed teleseismic records, and the other is the slip value of the near-surface fault and the coseismic displacements.

2. Fault Rupture Model

Study on seismogenic structure of the Wenchuan earthquake showed that from south to north, the dip angles of the Beichuan fault changed significantly, the dip increased from deep to shallow in the southwestern BCF, and the northeastern section of BCF had a greater dip angle than the southwestern [1, 6, 2427]. Based on the above features, a 3-D fault model for the Wenchuan earthquake was established (Figure 2). We parameterize this complex fault using six planes, one segment for PGF and five for the southwestern and northeastern BCF. For the southwestern BCF, from deep to shallow, the faults are labeled as BCF4, BCF3, BCF2, and BCF1, with dip angle 20°, 33°, 50°, and 65°. Northern BCF is BCF5 with dip angle 60° [21], and dip angle of PGF is 30°. BCF4 intersects BCF3 at a depth of 16.6 km, BCF2 intersects the PGF at a depth of 10 km, and BCF2 intersects BCF1 at a depth of 5.4 km. The length of aftershock distribution is obviously larger than the surface rupture and in the north and south ends exceed 50 km [6]. This suggested that there may be slipping without surface rupture at both ends. In order to find the possible rupture area in the deep region, the length of the fault is larger than the surface rupture. The PGF and southwest BCF are 132 km, northeast of BCF is 180 km, and the total length of BCF is 312 km. The strike of PGF and BCF is 224° according to the surface rupture direction [1]. The hypocenter is 30.986°N, 103.364°E, according to United States Geological Survey. The initial rupture point is on the BCF3 and the depth is 14 km (China earthquake networks center). The subfaults are 5 km long and 3 km wide, and total subfaults are 854 (Table 1).

3. Data

Two different data sets, teleseismic body waves and surface offsets, are used. The teleseismic waveform data are from the Incorporated Research Institutions for Seismology (IRIS) and are deconvolved from the instrument response, integrated to obtain displacements. We use the distinct P-wave displacement observed at 36 stations with epicentral distances between 30° and 90° (Table 2). The locations of these 36 stations are shown in Figure 3, and the uniform azimuthal coverage can limit slip on the fault. The records are band pass filtered with a zero-phase third-order Butterworth filter in the range between 0.02 and 0.5 Hz and resembled with a time step of 0.2 s. The PP-waves are not used.

Another important data set we used to constrain the slip is the coseismic displacements. Field investigation shows that surface breaks extend for about 240 km from YX to QC on BCF and about 70 km from the intersection of XYDF and PGF to HW on PGF, as the red curve shown in Figure 1. According to the results of Xu et al. [1, 3, 4], the surface rupture showed that the slip on the southwestern BCF was dominated by thrust slip; the northeastern segment is mainly the strike slip. The surface rupture of PGF is concentrated in the region between the intersection of PGF and XYDF and HW, with mainly dip slip. Since the length of the subfault of our model is 5 km, more than one observed surface offset may exist in one subfault, we averaged the measured surface slip along the 5 km length of each subfault and assigned that value to the shallowest subfault element. Figure 4 shows the total slip on BCF and PGF according to the results of Xu et al. [1, 3, 4]. Table 3 gives the surface offsets of the shallowest subfaults in our model, where the number 1 subfault is located at the southwestern tip of PGF, BCF1, and BCF5, and the numbers 26, 26, and 36 are at the northeastern terminal of the these faults, respectively. It is advantageous to incorporate the surface offset data to constrain slip on the shallowest area because we know where and how much the surface slipped [28].

4. Inversion Method

After the earthquake, the slip of the fault plane generates the seismic waves which are recorded by the stations. In order to obtain the slip values of the fault plane, the inversion method is usually utilized to achieve this goal. For the inversion, the input is the observed records. Meanwhile, the output is the slip vectors of the fault plane. We apply the nonnegative least-squares inversion method to calculate the slip on the fault according to Hartzell and Heaton [29]. The fault plane is divided into a series of subfaults with the same size. The multitime window method can simulate heterogeneity of amplitude, duration, and direction for slip in the fault, as well as the rupture velocity. Every subfault can slip in any time windows, and then a complicated rise time function can be constructed. We assume a constant rupture velocity; if the slip occurs in subsequent time windows, it is interpreted as a lower rupture velocity, otherwise equal to the set value. Here, the slip is divided into two directions, dip and strike, and the sum of these two values gives total slip and direction. In this approach, the inversion function is formed by the synthetic waveforms (matrix of ), solution vector of slip on each subfault (matrix of ), and the observed data vector (matrix of ):

The function can be solved by linear least-squares method, but the result is unstable because matrix of is an ill matrix. This leads to a phenomenon that the matrix of or has a small change, and the results of matrix show greater variations. The problem can be solved by adding damped and constrained to the above function.

Here represents a priori data covariance matrix which is used to set each record in the inversion accounts for the same weight. The matrix of for minimizing moments is obtained by letting . The smoothing constraints matrix of is used to minimize the slip between adjacent subfaults as well as adjacent windows for one subfault along strike and dip direction. The slip on the peripheral area with small values is achieved by the Matrix of .

The effect of Matrix is letting the sum of all time windows for surface subfaults to equal coseismic displacements (matrix of ) at corresponding position [28, 30]. Weights , , , and control the trade-off between satisfying these constraints and fitting the data, and the best results are obtained using the minimum value of waveform misfit, as follows:

In this expression, and denote observed and synthetic waveform data. The Green function is calculated using the ray theory method with the velocity structure model Ak135 and the filter band and sampling interval are the same as records. For Wenchuan earthquake, five isosceles triangle time windows are used; the duration is 2 s with total rise time 10 s and adjacent subfault with a delay of the same length. The dip and strike slip are 90° and 180° according to the rake angle of this earthquake.

5. Rupture Sequence

According to the fault model in this paper and previous corresponding research results [1416, 31, 32], in order to analyze the rupture sequences among BCF, PGF, and XYDF, three possible cases are used (Figure 5). In Case 1, the rupture was initiated at the low dip angle part of BCF3, and then it propagated to deep area BCF4 and to shallow area BCF2 and PGF. Rupture front propagated through BCF2 and caused the rupture of BCF1 (process 1). In Case 2, the rupture was initiated at the low dip angle part of BCF3, and then it propagated to deep area BCF4 and to shallow area PGF (process 1). The XYDF was triggered by the PGF, and the BCF1 was triggered by the XYDF. Then BCF1 and BCF2 produced bilateral rupture (process 2). In Case 3, the rupture was initiated at the low dip angle part of BCF3, then it propagated to deep area BCF4 and to shallow area BCF2 and BCF1, successively (process 1). In this sequence, the XYDF was triggered by the BCF1, and the PGF was triggered by the XYDF. Then PGF produced bilateral rupture (process 2).

6. Results

6.1. Rupture Velocity Sensitivity

The rupture velocity of Wenchuan earthquake was quite different with the range 2.5∼3.4 km/s [5, 9, 16, 17, 33, 34]. To search the optimal value, the maximum rupture velocity is 2.0∼3.4 km/s, about 0.6∼1 of the shear-wave velocity in the source region with interval 0.2 km/s, and the rupture velocity has 8 values. The teleseismic waveforms are used to study the rupture history for the three rupture sequences. The rupture velocity sensitivity for rupture 1, 2, and 3 is shown in Figure 6.

The result shows that the change of residual for the three rupture sequences is similar. The residual has little change with velocity less than 3.0 km/s and obviously increases with velocity greater than 3.0 km/s. The residuals of 2.2 km/s and 3.0 km/s are smaller. The slip of the area near Nanba on BCF5 contributes most to the end of synthetic waveforms with velocity 2.2 km/s; therefore, the slip near Nanba is very small, and the value is inconsistent with the observed data. By contrast, the slip near Nanba contributes to the synthetic records at about 60 s with large amplitude when the velocity is 3.0 km/s, so the slip is about 2.0 m and in good agreement with the observed value. Furthermore, the Wenchuan earthquake has a large rupture scale and complicated rupture process, and the range of rupture velocity may be larger. In order to simulate the rupture process, the velocity 3.0 km/s is set as the highest value. The residuals with velocity 3.0 km/s of rupture sequence 1, 2, and 3 are 0.271, 0.268, and 0.267, respectively, very close to each other. So, the teleseismic data alone are not capable of choosing the best rupture model without additional constraints. The same conclusion was obtained by Hartzell et al. [16].

6.2. Coseismic Displacements Constraints

Figure 7 shows the comparison between the surface offsets and inverted slip of the near-surface fault with velocity 3.0 km/s for three rupture sequences. On PGF, the inverted slip in the south of BL is quite small for sequences 1 and 2, which agree well with the observed value. However, the inverted slip value of sequence 3 is about 2 m and the maximum is 3.8 m, which is seriously overestimating the true value. The reason is that the rupture process of sequence 3 is different from 1 and 2 on PGF. For sequence 3, the PGF has a bilateral rupture from the point of intersection with the XYDF. Then, the south surface area of PGF has the same initial rupture time with the north area. The distance between the two regions is about 50 km, and the far field Green function is similar. For the teleseismic inversion, when the Green function, initial rupture time, and the travel time of seismic phase are close for different area, the slip cannot be distinguished. So the slip in the northern part of PGF is transferred to the southern part. The slip inversions using teleseismic data only are affected by a trade-off between rupture timing and slip location [16]. The slip near the surface area of southern part of PGF contributes to the synthetic waveforms at about 30 s, so the slip is larger. On the contrary, the slip of sequences 1 and 2 for the same area contributes to the synthetic records at the beginning of wave packet, so the slip is very small. There is a quite large surface offset between YX and HK areas on BCF. The result of sequence 2 is close to the observed value, but the inverted slip of sequences 1 and 3 is massively underestimated. At the point near HK with maximum surface offset, the value of sequence 2 is about 3/4 of the observed data. Nevertheless, the result of sequences 1 and 3 is only 1/7 of the observed value. The inversion results of the three sequences are similar, except the above two regions.

Through analysis of the difference between inversion results and coseismic displacements, it is shown that the surface area of the southern PGF of sequence 3 has a large slip value and the area from HK to YX on BCF of sequences 1 and 3 has a small slip value, which are contrary to the observation. So sequences 1 and 3 are not reasonable. For these two areas, the result of sequence 2 is in good agreement with the observation. From the surface offsets, it is considered that the rupture sequence 2 is in accord with the true rupture process of Wenchuan earthquake. For the rupture process of earthquake, the situation which the propagation is to the reverse direction of rupture is infrequent. However, the same phenomenon was proved to exist. The 1984 M6.2 Morgan Hill California earthquake, 1992 Mw 7.3 Landers earthquake, and 2010 Mw 7.2 EI Mayor-Cucapuh earthquake had the same rupture mode [3537]. This phenomenon is coincident with the numerical simulations of the rupture process [38]. The complex structure of seismogenic fault leads to such a rich rupture process.

6.3. Spatiotemporal Variation of the Rupture Process of the Wenchuan Earthquake

Comparison between observed teleseismic records and synthetic records of three sequences, as well as the slip distribution of rupture velocity 3.0 km/s, is shown in Figure 8. For the three rupture sequences, the slip is concentrated at 5 areas. The first, near the initial rupture point, the area is small but the thrust slip is large. The second, under LMS on BCF3 and BCF4, the slip block is large (A1). The slip is mainly the thrust slip on BCF3 and the right lateral strike slip on BCF4. The third, from YJS to GC on BCF1 and BCF2, there is a large slip area (A2, A2-1). At the surface area, it is mainly the right strike slip, while on the bottom area is dominated by the thrust slip. The fourth, on BCF5, the slip at surface area near BC is the large thrust slip (A4). The fifth, on BCF5 between NB and QC, the slip is thrust and strike (A5). On the other side, slip distribution is obviously different for three sequences at some areas. The first, on the high dip angle part between HK and YX of the southwestern BCF, sequences 1 and 3 have no slip. Nevertheless, for rupture sequence 2, a large slip patch is located in the region dominated by mainly the thrust slip. The second, on PGF, slip distribution is similar for sequences 1 and 2, a large slip patch is located at the bottom of the fault near the middle area, but nearly no slips exist on the southwestern PGF. In contrast, for sequence 3, the large slip patch shifts to the southwestern PGF. The synthetic records of three rupture sequences show good agreement with the observed ones. Because of the close spatial distribution of the faults and the similar steep takeoff angles, the synthetic teleseismic records are almost same for the three rupture sequences. This leads to the fact that teleseismic data alone cannot distinguish between different rupture sequences, but by combining the result of far-field record inversion and the surface offsets, the optimal rupture scenarios can be distinguished.

In order to limit the slip of near-surface subfaults, the surface offsets are involved in inversion. The result is shown in Figure 9, the Figure 9(a) shows the slip distribution, Figure 9(b) shows the moment rate function for the complete rupture, Figure 9(c) shows the rupture process snapshots with 4 s interval, and Figure 9(d) shows the comparison of seven-month aftershock distribution and the projection of the spatial-slip distribution. Table 4 comprises the parameters of asperities. The total seismic moment is 0.950 × 1021 N·m, which is close to the result of USGS and Harvard CMT solution. The result shows that the slip along the southwestern BCF is near the initial rupture point and encompassing the area between YX and LMS (asperity A2-2), YJS and GC (asperity A2-1), and the area under the LMS (asperity A1). The corresponding moment is 0.387 × 1021 N·m, accounting for 41% of total seismic moment. The slip on PGF is located near Bailu (asperity A3), and the corresponding moment is 0.115 × 1021 N·m, 12% of the total value. Generally, the slip on the southwestern BCF and PGF is dominated by the thrust. For BCF5, the slip is located near the surface of BC (asperity A4), as well as in the area between NB and QC (asperity A5). The moment of BCF5 is 0.446 × 1021 N·m, 47% of total moment.

Our inversion result shows that the slip is mainly distributed at the depth above 15 km. At the fault plane, there are 6 asperities, which indicates that the earthquake is composed of at least 6 subevents. The slip is mainly distributed on the BCF, indicating that the BCF is the main rupture fault. According to the moment rate function and rupture process snapshot, the whole rupture process about 100 s is divided into 4 stages. 0∼12 s is the first stage, with 3% of the total seismic moment, corresponding to the slipping near the initial rupture point. The second stage, from 12 to 40 s, is the process of concentrating release of energy, corresponding to the rupture of southern segment of BCF and PGF, releasing 50% of the total seismic moment. As shown in Figure 9(c), at about 12 s, the rupture front arrives at the area under LMS, and then asperities A1 and A3 begin to rupture. At 16 s, asperities A2-1 and A2-2 begin to produce a large slip; meanwhile, there is a bilateral rupture on the BCF from the point of intersection with the XYDF. The third stage, from 40 to 60 s, corresponding to the rupture near BC, released 18% of the total moment. 60∼88 s is the last stage, which releases 25% of the total moment, corresponding to the rupture of the region south of BC. The outstanding peculiarity of the moment rate function is that one-half of the energy is released in the second stage which lasted 28 s, a short duration for an earthquake of that size. In such a short period of time, a great deal of energy is released, causing serious damage.

Through the analysis of aftershock distribution and slip region, some useful results are obtained. One is that majority of the aftershocks occurred in the periphery of large slip. It shows that most of the stress was released in the area of asperities so that no aftershocks occurred after the main shock. The pattern of aftershock location in relation to the areas of large slip is coincident with the 1999 Mw 7.6 Chi-Chi earthquake and 2015 Nepal Mw 7.9 earthquake [39, 40]. There are seven aftershocks with magnitude greater than 6.0 around the fault. It is remarkable that all the aftershocks larger than 6.0 occurred near the large slip region. Two aftershocks were located near the epicenter, and three were located at the periphery of the northernmost region. The other is that the aftershocks tend to be distributed at the area of the projection of the deep part of the fault. From the research of Jia et al. [27], the most dense aftershocks were concentrated from 14 to 16 km. The result of Tong et al. [8] and Qi et al. [22] showed the same situation. It is suggested that for the Wenchuan earthquake, the slip is mainly distributed at the shallow region, while the aftershocks mainly occurred in deep regions, which led to fateful disasters.

7. Conclusion

The rupture process of Wenchuan earthquake is complicated including an obvious change in slip direction, which suggests that several faults participated in the rupture. Especially at the area of 40 km northeast of the epicenter, the XYDF intersects the BCF and PGF, forming complex fault geometry. The surface ruptures and results of dynamic studies indicated that an important transformation might take place at the intersection. Unfortunately, few research studies about inversion of source rupture process have focused on this issue. So, based on three fault models according to the results of previous studies, we investigated the possible rupture sequences of the earthquake from kinetic aspects, using teleseismic records and coseismic displacements. The synthetic records and the inversion slip values of shallowest subfaults are utilized to calibrate the model. For the synthetic records, three rupture sequences yielded similar waveforms. However, the slip values of the shallowest subfaults near Hongkou area on BCF and southern Bailu area on PGF are markedly different. By comparing differences of the three sequences, we obtained that only the high dip angle part of the BCF is triggered by the XYDF, and then the slip that corresponds to coseismic displacement will appear between Yingxiu and Hongkou area. So the possible rupture sequence is that the rupture was initiated at the low dip angle part of BCF, and then it propagated to shallow area PGF. The XYDF was triggered by the PGF, and the high dip angle of BCF was triggered by the XYDF. Then the high dip angle of BCF produced bilateral rupture. Hartzell et al. [16] obtained a similar bilateral rupture on the BCF fault from the point of intersection with the XYDF to synthesize the second pulse of near field stations at the southwestern end of the BCF.

The preferred model features a heterogeneous slip distribution with a total of 5 mainly discrete asperities along BCF and one in the PGF. The main rupture is from Yingxiu to Qingchuan along the BCF, shorter than the length of aftershocks, but longer than the mapped surface rupture. Meanwhile, the slip tends to be distributed in the shallow region (depth 0–15 km) rather than in the deep. The results of the research of Wang et al. (2011) showed that the maximization of slip in shallow areas was caused by the accumulation of strain energy left over from past blind earthquakes that did not rupture the surface. Furthermore, the inverted slip of shallowest subfaults and the field investigation are compared. The following results can be obtained without the constraint slip near the surface and can easily be overestimated [30]; the constraint has negligible impact on the deeper subfaults. The resolution of the teleseismic data is low; in contrast, the strong-motion data with higher resolution is ideal for determining the direction of rupture from directivity [16]. Whether the inversion of near field records can accurately identify the rupture sequence needs further study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This research is supported by the National Natural Science Foundation of China (Grant no. 51378479), the Natural Science Foundation for Colleges and Universities in Jiangsu Province (17KJB130003), the Project with Science and Technology of Huai’an City (HAG201606), National Natural Science Foundation of China (Grant no. 51708245), and the Natural Science Foundation of Jiangsu Province (Grant no. BK20160426).