Deep coseismic slip of the 2008 Wenchuan earthquake inferred from joint inversion of fault stress changes and GPS surface displacements

Journal of Geodynamics 87 (2015) 1–12

Contents lists available at ScienceDirect

Journal of Geodynamics journal homepage: http://www.elsevier.com/locate/jog

Deep coseismic slip of the 2008 Wenchuan earthquake inferred from joint inversion of fault stress changes and GPS surface displacements Qiang Chen a,∗ , Yinghui Yang a , Rong Luo a , Guoxiang Liu a , Kui Zhang b a b

Department of Remote Sensing and Geoinformation Engineering, Southwest Jiaotong University, Chengdu, China College of Communication Engineering, Chongqing University, Chongqing, China

a r t i c l e

i n f o

Article history: Received 2 July 2014 Received in revised form 4 March 2015 Accepted 4 March 2015 Available online 12 March 2015 Keywords: Deep coseismic slip Joint inversion GPS displacements Stress change Wenchuan earthquake

a b s t r a c t Geodetic data are increasingly being used to infer coseismic slip distribution due to its advantages of wide coverage and high accuracy. However, it is difficult to obtain a comprehensive rupture pattern at depth when a source model is only constrained by geodetic surface deformation. In this study, a joint inversion approach incorporating stress changes and GPS surface displacements is explored and applied to characterize the fault slip of the 2008 Mw 7.9 Wenchuan earthquake, China. The earthquake data for the 20-year period before the main quake, which are collected from the background seismicity catalogues, and one month of aftershock data are statistically analysed to determine the fault stress changes based on the Dieterich model. The coseismic surface deformation measurements from 158 GPS surveying sites are jointly used to constrain the solution. Our preferred rupture model reveals four high-slip concentrations on the Yingxiu–Beichuan fault and one on the subparallel PengGuan fault. The spatial distribution suggests that the coseismic slip occurs not only above the hypocentre but also with a significant thrusting motion, with a mean slip of 8.5 m and a maximum of 9.7 m at a depth of 10–16 km. A significant high-slip concentration is found for the first time in this study. The coseismic faulting extends toward ∼16 km southwest of the Yingxiu–Beichuan fault and has a dextral strike-slip with a mean displacement of 4.8 m at a depth of 7–19 km. The joint inversion model misfits (GPS: 1.7 cm, stress change: 0.02 MPa) exhibit a good compatibility between the two types of datasets. The derived slip model, which has an improved resolution at depth, explains 98% of the coseismic surface displacements and 93% of the fault stress changes. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction The near-fault coseismic surface deformation observed using geodetic techniques, such as the Global Positioning System (GPS) and/or Interferometric Synthetic Aperture Radar (InSAR), is commonly utilized as a significant constraint for fault slip model inversion (Okada, 1985; Massonnet et al., 1993; Feigl et al., 1995; Hao et al., 2009; Shen et al., 2009; Chen et al., 2010, 2012; Wang et al., 2011, 2014; Zhao et al., 2011; Qu et al., 2014). The inversion approach based on geodetic data can be approximately separated into two steps. First, a fault plane is discretized into many rectangular elements or boundary elements, and each of them is considered to be a dislocation source (Okada, 1992; Thomas, 1993; Sun and Okubo, 2002). Second, the inversion Green’s function is calculated based on a unit slip of each fault element. By treating the smoothing constraints as pseudo-observations in the inversion model,

∗ Corresponding author. Tel.: +86 2887601675. E-mail addresses: [email protected], [email protected] (Q. Chen). http://dx.doi.org/10.1016/j.jog.2015.03.001 0264-3707/© 2015 Elsevier Ltd. All rights reserved.

coseismic faulting distributions in an earthquake event can be solved with the use of non-negative least squares algorithms (Lawson and Hanson, 1974; Parker, 1994). If a source model inversion is only constrained by geodetic surface deformation, then the inferred fault slip is typically distributed within a depth of 15 km or less and is tapered below this depth. The model inversion, based on geodetic surface deformation, is actually a mathematical optimization, and the object function is to find the preferred slip model with the minimized misfits for geodetic observations. It can be solved by repeatedly computing a forward model of surface displacements and adjusting the fault slip parameters until the misfits between the modelled and the observed displacements are minimized (Fukahata and Wright, 2008; Weston et al., 2012). This can be mathematically achieved with some non-linear optimization algorithms (e.g., genetic algorithm, simulated annealing algorithm) that are used to vary the source parameters to find the solution with the best fit to the observed geodetic data (Funning et al., 2005). However, surface displacements produced by a unit slip at the shallow crust are more sensitive than those by a unit slip at depth, which means that the inversion process will easily adjust

2

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

the shallow slip to minimize the model residuals in the optimization. This often leads to a deep slip that may not to be effectively recognized in an inversion process due to the limitation of inverting geodetic data-only without using additional constraints. Existing investigations (Somerville and Pitarka, 2006; Song and Somerville, 2010) have found that the inferred slip depth from geodetic data-only is shallower than that from seismic data. The median source depth difference between geodetic data and seismic data can reach −5.7 km and −2.7 km for the strike-slip and thrust faulting mechanisms, respectively (Weston et al., 2012). Some studies (Somerville and Pitarka, 2006; Song et al., 2009) have also suggested that, once a fault rupture begins, it can propagate both above and below the seismogenic zone. A small seismic event may rarely break the surface because its slip is limited to the seismogenic depth. However, in contrast, a high magnitude earthquake can propagate both upward from the hypocentre to the surface and downward below the seismogenic zone, which may cause a large slip at depth (Song et al., 2009; Song and Somerville, 2010). In these cases, geodetic data-only constrained source models appear to be inconsistent with the real coseismic faulting pattern with both shallow and deep slips (Ji and Hayes, 2008; Shen et al., 2009; Zhao et al., 2010; Wang et al., 2011). King and Wesnousky (2007) have reported that, theoretically, each tapered shape slip model of depth distribution for a strike-slip earthquake can correspond to a box model, and these slip models share an almost identical surface displacement field. Song et al. (2009) and Song and Somerville (2009, 2010) performed geodetic slip inversion for a strike-slip event of the 1999 Kocaeli earthquake, Turkey, to test the resolving power of geodetic data, especially as a function of depth. The test results show that the resolution at depth does not arise significantly when there is an increase in the number of singular values of the slip inversion model. In addition, the spatial slip resolution decays rapidly with increasing depth, owing to the distance squared decay of the surface static displacement field from the source, which implies that the slip will be limited to the upper few kilometres. Therefore, a fault slip in the deep crust probably exists but may not be precisely detected by the geodetic data-only inversion. Likewise, this phenomenon exists in the case of the 2008 Wenchuan earthquake, which had a magnitude of Mw 7.9 and occurred on the western rim of the Chengdu basin in Sichuan Province, China on 12 May 2008 (Burchfiel et al., 2008; Huang and Li, 2009; Shen et al., 2009; Tong et al., 2010). The earthquake is the country’s largest seismic event over the past 50 years and caused more than 69,200 fatalities, 18,195 persons missing and 374,216 injured (Cui et al., 2009; Tang et al., 2011; Chen et al., 2014). The seismological record data indicates that the rupture initiated in the southwestern Longmen Shan fault zone, which is at the eastern margin of the Tibetan Plateau, and propagated unilaterally toward the northeast along a northwest dipping fault (Xu et al., 2009). The earthquake ruptured both the Yingxiu–Beichuan fault at a length of approximately 320 km and the PengGuan fault at a length of approximately 80 km (see Fig. 1), which are linked by a short northwest-striking rupture zone at the southern end of the PengGuan fault through a lateral ramp, known as the Xiaoyudong rupture zone (Liu et al., 2009; Chen et al., 2014). The geometry of the fault plane is northwest dipping. However, the dip angle shows strong variations; it is almost vertical (dip angle of ∼75◦ ) near the ground surface, then tends to decrease with the depth, finally converging into the same ramp system with another subparallel PengGuan fault in the mid-crust (Li et al., 2009; Wang et al., 2011). The faulting model inferred from GPS and/or InSAR deformation inversion reveals three high-slip concentrations, which have, respectively, 5 m of thrust slip from Yingxiu to Xiaoyudong at a depth of 0–10 km, 5.2 m and 4.8 m of thrust and dextral slip

near Beichuan at a depth of 0–7 km, and 3.0 m of strike-slip near Qingchuan at 0–5 km (Hao et al., 2009; Shen et al., 2009; Tong et al., 2010). However, coseismic faulting below 15 km is not recognized in the previous studies using only geodetic data (Hao et al., 2009; Shen et al., 2009; Diao et al., 2010; Tong et al., 2010; Wang et al., 2011, Yang et al., 2014), while the inversions based on seismic data indicate deeper slip (Ji and Hayes, 2008; Zhao et al., 2010). Therefore, it is instructive to incorporate additional data sets into coseismic faulting inversion to obtain a comprehensive slip distribution, especially at depth. In this study, fault stress changes caused by the main shock, which are highly sensitive to fault slip at depth, are used as the complementary data for the slip inversion. The stress change values can impose direct constraints on the fault plane, implying that they have the potential to provide an improved resolving power to map slips at depth. By using the recorded earthquake occurrence rate data over the past 20 years before the main quake and one month of aftershock data, we estimated the seismic fault stress changes based on Dieterich’s model. The estimated fault stress changes and observed GPS surface displacements are jointly used to constrain the source model, thereby inferring the comprehensive faulting pattern for the event. 2. Estimate of fault stress changes 2.1. Dieterich’s aftershock model High-magnitude earthquakes can cause changes in fault stress distribution and a series of aftershocks (Hanks, 1977; Zhang et al., 2003; Wu et al., 2006). The short-term aftershock occurrence rate is much higher than the pre-seismic long-term average earthquake rates. Some studies (Dieterich, 1994; Dieterich et al. 2000) have indicated that there exists a widespread correlation between stress changes and seismic activity changes before and after a main shock. Dieterich et al. suggested that data statistics analysis of earthquake occurrence rates can provide an available remotesensing metre of fault stress (Dieterich, 1994; Dieterich et al. 2000). Stress changes can be quantitatively estimated from the earthquake occurrence rates recorded in seismicity catalogues. Hence, in this study, Dieterich’s aftershocks model (Dieterich et al., 2000), which employs a constitutive friction coefficient that depends on the logarithm of the sliding speed, and the relation between coseismic stress changes and instantaneous aftershock rate changes are adopted to estimate stress changes. It can be expressed as follows. ı = Aln(R+ /R− )

(1)

where, ı denotes the coseismic stress change, A is a constitutive parameter that weights the dependence of the friction coefficient on the logarithm of slip rate,  indicates the effective normal stress that will be set equal to the lithostatic pressure (Dieterich, 1994; Dieterich et al. 2000; Ziv, 2012), and R+ and R− are the short-term aftershocks rate and pre-seismic long-term earthquake occurrence rate, respectively. The sign of division means that it is resolvable only for fault areas that were seismically active both before and after the main shock. These fault areas can remain as a similar spatial distribution of preseismic small earthquakes and aftershock activities. The Wenchuan earthquake, which occurred in the Longmen Shan fault zone with active seismic events, provides a good example for the application of stress change estimation. Dieterich’s aftershock model has been used to assess probabilistic seismic hazards (Toda et al., 1998, 2002), estimate the tectonic stressing rate from the spatiotemporal evolution of Loma Prieta aftershocks (Gross and Burgmann, 1998) and infer stress change due to dike intrusion in the Kilauea volcano in Hawaii (Dieterich et al., 2000). These studies have demonstrated that this model has

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

3

Fig. 1. Main ruptured faults and distribution of seismic activities before and after the main shock of the Wenchuan earthquake over the Longmen Shan fault zone, mapped onto a shaded relief of the SRTM digital terrain model. The blue dots denote the background earthquake events from May 1988 to May 2008, while the red dots indicate one month aftershocks (Huang et al., 2008; Zhu et al., 2008; Zhao et al., 2010). The black lines indicate the surface projection of the Yingxiu–Beichuan and PengGuan faults (Li et al., 2009; Liu et al., 2009). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

a good potential to infer fault stress changes associated with strong earthquakes (Dieterich and Kilgore, 1996; Toda et al., 2002; Ziv, 2012). In this study, fault stress changes can be considered as a significant complement (i.e., pseudo-data) to constrain the source model. Noting this, we incorporated the fault stress change data estimated from the Dieterich aftershock model into faulting inversion to reveal a comprehensive slip pattern associated with the Wenchuan earthquake. The fault plane was first discretized into many elements before projecting the earthquake events onto it. The investigations by Zhang et al. (2011) show that the slip distributions, respectively derived from the different discretization sizes of 4 km, 2 km and 1 km, agree well with each other for the Wenchuan earthquake (Zhang et al., 2011). Therefore, considering the huge fault length of almost 320 km associated with the main event and the computation load, we divided the Yingxiu–Beichuan fault into five segments along the strike and discretized the fault plane to a grid of 2840 elements (sub-fault blocks); the dimension of each element is 4 km by 2 km along the strike and dip directions, respectively. We collected the earthquake occurrence rate data for the 20 year period before the main quake from the background seismicity catalogues. As for the calculation of R+ related to the aftershock data, once an aftershock sequence begins, each aftershock is considered to be acting locally as a small main shock by changing the stress field in its surrounding areas (Dieterich, 1994; Ziv, 2012). The aftershocks that happened long after the event may not be triggered by the main shock but, rather, earlier nearby aftershocks. Therefore, the collection time of the aftershock data was first determined to calculate the ratio of R+ /R− and stress changes. To test the stability and reliability of the ratio, the aftershocks from 5 days to 200 days after the main event are used to calculate the earthquake occurrence rate ratio R+ /R− . The resulting time series

ratio of R+ /R− is shown in Fig. 2. It can be seen that the ratio has a trend of decreasing with time, varying from 5 to 200 days after the main quake. The spatial distribution pattern and magnitude of R+ /R− show good agreement from Fig. 2(c) to (f), i.e., from 15 days to 30 days after the main event, which means that the aftershock occurrence rates during the phase reach a relatively stable status. After the phase, the ratio continues to decrease until the almost preseismic activity level. Hence, we choose one month of aftershock data to calculate the ratio of earthquake occurrence rates. It should be noted that the determination of stress changes is a logarithmic relationship with respect to the ratio, thus producing the relatively stable stress change values. The ratio R+ /R− can be calculated after the projection of the earthquake occurrence rate data onto the fault plane. The projection results show that the fault patches, which were seismically active both before and after the main shock, are almost located on the causal fault, i.e., the Yingxiu–Beichuan fault. Fig. 3(a) shows the resulting earthquake rate change (R+ /R− ), where the green dots indicate the location of fault patches and the contours indicate the value of R+ /R− , which details the distribution of seismic activity change on the fault patches (Huang et al., 2008; Zhu et al., 2008). Fig. 3(b) shows the bicubic interpolation result based on the midpoints of the fault grid. 2.2. Determination of friction parameter To obtain a reasonable and reliable friction parameter, namely, the value of sign A in Eq. (1), we performed the inversion testing with different values of the friction parameter within the interval (0.001, 0.016) of typical laboratory values (Dieterich, 1994; Ziv, 2012). Fig. 4 shows the stress change misfit percentage between the modelled and observed, as well as the GPS RMS versus the

4

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

Fig. 2. Ratios of earthquake occurrence rates from different aftershocks collection times. The green dots indicate the double-differenced location of earthquake events on fault patches, which were seismically active both before and after the main shock. The contours denote the value of R+ /R− . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

change of friction parameter A. It should be noted that the existing studies (Fukuda et al., 2009, Foster et al., 2013) have proposed that the friction parameter can be determined based on the slow slip and afterslip derived from a long time series of GPS data. Unfortunately, the long time series of GPS data are not available over the Longmen Shan fault zone because of the rough terrain and geographical impediments in the study area. Therefore, the optimal forward testing with different values was here used to determine the best friction parameter in the study. It can be seen from Fig. 4 that the stress change misfit decreases rapidly with the increase of the friction parameter within the interval [0.001, 0.005], then gradually reaches a stable level when the parameter is within the interval [0.005, 0.016]. The high friction parameter does not improve the misfit much. In addition, the influence of the friction parameters on the GPS RMS is very slight, and the largest discrepancy of the GPS RMS is less than 0.2 mm. Therefore, we pick the preferred friction parameter of 0.005. Then, the stress changes can be calculated based on Eq. (1) and the earthquake rate change data in Fig. 3(b) after determining the friction parameter. Fig. 3(c) shows the resulting stress changes on the fault plane.

The stress changes estimated from the seismic activity change can be considered as the pseudo-observations in the following source model inversion. Green’s function of stress change is constructed based on the Okada elastic dislocation model (Okada, 1992) with Poisson’s ratio and Young’s modulus of 0.25 and 70 GPa (Perrier and Ruegg, 1973; Dziewonski and Anderson, 1981; Ziv, 2012), respectively. The relationship between fault slip and stress change can be expressed as follows. s Kij uj = s ıi

(2)

where,  s indicates the model weight, Kij is Okada’s elastic kernel for dislocation, related to a unit slip on dislocation j with a stress change on site i, while i, j denote the indexes of the seismically active fault elements for which the R+ /R− is resolvable, uj is a permanent slip on the j fault patch, ıi is the stress changes calculated from Eq. (1). 3. Joint inversion with stress changes and GPS data The coseismic surface displacements associated with the Wenchuan earthquake were monitored by the Crustal Motion

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

5

Fig. 3. Earthquake occurrence rate ratio and estimated stress change along the Yingxiu–Beichuan fault. (a) The ratio between one month of aftershocks and 20 years of background seismicity. The red star indicates the hypocentre. (b) Earthquake rate change diagram from the bicubic interpolation on the midpoints of the fault grid. (c) Estimated stress change map on the fault plane, calculated from seismic activity change data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4. Stress changes misfits and GPS residuals versus the given friction parameters.

Observation Network of China Project, which consists of 158 GPS stations over the Longmen Shan fault zone near the eastern edge of the Tibetan Plateau. Table 1 shows the collection time of GPS data after the main quake. It can be seen that 99.4% of GPS data were collected within 30 days after the quake, which can match well with the collection time of the used aftershocks data. These GPS observations were acquired with a data sampling interval of 30 s and 2∼3 observation sessions. Each session is composed of almost one day of data acquisition. The continuous GPS data were processed to derive daily solutions based on the double difference model and precision ephemeris distributed by the IGS. As a result, the coseismic displacements at the 158 GPS stations were determined. Unfortunately, due to the uncertainties of antenna setups and/or low signal-to-noise ratio at some campaign-surveyed GPS stations, the vertical component is only available at 44 stations with a relatively low precision (Working Group of the Crustal Motion Observation Network of China Project, 2008). Therefore, only the

Table 1 Collection time of GPS data. Days after main quake

5

10

15

20

25

30

40

GPS number Percentage (%)

45 28.5

72 45.6

93 58.9

133 84.2

147 93.0

157 99.4

158 100.0

6

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

Fig. 5. Coseismic displacements at the GPS observation stations. Unit lengths, indicated by brown and blue vectors, stand for 40 cm and 10 cm horizontal displacements, respectively. Red dots denote aftershocks greater than Ms 4.0. Black lines depict the traces of the Yingxiu–Beichuan and PengGuan faults. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

horizontal components of coseismic displacements were used in the inversion process. Fig. 5 shows clearly the displacement vectors derived from the GPS observations. The significant surface displacements, with the values of 1.5∼3.0 m, can be seen over the Dujiangyan, Beichuan and Qingchuan areas, which are located near the fault junctions and subject to the most severe structural damage and life loss (Cui et al., 2009; Tang et al., 2011). For the purpose of comparison, we first conducted a GPS-only fault slip inversion. Eq. (3) shows the corresponding coseismic slip inversion model.

⎡

GN1

GN2

⎢ ⎢ GE1 GE2 ⎢ ⎢ L 0 ⎣ 1 0

L2

⎤

⎡ ⎤ dN  ⎥ ⎢ ⎥ ⎥ U1slip ⎢ dE ⎥ ⎥ = ⎢ ⎥ ⎥ U ⎣0 ⎦ ⎦ 2slip

(3)

0

where, GN and GE indicate Green’s functions with respect to the north and east displacement components, respectively, Li denotes the Laplace smoothing matrix,  is the smoothing weight, U1slip and U2slip denote the strike-slip and thrust-slip components to be solved, and dN and dE are the observed displacements in the north and east directions, respectively. Fig. 6 shows the GPS-only inversion result of coseismic slip on the Yingxiu–Beichuan fault and the PengGuan fault. The derived fault slip model suggests three high-slip concentrations on the Yingxiu–Beichuan fault. First, the slip is mostly thrust with an average of 5 m and a maximum of 9.0 m in the southwestern fault segment and occurred mostly in the upper crust within a depth of 10 km. Then, the fault rupture progressively transformed to thrust and dextral strike slip with an average of 5.2 and 4.8 m in the middle fault segment at a depth of 0–7 km. Finally, the rupture pattern becomes the dominant dextral strike slip of 3.0 m in the northeastern fault segment at a depth of 0–5 km.

The spatial distribution of fault slip from the GPS-only inversion agrees well with the geodetic data obtained in previous studies, including GPS and/or InSAR observations (Hao et al., 2009; Shen et al., 2009; Diao et al., 2010; Tong et al., 2010; Wang et al., 2011, Yang et al., 2014). The coseismic slips on fault planes are with varying characteristics along the strike direction. The significant high-slip concentrations occurred on the surface near the Dujiangyan, Beichuan and Qingchuan areas located in the southwest, middle and northeast segments of the YingXiu–Beichuan fault, respectively. However, the derived coseismic slip is mostly concentrated within a shallow depth of 10 km or less, and the slip magnitude decreases rapidly with the depth. Moreover, the coseismic slip is hardly found below the hypocentre and even around the source, which may be inconsistent with the real faulting pattern. To overcome the problem of weak sensitivity to coseismic slips at depth inverting only geodetic data, we incorporated the estimated fault stress changes and GPS data to jointly constrain the source model. The stress changes data are considered as a type of pseudo-observations in the joint inversion, while the GPS surface displacements are considered as the real observations. Green’s functions of stress change and surface deformation can be calculated from Eqs. (2) and (3), respectively. Eq. (4) shows the joint inversion model constrained by the GPS displacements and stress changes.

⎡ ˛G

N1

˛GN2

⎢ ˛G ⎢ E1 ˛GE2 ⎢ ⎢ ˇKS1 ˇKS2 ⎢ ⎢ ⎣ L1 0 0

L2

⎤

⎡

˛dN

⎤

⎥ ⎢ ˛d ⎥ ⎥ ⎢ E ⎥ ⎥ U1slip ⎢ ⎥ ⎥ = ⎢ ˇı ⎥ ⎥ U ⎢ ⎥ ⎥ 2slip ⎣0 ⎦ ⎦

(4)

0

where, KS1 and ı indicate Green’s function and stress changes data, respectively, ˛ and ˇ are the relative model weights, and the other parameters have been defined previously.

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

7

Fig. 6. Coseismic slip derived from GPS-only inversion with planar view (top panel) and three-dimensional view (bottom panel) on the Yingxiu–Beichuan and PengGuan faults.

The model weights regarding the two types of data sets have to be considered in the joint inversion. It is believed that the weights will be adjusted to achieve a slip distribution, resulting in a similar sensitivity and an equal normalized misfit for the two data sets (Ziv, 2012). The normalized misfit of the GPS observations and stress change data are defined as the following Eq. (5).



 n 

Wi diobs − diest    ε= n × Wi dobs  i=1

(5)

i

where, W is the model weight factor, d indicates either GPS surface displacements or fault stress changes, the superscripts “obs” and “est” signify “observed” and “estimated”, respectively, n denotes the number of observations, and the sign “|| ||” indicates L2 norm operator (Ziv, 2012). Given an initial smoothing factor , the trial-and-error algorithm is used to determine the equal normalized misfit and the corresponding model weights. The smoothing factor  can be redefined after determining the model weights. The model misfit of the two observations is statistically analysed from the continuous smoothing factors. Fig. 7 shows the trade-off relationship between the misfits of two data sets and solution roughness. We pick a roughness of 3.0 cm/km for the solution, as lower roughness results in a significantly worse misfit but higher roughness does not improve the misfit much. Then, we fix the new smoothing factor and reuse the trial-and-error algorithm to obtain new model weights. Finally, the stable smoothing factor of 13.0 and model weights of 1.2 and 1.0 for GPS and stress change data, respectively, can be obtained through repeating the above process.

Fig. 7. Trade-off relationship with different misfits of two observations as a function of average solution roughness.

4. Results and discussions 4.1. Slip distribution from joint inversion The joint inversion through Eq. (4) can be solved in a nonnegative least squares sense (Lawson and Hanson, 1974; Parker, 1994). Fig. 8 shows the derived spatial distribution of fault slip related to the 2008 Wenchuan earthquake. Four high-slip concentration areas are found along the strike of the Yingxiu–Beichuan fault and one can be seen along the subparallel PengGuan fault. Overall, the faulting pattern retains a similar spatial distribution to that from geodetic

8

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

Fig. 8. Coseismic slip from joint inversion with planar view (top panel) and three-dimensional view (bottom panel) on the Yingxiu–Beichuan and PengGuan faults.

data-only inversion within the shallow crustal range of 10 km. However, the slip distribution below a depth of 10 km can be further mapped from the joint inversion. The faulting pattern varies along the strike direction of the Yingxiu–Beichuan fault. First, the fault rupture occurred with mostly thrusting in the southwest segment of the Yingxiu–Beichuan fault from the ground surface to a depth of 17 km, then it progressively transformed to a mixed mode of dextral slip and thrust motion in the middle segment, finally becoming a dominant right-lateral slip in the northeast segment. The maximum and mean values of the coseismic slip reach 9.7 m and 8.5 m near the initial rupture area, respectively. However, compared to the GPS data-only inversion results (See Fig. 6), the coseismic slip below the crustal depth of 10 km is obvious from the joint inversion. Two high-slip concentration areas, marked by “Patch-1” and “Patch-2” in Fig. 8(a), respectively, and located in the northeast 20 km of the hypocentre, have to be noted. The slip concentration “Patch-1” is essentially consistent with the previous results at a depth of 0–5 km (Shen et al., 2009; Tong et al., 2010; Wang et al., 2011; Yang et al., 2014). However, “Patch-2” is located in a depth ranging from 10 to 16 km, which is a significant slip area inaccessible to inverting GPS data only. It seems that the single slip event is divided into two high-slip concentrations, i.e., “Patch-1” and “Patch-2”. The slip “Patch-2” covers more than 70 km2 with a maximum and average slip of 9.7 m and 8.5 m, respectively. In addition, the joint inversion model contains a slip distribution labelled “Patch-6”, which is at a depth of 13.5 km in the middle fault segment near the Beichuan area, that is deeper than the GPS-only one. Another new high-slip concentration, which extended 16 km southwest to the hypocentre along the Yingxiu–Beichuan fault and occurred mostly at the depth ranging from 7 to 19 km, is recognized, with the right-lateral strike slip mode marked by “Patch-3” in Fig. 8(a). The slip area reaches ∼320 km2 and has an average slip of 4.8 m. However, this fault segment produces a small coseismic

surface displacement of less than 1.0 m and does not form a significant surface break trace, which agrees well with the field geological investigations (Liu et al., 2009). In addition, the slip area of “Patch3” extends to a decollement at a depth of more than 20 km, with two relatively slight slips corresponding to “Patch-4” and “Patch5” (See Fig. 8(c)). The slip area “Patch-4” is mainly characterized by thrust motion with an average slip of 3.8 m, while “Patch-5” shows a dominant right-strike slip with an average value of 3.1 m. To validate the reliability of slip inversion using geodetic data and stress changes, the source model derived from the seismic data inversion by Zhao et al. (2010) is introduced here as the reference, which is shown in Fig. 9. It reveals that the Wenchuan earthquake is composed of five distinguishable sub-events during the whole rupture process, which are marked as the five black boxes and named “a–e” in Fig. 9. The five rupture sub-events include one right-strike

Fig. 9. Fault slip model derived from seismic data for Wenchuan earthquake by Zhao et al. (2010). The red star indicates the hypocentre. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

9

Fig. 10. Discrepancies of strike slip and thrust motion between joint inversion and GPS-only inversion. (a) Strike slip component. (b) Thrust motion component.

slip near the source, two thrust motions over the Dujiangyan area located in the southwest fault segment, and two ruptures composed of right-strike slip and thrust motion over the Beichuan and Qingchuan area. Analysing the source model derived from the joint inversion (See Fig. 8), the inferred slip “Patch-3” in the study agrees well with the right-strike slip event near the source labelled as box “a”. “Patch-1” and “Patch-2” correspond to the two thrust motions (i.e., boxes “b” and “c”) in the southwest fault segment. “Patch-6” and “Patch-7” agree with the right-strike slip and thrust motion in the Beichuan and Qingchuan area (i.e., boxes “d” and “e”). However, the GPS-only inversion slip distribution lacks the deep slip concentrations corresponding to boxes “a” and “c” in Fig. 9. Overall, the slip distribution from the joint inversion shows better agreements with the seismic result than that from the GPS-only inversion. In particular, the two rupture sub-events, marked as “Patch-2” and “Patch-3” near the source in the southwest fault segment, are clearly recognized in the joint inversion. 4.2. Slip difference between GPS-only and joint inversions To quantify the slip difference between the GPS-only and joint inversions, a differential operation was conducted on the two slip models. Fig. 10(a) and (b) show the discrepancies along the strike slip and thrust motion directions between the two models, respectively. It can be seen from Fig. 10(a) that the right-lateral strike slip distribution in the upper crust of 10 km is generally consistent between the two models, with a mean absolute deviation (MAD) of less than 0.39 m. The main differences are observed in the areas of “Patch-3” and “Patch-5” at the depth range of 10–20 km. In

particular, the area “Patch-3” is a newly found high-slip concentration in the joint inversion, and its MAD exceeds 3.4 m. The slip difference on the PengGuan fault is generally slight, except for one asperity, marked by “Patch-8”, with a MAD of 0.8 m. It can be seen from Fig. 10(b) that the slip difference in the thrust motion direction is mainly located in the depth ranging from 10 to 20 km along the Yingxiu–Beichuan fault, concentrating particularly on the areas marked by “Patch-2” and “Patch-4”. The area “Patch-2” is around the source and has a MAD exceeding 3.7 m, while “Patch4” is at the decollement and has a MAD of 1.2 m. There is a slight difference above the upper crust of 10 km with a MAD of 0.67 m. It can be found from Fig. 10 that the slip magnitude from the GPSonly inversion is generally larger than that from the joint inversion in the upper crust at a depth of 10 km, but opposite in the depth range of 10–20 km. This balance and readjustment of slip distribution probably results from the poor resolution of the GPS-only inversion at depth, which generates a larger shallow slip for the minimization of model misfit through the removal of deep slip. Sharing an almost identical surface displacement field, the shallow slip magnitude from the GPS-only inversion is artificially enlarged with respect to that from the joint inversion. Conversely, the joint inversion model, which has an improved resolution for deep slip, contains a relatively smaller slip in the shallow crust. The difference between the two slip patterns probably resulted from whether the additional constraints (e.g., stress changes) on the model solution in the mathematical optimization of the inversion process were used. This may be the main limitation of the geodetic data-only inversion without additional constraints used, which made the deep slip difficult to derive in the inversion process.

10

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

Fig. 11. Fault stress change maps from the model prediction (a) and residuals with respect to estimates from Dieterich’s model and seismic activity change (b). The red star indicates the hypocentre. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

4.3. Model misfits of joint inversion To assess the compatibility between GPS observables and stress changes in the model inversion, we calculated the residuals of the two types of constraint data. Fig. 11(a) shows the predicted fault stress changes, which are calculated from the coseismic slip models of joint inversion. Fig. 11(b) shows the stress change difference between the observed and predicted data. The corresponding MAD is 0.02 MPa, which is less than 10% of the mean absolute value of observations (equal to 0.24 MPa). The statistics suggest that the slip model from the joint inversion explains 93% of the stress changes.

The residuals between the GPS displacement observations and the model prediction are shown in Fig. 12. The mean GPS misfits from the joint inversion reach 1.6 cm and 1.8 cm for the north and east components, respectively, which is significantly less than the results from the inversion using only GPS and/or InSAR data (Hao et al., 2009; Shen et al., 2009; Diao et al., 2010; Tong et al., 2010; Wang et al., 2011). Most model misfits are less than 1.0 cm, and more than 98% of the GPS displacements can be well explained in the joint inversion slip model. In addition, few residuals greater than 2.0 cm mainly concentrate over the elliptical area shown in Fig. 12. It is believed that the residuals probably result from

Fig. 12. Spatial distribution of GPS residuals from the source model. Mean misfits of GPS data from the joint inversion are 1.6 cm and 1.8 cm for the north and east components, respectively. Unit lengths, indicated by the brown and blue vectors, stand for 5 cm and 1 cm horizontal displacements, respectively. The red lines indicate the surface projection of the Yingxiu–Beichuan and PengGuan faults. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

a small unmodelled rupture occurring on what is referred to as the Xiaoyudong fault, which is approximately perpendicular to the causal fault. The total geodetic moment estimated from our preferred slip model is 9.4 × 1020 Nm (equivalent to a moment magnitude of 7.96), which is slightly greater than the 9.2 × 1020 Nm (equivalent to a moment magnitude of 7.91) derived from the GPSonly inversion slip model. The moment magnitude of our preferred slip model is close to the existing results from the geodetic inversion using GPS and/or InSAR data but also larger than the seismic moment magnitude of 7.8 from Zhao et al. (2010), which may result from the additional surface displacements caused by the few aftershocks in the period of used data (Shen et al., 2009; Diao et al., 2010; Tong et al., 2010; Wang et al., 2011; Zhao et al., 2010).

5. Conclusions The geodetic surface displacements have the limited constraint capability of resolving deep coseismic slip because the spatial resolution decays rapidly with depth. In this study, we utilize an integrated inversion approach based on GPS displacement observations and stress changes to jointly constrain the source model solution. The changes in seismic activities before and after the main quake are used to estimate fault stress changes based on Dieterich’s aftershock model. The model weights between the two data sets are determined through an optimization process of fault roughness. Green’s functions of stress changes and surface deformation are constructed and the slip model can be solved in a nonnegative least squares sense. The joint inversion approach is applied to infer the source model of the 2008 Mw 7.9 Wenchuan earthquake. The historical earthquake occurrence rates during the 20 years before the main shock and one month of aftershock data were statistically analysed to estimate fault stress changes. The coseismic surface displacements observed from the 158 GPS surveying sites over the Longmen Shan fault zone were used as the constraints in the joint inversion. Our preferred slip model suggests that the coseismic slip occurs not only within the upper crust at a depth of 10 km above the hypocentre but also with a significant thrusting motion of an average offset 8.5 m at a depth of 10–16 km, which is inaccessible to the geodetic data-only inversion. A significant high-slip concentration is newly recognized in this study. It reveals that the coseismic faulting extends toward ∼16 km southwest of the Yingxiu–Beichuan fault with a dextral strike-slip of mean 4.8 m at a depth of 7∼19 km. The model misfits of joint inversion show good compatibility between the two types of data sets, with GPS residuals of 1.7 cm and stress changes residuals of 0.02 MPa. The derived slip model explains well 98% of the coseismic surface displacements and 93% of the fault stress changes. This study also indicates that the data used to constrain the coseismic faulting model could be enriched both pre-seismically and post-seismically (Ziv, 2012). The joint inversion, incorporating geodetic data and seismic activity data, not only has the advantage of improved resolution at depth but also provides mechanical insight on the fault motion properties.

Acknowledgements This work is supported by the National Natural Science Foundation of China (nos. 41472255, 51178404). We are grateful to the Crustal Motion Observation Network of China Project for providing the GPS data and the China Earthquake Networks Center for providing the seismicity catalogue. Most of the figures were generated using the Generic Mapping Tool (GMT) software (Wessel and Smith, 1998).

11

References Burchfiel, B.C., Royden, L.H., Hilst, R.D., Hager, B.H., Chen, Z., King, R.W., Li, C., Lü, J., Yao, H., Kirby, E., 2008. A geological and geophysical context for the Wenchuan earthquake of 12 May 2008, Sichuan, People’s Republic of China. Geol. Soc. Am. Today 18 (7), 4–11, 10.1130/GSATG18A.1. Cui, P., Zhu, Y.Y., Han, Y.S., 2009. The 12 May Wenchuan earthquake-induced landslide lakes: distribution and preliminary risk evaluation. Landslides 6, 209–223, 10.1007/s10346-009-0160-9. Chen, Q., Liu, G.X., Ding, X.L., Hu, J.C., Yuan, L.G., Zhong, P., Omura, M., 2010. Tight integration of GPS observations and persistent scatterer InSAR for detecting vertical ground motion in Hong Kong. Int. J. Appl. Earth Obs. Geoinf. 12, 477– 486. Chen, Q., Liu, G.X., Hu, J.C., Ding, X.L., Yang, Y.H., 2012. Mapping ground 3-D displacement with GPS and PS-InSAR networking in the Pingtung area, southwestern Taiwan. Chin. J. Geophys. 55 (10), 3248–3258 (in Chinese). Chen, Q., Cheng, H.Q., Yang, Y.H., Liu, G.X., Liu, L.Y., 2014. Quantification of mass wasting volume associated with the giant landslide Daguangbao induced by the 2008 Wenchuan earthquake from persistent scatterer InSAR. Remote Sens. Environ. 152, 125–135. Diao, F.Q., Xiong, X., Wang, R.J., Zheng, Y., Houtse, H., 2010. Slip model of the 2008 Mw7.9 Wenchuan (China) earthquake derived from co-seismic GPS data. Earth Planets Space 62, 869–874, 10.5047/eps.2009.05.003. Dieterich, J.H., 1994. A constitutive law for rate of earthquake production and its application to earthquake clustering. J. Geophys. Res. 99, 2601–2618, 10.1029/ 93JB02581. Dieterich, J.H., Kilgore, B., 1996. Implications of fault constitutive properties for earthquake prediction. Proc. Natl. Acad. Sci. USA 93, 3787–3794, 10.1073/pnas. 93.9.3787. Dieterich, J.H., Cayol, V., Okubo, P., 2000. The use of earthquake rate changes as a stress meter at Kilauea volcano. Nature 408, 457–460, 10.1038/35044054. Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference earth model. Phys. Earth Planet. Inter. 25, 297–356, 10.1016/0031-9201(81)90046-7. Feigl, K.L., Sergent, A., Jacq, D., 1995. Estimation of an earthquake focal mechanism from a satellite radar interferogram: application to the 4 December 1992 landers aftershock. Geophys. Res. Lett. 22 (9), 1037–1048, 10.1029/94GL03212. Foster, J.H., Lowry, A.R., Brooks, B.A., 2013. Fault frictional parameters and material properties revealed by slow slip events at Kilauea volcano, Hawaii. Geophys. Res. Lett. 40 (23), 6059–6063, 10.1002/2013GL058234. Fukuda, J., Johnson, K.M., Larson, K.M., et al., 2009. Fault friction parameters inferred from the early stages of afterslip following the 2003 Tokachi-oki earthquake. J. Geophys. Res. 114, B04412, 10.1029/2008JB006166. Fukahata, Y., Wright, T.J., 2008. A non-linear geodetic data inversion using ABIC for slip distribution on a fault with an unknown dip angle. Geophys. J. Int. 173 (2), 353–364, 10.1111/j.1365-246X.2007.03713.x. Funning, G.J., Parsons, B., Wright, T.J., Jackson, J.A., Fielding, E.J., 2005. Surface displacements and source parameters of the 2003 Bam (Iran) earthquake from Envisat advanced synthetic aperture radar imagery. J. Geophys. Res.: Solid Earth 110, B09406, 10.1029/2004JB003338. Gross, S., Burgmann, R., 1998. Rate and state of background stress estimated from the aftershocks of the 1989 Loma Prieta, California, earthquake. J. Geophys. Res. 103, 4915–4928, 10.1029/97JB03010. Hanks, T.H., 1977. Earthquake stress-drops, ambient tectonic stresses, and the stresses that drive plates. Pure Appl. Geophys. 115, 441–558, 10.1007/ BF01637120. Hao, K.X., Si, H.J., Fujiwara, H., Ozawa, T., 2009. Coseismic surface-ruptures and crustal deformations of the 2008 Wenchuan earthquake Mw7.9, China. Geophys. Res. Lett. 36, L11303, 10.1029/2009GL037971. Huang, Y., Wu, J.P., Zhang, T.Z., Zhang, D.N., 2008. Relocation of the M8.0 Wenchuan earthquake and its aftershock sequence. Sci. China Ser. D 51 (12), 1703–1711, 10.1007/s11430-008-0135-z. Huang, R.Q., Li, W.L., 2009. Analysis of the geo-hazards triggered by the 12 May 2008 Wenchuan earthquake, China. Bull. Eng. Geol. Environ. 68, 363–371, 10.1007/s10064-009-0207-0. Ji, C., Hayes, G., 2008. Preliminary result of the May 12, 2008 Mw7.9 eastern Sichuan, China earthquake, U.S. Geol. Surv., Reston, Va. http://earthquake.usgs.gov/ eqcenter/eqinthenews/2008/us2008ryan/finite fault.php. King, G.C., Wesnousky, S.G., 2007. Scaling of earthquake fault parameters. Bull. Seismol. Soc. Am. 97, 1833–1840. Lawson, C.L., Hanson, R.J., 1974. Solving Least Squares Problems. Prentice Hall, Englewood Cliffs, N. J. Li, Y., Huang, R.Q., Densmore, A.L., Zhou Rongjun, Shuyou, 2009. Basic features and research progresses of Wenchuan Ms 8. 0 earthquake. J. Sichuan Univ. 41 (3), 07–25 (in Chinese). Liu, Z.J., et al., 2009. Co-seismic ruptures of the 12 May 2008, Ms 8.0 Wenchuan earthquake, Sichuan: east–west crustal shortening on oblique, parallel thrusts along the eastern edge of Tibet. Earth Planet Sci. Lett. 286 (3–4), 355–370, 10.1016/j.epsl.2009.07.017. Massonnet, D., Rossi, M., Carmona, C., et al., 1993. The displacement field of the landers earthquake mapped by raddar interferometry. Nature 364, 138–142, 10.1038/364138a0. Okada, Y., 1985. Surface deformation to shear and tensile fault in a half-space. Bull. Seismol. Soc. Am. 75 (4), 1135–1154. Okada, Y., 1992. Internal deformation due to shear and tensile faults in a half-space. Bull. Seismol. Soc.Am. 82, 1018–1040. Parker, R.L., 1994. Geophysical Inverse Theory. Princeton Univ Press, Princeton, N. J.

12

Q. Chen et al. / Journal of Geodynamics 87 (2015) 1–12

Perrier, G., Ruegg, J.C., 1973. Structure profonde du Massif Central Francais. Ann. Geophys. 29 (4), 435–502. Qu, W., Lu, Z., Zhang, Q., Li, Z.H., Peng, J.B., Wang, Q.L., Drummond, J., Zhang, M., 2014. Kinematic model of crustal deformation of Fenwei basin, China based on GPS observations. J. Geodyn. 75, 1–8. Sun, W.K., Okubo, S., 2002. Effect of earth’s spherical curvature and radial heterogeneity in dislocation studies for a point dislocation. Geophys. Res. Lett. 29 (12) (46-1-46-4) doi: 10.1029/2001GL014497. Shen, Z.K., et al., 2009. Slip maxima at fault junctions and rupturing of barriers during the 2008 Wenchuan earthquake. Nat. Geosci. 010, 718–724, 10.1038/ngeo636. Somerville, P.G., Pitarka, A., 2006. Differences in earthquake source and ground motion characteristics between surface and buried earthquakes. In: Proceedings of the Eighth National Conference on Earthquake Engineering, San Francisco, California. Song, S., Pitarka, A., Somerville, P., 2009. Exploring spatial coherence between earthquake source parameters. Bull. Seismol. Soc. Am. 99 (2), 564–2571. Song, S.G., Somerville, P., 2009. Depth of faulting in large strike-slip earthquakes, final report. http://earthquake.usgs.gov/research/external/reports/ 08HQGR0084.pdf. Song, S., Somerville, P., 2010. Physics-based earthquake source characterization and modeling with geostatistics. Bull. Seismol. Soc. Am. 100, 482–496. Tang, C., Zhu, J., Qi, X., Ding, J., 2011. Landslides induced by the Wenchuan earthquake and the subsequent strong rainfall event: a case study in the Beichuan area of China. Eng. Geol. 122, 22–33, 10.1016/j.enggeo.2011.03.013. Thomas, A.L., 1993. Poly3D: a three-dimensional, polygonal element, displacement discontinuity boundary element computer program with applications to fractures, faults, and cavaties in the earth’s crust (M.S. thesis). Stanford University, pp. 221. Toda, S., Stein, R.S., Reasenberg, P.A., Dieterich, J.H., Yoshida, A., 1998. Stress transferred by the 1995 Mw = 6.9 Kobe, Japan, shock: effect on aftershocks and future earthquake probabilities. J. Geophys. Res. 103 (24), 543–624 (565) 10.1029/98JB00765. Toda, S., Stein, R.S., Sagiya, T., 2002. Evidence from the AD 2000 Izu Islands swarm that stressing rate governs seismicity. Nature 419, 58–61, 10.1038/nature00997. Tong, X.P., Sandwell, D.T., Yuri, F., 2010. Coseismic slip model of the 2008 Wenchuan earthquake derived from joint inversion of interferometric synthetic aperture radar, GPS, and field datam. J. Geophys. Res. 115, B04314, 10.1029/ 2009JB006625. Xu, X., et al., 2009. Coseismic reverse- and oblique-slip surface faulting generated by the 2008 Mw7.9 Wenchuan earthquake, China. Geology 37, 515–518, 10.1130/G25462A.1.

Wang, Qi, et al., 2011. Rupture of deep faults in the 2008 Wenchuan earthquake and uplift of the Longmen Shan. Nat. Geosci. 4 (9), 634–640, 10.1038/ NGEO1210. Wang, X.W., Liu, G.X., Yu, B., Dai, K.R., Zhang, R., Chen, Q., Li, Z.L., 2014. 3D coseismic deformations and source parameters of the 2010 Yushu earthquake (China) inferred from DInSAR and multiple-aperture InSAR measurements. Remote Sens. Environ. 152, 174–189. Wessel, P., Smith, W.H.F., 1998. New, improved version of the generic mapping tools released. Eos Trans. Am. Geophys. Union 79, 579, 10.1029/98EO00426. Weston, J., Ferreira, A.M., Funning, G.J., 2012. Systematic comparisons of earthquake source models determined using InSAR and seismic data. Tectonophysics 532–535, 61–81. Working Group of the Crustal Motion Observation Network of China Project, 2008. The coseismic displacement field of 2008 Ms8.0 Wenchuan earthquake measured by GPS. Sci. China Ser. D 38 (10), 1195–1206 (in Chinese). Wu, J.C., Tang, H.W., Chen, Y.Q., Li, Y.X., 2006. The current strain distribution in the North China basin of eastern China by least-squares collocation. J. Geodyn. 41, 462–470. Yang, Y.H., Chen, Q., Liu, G.X., et al., 2014. Correction of coseismic deformation field associated with Wenchuan earthquake with GPS observables and InSAR adjacent track smoothing and fault slip inversion. Chin. J. Geophys. 57 (5), 1462–1476, 10.6038/cjg20140511 (in Chinese). Zhang, Y.Q., Ma, Y.S., Yang, N., Shi, W., Dong, S.W., 2003. Cenozoic extensional stress evolution in North China. J. Geodyn. 36 (5), 591–661. Zhang, G.H., Qu, C.Y., Shan, X.J., et al., 2011. Slip distribution of the 2008 Wenchuan Ms 7.9 earthquake by joint inversion from GPS and InSAR measurements: a resolution test study. Geophys. J. Int. 186, 207–220. Zhao, C.P., Chen, Z.L., Zhou, L.Q., Li, Z.X., Kang, Y., 2010. Rupture process of the 8.0 Wenchuan earthquake of Sichuan, China: the segmentation feature. Chin. Sci. Bull. 55 (3), 284–292. Zhao, C.Y., Zhang, Q., Yang, C.S., Zou, W.B., 2011. Integration of MODIS data and short baseline subset (SBAS) technique for land subsidence monitoring in Datong, China. J. Geodyn. 52, 16–23. Zhu, A.L., Xu, X.W., Diao, G.L., Su, J.R., Feng, X.D., Sun, Q., Wang, Y.L., 2008. Relocation of the Ms8.0 Wenchuan earthquake sequence in part: preliminary seismotectonic analysis. Seismol. Geol. 30 (3), 759–767 (in Chinese). Ziv, A., 2012. Inference of coseismic slip via joint inversion of GPS and aftershock data: the 2004 Parkfield example. J. Geophys. Res. 117, B03307, 10.1029/ 2011JB008400.