In order to reproduce the rupture history of the 2008 Mw8.0 Wenchuan earthquake, the teleseismic and strong-motion records are adopted. Based on a multiple-segment, variable-slip model, the finite fault inversion method is utilized to recover the rupture process. The results are as follows: (1) the rupture duration of the Wenchuan earthquake is about 100 s, and the released seismic moment is 1.24 × 1021 N·m, equal to the moment magnitude Mw8.0. There are 5 asperities on the fault plane, indicating that the earthquake is composed of at least 5 subevents. (2) The slip is mainly distributed on the Beichuan fault, indicating that the Beichuan fault is the main rupture fault. On the southern part of the Beichuan fault, the dislocation underside the Longmenshan area and Hongkou-Yingxiu near-surface area is dominated by thrust, and the maximum slip is 11.8 m. Slip between the Yuejiashan and Qingping area is dominated by thrust. On the northern part of the Beichuan fault, the area under Beichuan is dominated by thrust, the slip under Nanba is thrust and strike, near Qingchuan, the slip turns into the strike slip, and the maximum slip is 13.1 m. The dislocation under Bailu is also dominated by thrust, with maximum slip 8.9 m. (3) The rupture of the Wenchuan earthquake is mainly a unilateral rupture to the northeast. The rupture started at the low dip angle part of the Beichuan fault, and after 3 s, it propagated to the Pengguan fault. After 10 s, the largest asperity under Longmenshan in the south section of the Beichuan fault began to break, lasting for about 24 s. Then, the Xiaoyudong fault was triggered by the Pengguan fault, and the bilateral rupture of the high dip angle part of the Beichuan fault started at about 6 s. South section of the Beichuan fault began to break at about 35 s, and at 43 s, 63 s, and 80 s, the rupture extended to Beichuan, Nanba, and Qingchuan areas.

1. Introduction

The 2008 Mw8.0 Wenchuan earthquake (also reported as Mw7.9 by some researchers [13]. In this study, the result of released seismic moment is 1.24 × 1021 N·m which is equal to moment magnitude Mw8.0 based on the relation between seismic moment and moment magnitude defined by Hanks and Kanamori [4]) consists of multiple fault segments and produced a complex system of faults along the edge of the Tibetan Plateau in Longmenshan, causing a large number of casualties and property losses. The Longmenshan fault zone, where the seismogenic fault is located, is a complex structure, which consists of four reverse faults [5]. In the process of propagation from the epicenter to the northeast, the rupture passes through several stages, with a fault length of 300 km. The two main rupture zones, which are close to parallel in space, have surface fractures in the Beichuan-Yingxiu fault (BCF) and Pengguan fault (PGF). The former is 240 km long, while the latter is 72 km long [57]. It is worth noting that there is also a Xiaoyudong fault (XYDF) about 30 km northeast of the epicenter [8], which is sandwiched between the BCF and PGF, about 6 km long (Figure 1). The complexity of the fault is particularly evident in the XYDF area, where obvious bends and discontinuity of the BCF and PGF existed at the intersection. At the southwest of the intersection, the distance between two faults was 5∼7 km, while at the northeast, it was 7∼15 km [5]. The result of aftershock relocation shows that there is a nearly vertical aftershock distribution belt about 50∼60 km long in the northwest extension direction of the XYDF [8, 9]. The focal mechanism solution of aftershock shows that it is a left-handed strike slip [10], which is obviously different from the thrust quality of aftershock in the south of the BCF. All these phenomena indicate that some important transformations may have taken place at the intersection.

For the Wenchuan earthquake, many researchers have given the slip distribution [13, 1122]. Different data are used in inversion, and the information of the fault slip distribution revealed is different. The inversion based on single data has both advantages and disadvantages. For the teleseismic data with an epicenter distance of 30°∼90°, the propagation path of the seismic phase is simple, the time of calculation of Green’s function is short, the observed record can be obtained quickly after an earthquake, and it can reflect the slip of the deep fault area. However, due to the low resolution of teleseismic Green’s function, it is not easy to obtain the fine rupture process for the Wenchuan earthquake, which is very complex. However, the strong-motion data can reflect the details of the fault rupture process and recognize the slip in the shallow area of the fault, but they are not sensitive to the slip in the deep area, and the calculation of Green’s function is greatly affected by the velocity structure. In order to overcome the shortcomings of the single-data inversion, the rupture process of the Wenchuan earthquake is inverted by combining teleseismic and strong-motion data.

2. Data

The teleseismic records are selected from the Incorporated Research Institutions for Seismology (IRIS) (http://ds.iris.edu/wilber3). Vertical P wave (BHZ) records of 36 stations with uniform coverage of azimuth in the range of epicenter distance 30°∼90° are selected. Figure 2 shows the station location. The azimuthal coverage of the station is relatively uniform, which is very important for inversion of a reasonable fracture process. To ensure the accuracy of arrival time of seismic phases, the third-order Butterworth band-pass with zero-phase shifting filter is adopted, and the frequency band is 0.02∼0.5 Hz. The record resampling interval is 0.2 s, and the calculated far-field Green’s function does not contain the information after the PP seismic phase, so the waveform between P and PP seismic phases is taken as the inversion data. The arrival time of P and PP phases is the time when the waveform at the source propagates to the station. Due to different epicenter distance, arrival times between P and PP seismic phases are different, and the length of each record is different. The waveforms of stations with different azimuth have high similarity. Most of the waveforms can be divided into four sections. Here is the waveform of station KMBO (Figure 3). The first segment lasts for 10 s from 0 to 10 s, with small amplitude. The second segment lasts for 50 s from 10 to 58 s, with the largest amplitude. The third segment lasts for 20 s from 58 to 78 s, with small amplitude. The fourth segment is after 78 s, with long duration and small amplitude.

Three component acceleration records of 43 stations within 163 km of fault distance were selected from the “National Strong-Motion Network Center” of Institute of Engineering Mechanics, China Earthquake Administration (Figure 1). The acceleration record is integrated into the velocity record, and the band-pass nonphase shift filter is used for 0.08∼0.5 Hz. The interval of resampling is 0.2 s, and the length is 120 s. The full waveform is taken as the inversion data, including the P wave, S wave, and surface wave.

Except for 8 bedrock stations such as 51WCW, 51BXZ, 51MXT, 51SPT, 62WIX, 51PXZ, 51XJL, and 51CDZ, the rest are all soil stations, most of which are located in class II sites with the average shear wave velocity in the range of 500∼800 m/s. According to drilling information given by Yu [23], predominant frequency of the site of these stations can be calculated by f = Vs/4H. Vs is the shear wave velocity, and H is the thickness of the soil layer. The predominant frequencies are all above 2.0 Hz, which is higher than the upper limit of 0.5 Hz, so the influence of topsoil cannot be considered in the calculation of Green’s function. Station azimuth coverage is uniform in the PGF and the south section of the BCF. There are 15 stations with fault distance less than 60 km and some stations such as 51WCW, 51SFB, and 51MXN, which are close to the fault. Station azimuth coverage is relatively poor in BC5. There are only 51JYH, 51JYC, 51JYD, and 51GYZ with fault distance less than 60 km, and there are no stations at the northeast end of the fault. The distribution of stations in the south section of the Beichuan fault is better than that in the north section, and the inversion results in the north section may not be well controlled.

Meanwhile, the coseismic displacements are utilized to limit slip values of the shallowest subfaults. According to the field investigation of Xu et al. [5], the values of surface offsets of the shallowest subfaults are used in our inversion [21].

3. Fault Rupture Model and Rupture Sequence

Wenchuan earthquake occurred on the imbricated curved fault. Some research studies have analyzed the seismogenic structure [5, 17, 2427]. Based on the research results, the change of dip angle in the south section of the Beichuan fault and the difference of dip angle between the north and south section are considered. A 3D fault model is established with fault strike 224°, including the BCF and PGF (Figure 4). The southern part of the Beichuan fault has a dip angle of 20° (BCF4), 33° (BCF3), 50° (BCF2), and 65° (BCF1) from deep to shallow. The deepest point of BCF4 is 21.7 km, BCF4 and BCF3 intersect at a depth of 16.6 km, BCF2 and PGF intersect at a depth of 10 km, the PGF dip angle is the same as BCF3, BCF2 and BCF1 intersect at a depth of 5.4 km, and BCF1 and PGF intersect with ground surfaces at a distance of 9.06 km.

The length of the aftershock distribution is obviously larger than the scale of the surface rupture, and the length of the north and south ends exceeds about 50 km. It indicates that there may be slip at the north and south ends without the surface rupture. In order to reflect the possible rupture area, the fault length is taken as greater than the surface rupture length. PGF and the south section of the Beichuan fault are 132 km long, while the north section of the Beichuan fault is 180 km long. The total length of the Beichuan fault is 312 km.

A single fault model is easy to determine the fracture mode. For the Wenchuan earthquake, the Longmenshan fault zone is complex, and the rupture sequence between the faults has not reached the same conclusion [1, 11, 20, 21, 28, 29]. Yin et al. [20, 21] adopted three possible cases to determine rupture sequences among BCF, PGF, and XYDF and obtained the possible rupture sequence as follows: the rupture starts at the low dip angle in the south part of the BCF, which causes the rupture of the PGF and then causes the rupture of the XYDF which triggers the bilateral rupture of the high dip angle in the south part of the BCF. The north section of the BCF is triggered by the south section of the BCF. In this paper, we use the same rupture pattern.

4. Weighting Factor

Combining the teleseismic and strong-motion records for inversion, the number of equations is 108,099, which is about 13 times bigger than the parameters. We need to determine the weight of different data in the inversion. There are different methods for selecting the weight. Hartzell and Heaton [30] took into account the whole record of teleseismic records has a larger amplitude after normalization, while a larger amplitude only accounts for a part of the whole record after normalization for strong-motion records, so the maximum value of far-field record is adjusted to 0.5 during joint inversion. Hartzell et al. [1] used the teleseismic record, near-field record, and GPS data to perform joint inversion of the Wenchuan earthquake rupture process, setting different data with the same weight. The data used in this paper are similar to those used in [1]; we refer to their calculation method. First, the sum of the absolute values of the data in Green’s function matrix Sum (Gall) is calculated, and then the sum of the absolute values of Green’s functions for teleseismic and strongmotion records Gsi is calculated. The weighting factor is Wi = Sum (Gall)/Gsi. The weighting factor for teleseismic and strong-motion data is 9.2 and 1.1. The weight factor for teleseismic data is about 9 times of the weight factor for strong-motion data, and the inversion results are more consistent with the results of teleseismic records. Considering the number of data used in this paper, the strong-motion records are about 4 times of the teleseismic records. The strong-motion records should occupy the main weight in the inversion. Meanwhile, in the inversion, we find that the residuals of strong-motion data are more sensitive to weighting factors than teleseismic data. If we set the weighting factor for teleseismic and strong-motion data as 9.2 and 1.1, the residual of strong-motion data is large. The misfit of strong-motion data is obvious. However, with the decrease of the weighting factor for teleseismic data, the decrease of misfit for teleseismic data is slow. In order to get acceptable misfits for the strong-motion and teleseismic data, the factor of strong-motion data is improved, and the factor of teleseismic data is reduced. The weight factors of teleseismic records and strong-motion records are 2.1 and 1.9 by trial and error.

5. Methods and Parameters

The nonnegative least square method and multitime window method are used to invert the fault slip history. Refer to Hartzell and Heaton [30] for details. The fault plane is divided into uniform subfaults. Each subfault contains a series of windows. Every subfault can slip in any time windows. In order to recover the slip direction of these subfaults, the direction of each subfault is divided into two mutually vertical directions. The vector sum of the two vertical directions’ amplitudes is the final slip value. According to the CMT solution of Harvard University, the slip angle of this earthquake is about 138°. So, two orthogonal directions are 90° and 180°. Slip values of each subfault for the two mutually vertical directions are the solution vector x. Green’s functions are matrix A. The observed records are matrix b. To obtain a stable solution, the linear stability constraints are added into the least square method. The inversion equation is as follows:

represents a apriori data covariance matrix which can be used to set each record has the same weight in the inversion. Matrices M and S are utilized to minimize the moment and slip between adjacent subfaults as well as adjacent windows. Matrix B can restrain the slip on the peripheral area. Matrix F is used to let the slip near the surface to equal the coseismic displacements. The purpose of weights λ1, λ2, λ3, and λ4 are to control tradeoff between satisfying the above constraints and fitting the data. Here, we use a trial-and-error approach to determine the constraint weight. Based on equation (2), the best results are obtained. Here, xi and yi denote observed and synthetic waveform data, and n is the number of total of stations.

The form of the slip time function is an isosceles triangle, lasting for 2.0 s. The subfault contains 5 time windows, the rupture delay of adjacent windows is 2.0 s, and the maximum rise time of the subfault is 10 s. The size of the subfault along strike and dip directions is 5 km and 3 km. For teleseismic records, Green’s function is calculated using the ray theory method [31] based on velocity structure model Ak135 [32]. The discrete wave number method [33] is used to calculate Green’s function of strong-motion data. The thickness and medium velocity structures of the Tibet Plateau and Sichuan Basin on both sides of the Longmenshan fault zone are different. A 1D crustal velocity model (Table 1) is selected according to Hartzell et al. [1]. Here, the teleseismic records are used to determine the suitable rupture velocity. The rupture velocity is 2.0∼3.4 km/s with interval 0.2 km/s. Through the comparison of these results, the reasonable rupture velocity is 3.0 km/s [20].

6. Results

6.1. Slip Distribution

The inversion results are shown in Figures 59 and Table 2. The residuals of strong-motion and teleseismic records are 0.50 and 0.34. The main waveform characteristics of these observation records are well explained, and 63% and 70% of the strong-motion and teleseismic stations have a correlation coefficient greater than 0.7 and 0.8. The correlation coefficient is as follows. X and Y are the observed and synthetic records. Cov is the covariance. Var is the variance.

The dislocation distribution shows that there are five asperities on the fault plane marked on Figure 5 by A1, A2, A3, A4, and A5. An asperity can be defined as enclosing fault elements that have a slip 1.5 times larger than the average slip over the fault [34]. The slip of the south section of the BCF is mainly concentrated in two regions. The first asperity A1 is located in the low dip area under Longmenshan (LMS), with the strike slip dislocation in the south and thrust dislocation in the north and shallow areas. The maximum value of the slip is 11.8 m, and the size is about 1400 km2. The second area A2 is in the high dip angle from Yingxiu (YX) to Hongkou (HK), which is dominated by the thrust slip. The peak value of the slip is about 4.5 m, and the size is 450 km2. Besides, there is also some slip near the surface of Yuejiashan (YJS) and at the bottom of the north side of the fault. The moment is 0.47 × 1021 N·m, accounting for 40% of the total moment. The rupture duration of the south section of the BCF is about 40 s. For the high dip angle faults BCF1 and BCF2, the discreteness of strong-motion results (Figure 10(b)) is reduced, and the slip distribution of joint inversion (Figure 5) is much smaller than that of teleseismic records (Figure 11(b)). The resolution of joint inversion is higher than that of single data.

The slip of PGF A3 is concentrated in the bottom area of Bailu (BL) and the near-surface area on both sides of BL. Generally, it is mainly thrust dislocation. The peak slip value is 8.9 m, and the size is 600 km2. The rupture lasted for 29 s. The released seismic moment is 0.14 × 1021 N·m, accounting for 11% of the total moment. The slip of BCF5 is also distributed in two areas. A4 is concentrated in the area under Beichuan (BC); A5 is located from Nanba (NB) to Qingchuan (QC). The area under BC is dominated by thrust, with a peak slip value of 10.5 m. The slip under NB is thrust and strike, and near QC, the slip turns into strike slip. The rupture duration is about 60 s. The released seismic moment is 0.75 × 1021 N·m, accounting for 49% of the total seismic moment.

6.2. Rupture Evolution

According to the moment release rate (Figure 8) and the rupture process (Figure 9), the rupture of this earthquake can be divided into four stages. In the first stage, 3% of the total seismic moment is released from 0 to 12 s, corresponding to the slip near the initial rupture point in the south part of the BCF. The second stage, from 12 to 40 s, is the stage of energy concentrated release, corresponding to the rupture of the south section of the BCF and PGF, releasing 48% of the total seismic moment. In the third stage, from 40 to 60 s, 18% of the total seismic moment was released corresponding to the rupture near BC. The fourth stage is from 60 to 88 s, corresponding to the rupture from NB to QC, releasing 25% of the total seismic moment.

According to the rupture process snapshots with 4 s interval (Figure 9), we can get the propagation process of this earthquake. The rupture started on BCF3, and the rupture front reached the lower part of LMS at about 12 s. At this time, the main slip area on the south part of the BCF began to break. At 16 s, the slip on the south part of the BCF and PGF started to produce large slip, while the high dip angle part of BCF1 and BCF2 also began to break the bilateral rupture at about 20 s. On the south section of the BCF, the rupture front reached QP and GC at about 30 s and 35 s. On BCF1 and BCF2, at about 20 s and 28 s, the rupture front reached HK and YX. At about 44 s, the rupture front reached BC. At about 64 s, the rupture front reached NB.

6.3. Comparison with Other Research Results

Through the comparison of inversion results of different fault models, it can be seen that although the fault models are different, the areas with larger slip on the fault plane are located in YX, YJS, and BC areas [13, 1119, 22]. The forms of dislocations are different. Compared with most studies, the rupture modes of BCF3, BCF4, PGF, and BCF5 in this paper are similar, and the general characteristics of slip distribution on these fault planes are similar, but the location, maximum slip, and dislocation form of the slip are different. However, the slip distribution in the area from HK to YX is obviously different. In teleseismic record inversion, most studies obtained that the high dip angle part of the south section of the BCF is triggered by the low dip angle part. The rupture propagates from south to north, and there is no slip in this area [14, 22]. The results show that only a bilateral rupture on the BCF from the point of intersection with the XYDF, the surface value between HK and YX is consistent with the observation results. Wang et al. [2], Nakamura et al. [15], and Hartzell et al. [1] obtained that there is a near-surface slip in this area, but the slip type and amount are quite different from the observation results. By comparing the results of this paper with the results of geodetic data, the strengths and limitations of teleseismic data, strong-motion data, and geodetic data are clearly demonstrated. Geodetic data can better locate the slip areas in the shallow part. The situation can be seen from the research results of Tong et al. [18], Xu et al. [19], Hartzell et al. [1], and Chen et al. [35]. However, the seismic data can distinguish shallow from deep slip.

The comparison between the synthetic and observation records of stations located behind the rupture direction shows that the bilateral rupture of BCF1 and BCF2 is better than the unilateral rupture of BCF1 and BCF2. The amplitude and phase information of the second band of the synthetic record are closer to the observation record. Hartzell et al. [1] combined the strong-motion record, teleseismic record, and GPS data to inverse the rupture history of this earthquake and obtained the bilateral rupture. From the results of this paper, the inversion results of teleseismic and strong-motion records also indicate the existence of the bilateral rupture. Therefore, in order to satisfy the second wave band of synthetic records of these stations located behind the rupture direction (Figure 1, stations with blue triangles), the BCF should have reverse rupture. It may be related to the complex seismogenic structure of the Wenchuan earthquake [29]. This shows the complexity of the Wenchuan earthquake rupture.

The strong-motion record and teleseismic record utilized by Hartzell et al. [1] are similar to the record used in this paper. The inversion results of only strong-motion data or teleseismic data as well as the results of this article are shown in Figures 10 and 11. Compared with the single strong-motion or teleseismic data inversion result, the area and dislocation of slip under LMS on BCF3 and BCF4 are similar to the result of the teleseismic record inversion. The distribution and dislocation form of the slip from YJS to the northern end of the fault on BCF3 and BCF4 are similar to the results of the strong-motion data inversion. The dislocations in the northern half of the PGF synthesize the results of strong-motion and teleseismic data inversion results. For the southern half of the PGF, the dislocation obtained by the joint inversion is significantly lower than that of the strong-motion data inversion result. The combination of these two records can better limit the slip distribution in this region. In addition, the inversion results of strong-motion records on BCF5 show that the dislocation distribution is very discrete, and it is not easy to distinguish the asperity (Figure 10). However, combining the results of the two data inversions, we found that the slip concentrated area appeared in the lower side of Beichuan, and the form of dislocation was similar to that of the teleseismic inversion. Compared with the inversion results of strong-motion records, the joint inversion results improve the recognition ability of the dislocation distribution near Beichuan on BCF5, but they are still not suitable to identify the location of asperity in the north of Beichuan.

7. Conclusion

In order to overcome the defects of single data with low resolution, combining the teleseismic and strong-motion data, the rupture process of the Wenchuan earthquake is inversed. Combining teleseismic records and strong-motion records has improved the identification ability of the slip distribution in the high dip angle area on the southwestern part of the BCF, the southern part of the PGF, and the area around BC on the northeastern part of the BCF.

The areas of slip concentration on the fault plane are as follows: low dip area under LMS in the south section of the BCF, near-surface area from HK to YX, near-surface area from YJS to QP, bottom area of the fault under Bailu in the PGF, near BC in the south section of the BCF, and area from NB to QC. In the whole process, the release seismic moment is 1.29 × 1021 N·m, and the rupture duration is about 100 s. The slip mainly occurs on the BCF, indicating that the BCF is the main slip plane. The south section of the BCF is dominated by thrust dislocation, and the area around BC has the same slip pattern. There are obvious strike slip dislocations in the area north of Nanba.

Data Availability

The teleseismic records used for the inversion are from the GSN of Incorporated Research Institutions for Seismology (IRIS). The website is http://ds.iris.edu/wilber3/find_event. The strong-motion data are from the “National Strong-Motion Network Center” of Institute of Engineering Mechanics, China Earthquake Administration.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Yin Deyu collected and processed the inversion data, did the work of inversion, and analyzed the inversion result. Meanwhile, YDY wrote the manuscript. Dong Yun and Liu Qifang participated in the work of writing the computer program for the inversion method and analyzed the inversion result. She Yuexin, Wu Jingke, and Zhang Jihua helped to draft the manuscript and participated in the design of this inversion.


This work was supported by the Huaiyin Natural Science Research Plan (HAB202060), Jiangsu Province Industry University Research Cooperation Project in 2020 (BY2020327), the Open Fund for Jiangsu Engineering Laboratory of Structure Assembly Technology on Urban and Rural Residence, Huaiyin Institute of Technology (JSZP201903), and the National Natural Science Foundation of China (51978434, 51804129, 51904112, 51904113, and 11902123) provided the data for this inversion.