Slip distribution of the 2008 Wenchuan M s 7 . 9 earthquake by joint inversion from GPS and InSAR measurements : a resolution test study

SUMMARY 
 
We investigate the slip distribution of the 2008 May 12 Wenchuan Ms 7.9 earthquake using GPS data and InSAR measurements under a linear inversion scheme, with emphasis on the effect of three factors, including constraint on rake, different discretizations, and layered elastic model. Within our inversion parameterization context, we find the most influential factor would be constraint on rake. Without constraint on rake, the slip model seems physically wrong under the depth of 15 km, due to the limited depth resolution of the geodetic data used, especially the one orbit of InSAR measurements. Thus it is necessary to add a priori to the slip rake to obtain a reasonable fault source model. Different discretizations of the subfault patches have a notable impact on the slip distribution. Also, the layered elastic model predicts more slip at depth than does the half-space model, by about 15–20 per cent. The characteristics of slip distribution established through our inversions include some points as follows: (1) the most reliable results would be expected at depth range of 0–15 km; (2) three peak-slip asperities are inverted, at Yingxiu county (near epicentre), Yuejiashan county, and Beichuan county, respectively; (3) the inverted rake along the Yingxiu-Beichuan fault changes from dominantly thrusting motion at the southwest segments to dominant right-lateral or even pure right-lateral strike slip at the northeast segments, and (4) only thrust slip is occurred on the Guanxian-Jiangyou fault.


I N T RO D U C T I O N
The 2008 May 12 M s 7.9 Wenchuan earthquake is one of the most devastating natural disasters that ever occurred in Chinese history and is another great earthquake since the 2001 M s 8.1 Kokoxili earthquake in China mainland. The event occurs at a northeast-southwest striking and northwest dipping fault zone, known as the Longmenshan fault (LMSF) zone, which is on the margin between the Tibetan plateau and the Sichuan Basin (Burchfiel et al. 2008). Along the Yingxiu-Beichuan Fault (YBF), known as the central fault of LMSF zone, an over 240-km-long surface rupture is produced by the main shock; another over 70-km-long surface rupture appears on the Guanxian-Jiangyou Fault (GJF), the frontal fault of LMSF zone Xu et al. 2009). Seismological studies exhibit a dominant thrust slip with some right lateral strike slip of the main shock, of which the rupture processes are composed of two to four subevents, lasting for about 80-120 s and propagating unilaterally from the southwestern end to the northeastern end of the ruptures (Table 1; Wang et al. 2008;Zhang et al. 2009). Unlike the 1999 Chi-chi, Taiwan earthquake, which occurred on a slow dipping (dipping 20-25 • ) reverse zone (Jonhson et al. 2001), the 2008 Wenchuan earthquake happened on a high dipping (∼50 • ), oblique, listric-reverse faulting zone ). These previous studies indicate an unusual and very complicated source mechanism of the Wenchuan earthquake, of which the rupture processes and the slip distribution may vary significantly in time and space.
The fault slip distribution is fundamental for further studies of the Wenchuan earthquake and the LMSF zone as well, such as: (i) as a constraint for the time-dependent inversions using seismic data (Delouis et al. 2002); (ii) for estimating the coseismic stress changes and stress transfer (Parsons et al. 2008;Toda et al. 2008); (iii) characterizing the seismic cycle (Shen et al. 2009). After the Wenchuan earthquake, numerous studies have been done in terms of rupture process and slip distribution inversion. Some seismologists have presented quick but preliminary slip models based on a single fault model (e.g. Ji & Hayes 2008;Chen et al. 2008;Zhang et al. 2009). More delicate seismic data analysis is done by Nakamura et al. (2009). Geodetic data inversion provides another independent way to explore the source model. Hao et al. (2009) has used the ALOS PALSAR data to achieve the slip distribution along a fault model composed of two fault planes. Shen et al. (2009) has used GPS and InSAR data, including the ALOS PALSAR data and Envisat ASAR images, to achieve simultaneously the fault geometry and the coseismic slip distribution within a nonlinear inversion  scheme. Recent studies done by Feng et al (2010) and Tong et al (2010) have also used the GPS and ALOS PALSAR data to achieve the slip distribution. Table 2 gives a brief summery of these studies. It is well known that an inversion problem involves some objective decision makings, which may introduce some uncertainties to the inversion results (Beresnev 2003). Then it seems necessary to do a resolution test in order to see how robust the inversion could be. In this article, we apply a scheme of linear inversion to invert the available geodetic data (GPS and InSAR), but put our emphasis on evaluation of three aspects, including the effect of boundary constraints on the slip rake, the effect of different discretizations of the subfault patches, and the effect of a layered elastic model.
We obtain the coseismic deformation measured through interferometry of SAR images acquired before and after the Wenchuan event. The earthquake rupture is then simulated using five subfaults. The segmentation of the fault geometry is simplified by analysis of the ALOS PALSAR offset map. A grid searching method is then applied to explore the possible variability of dip angle among each segment. After the determination of the fault geometry, we apply our first resolution inversion, focusing on the effect of a boundary constraint on rake. In our second inversion test we try to examine the effect of different discretization schemes, that is, making the subfault element uniform patch or resolution-based patches. Further comparison is carried out between a homogeneous elastic half-space model and a horizontally layered elastic model, to see the effect of assumed layered elastic model.

InSAR measurements
InSAR is a widely used technique for measuring large-scale crustal displacement and has been successfully applied to many earthquakes (e.g. Massonnet et al. 1993;Feigl et al. 1995Feigl et al. , 2002Funning et al. 2007). The Advanced Land Observing Satellite (ALOS) PALSAR sensor uses the microwave in L band with a wavelength of about 23.6 cm, which could ensure a fine coherence of SAR images, even in heavily vegetated areas (e.g., Sandwell et al. 2008;Shimada e al. 2008). We choose PALSAR raw data in seven tracks, including track 471,472,473,474,475,476 and 477, which could cover almost the entire coseismic deformed area, including Yingxiu (macro epicentre of the main shock), Dujiangyan, Wenchuan (microepicentre of the main shock), Maoxian, Beichuan, Pingwu and Qingchuan (Fig. 1). Each track contains 16 frames of data, of which eight frames were acquired before the Wenchuan earthquake and the other eight frames after the earthquake.
The conventional two-pass differential InSAR method and the Gamma software are used to process the ALOS PALSAR data. The frames of raw data in each track are connected above all and then transform into single-look images (SLCs). Then seven interferograms are generated through interferometry of the SLC images ( Fig. 1). We use the SRTM 3 s data as the external DEM. Topographic and orbital geometric phases are then estimated and erased from interferograms. We use the minimum cost flow algorithm (MCF, Werner et al. 2002) to unwrap the interferograms. Considering the fact that in track 472, 473, 474, 475 the phases on each side of the fault are separated by an incoherent zone, we set two start points at the northernmost and southernmost edges of the interferograms, respectively and apply the MCF unwrapping separately. For the rest three tracks of interferograms (track 471, 476 and 477), which have continuous phases, the start point is set at a highly coherent point at the southernmost edge of each interferogram. Low coherence area (coherence under 0.35) is masked out and will not be unwrapped. Finally, seven tracks of displacement images are geocoded at 90 m spacing and jointed into a large-scale displacement map.
Inversion of the InSAR data at its full resolution is an impractical and unnecessary task (Simons et al. 2002;Lohman et al. 2005). We take two steps to reduce the data and get a more manageable data for inversion. A three-by-three down-sampling is used to reduce the data points (geocoded pixel size of 270 m) first. After this averaging, the data points are decreasing from billions to millions. Further averaging will inevitably hurt the resolution too much, especially to the near field (Simons et al. 2002). We then employ the quad tree decomposition algorithm (Jónsson et al. 2002), which focus on the gradient of the LOS displacements. A 0.1 gradient value has been chosen, i.e. if a block in which the displacement gradient is bigger than 0.1, the block will be divided into four blocks until the gradient value is smaller than 0.1. In a final undivided block, the data point value is the average of all the data points included. This averaging will help to control noises to some extent. At the same time, a masking method has been applied. Those data points with a lower coherence than 0.5 will not be resampled and hence not be used in the process of inversion. Finally, we obtain 3500-4500 of data points for each track, with a data set of almost 25 000 data points totally (Fig. 2).

Slip distribution of 2008 Wenchuan earthquake
3 Figure 1. Local settings and coseismic deformation of the Wenchuan earthquake derived from InSAR data. The red star and red lines are the epicentre of the Wenchuan earthquake and regional faults, respectively. The upward-pointing arrow shows the flight direction of the ALOS satellite, and the rightward-pointing one is the line of sight (LOS). Background topographic data is from the SRTM 3. LMSF denotes the Longmenshan fault zone; EKLF is the Eastern Kunlun fault; MJF is the Minjiang fault; XSHF is the Xianshuihe fault.
We also use offset tracking technique to estimate displacement in both the range and azimuth directions. The offset fields are generated with a normalized cross-correlation of image patches of detected real-valued SAR intensity images. A pixel searching window of 256 is selected in both the range and azimuth directions and is then moved by 128 pixels. Finally, we produce maps of range and azimuth offset (Fig. 3). Offset tracking can provide the displacements in decorrelated regions, but with a lower precision of 15 cm or worse, compared the precision of 2-3 cm of InSAR measurements (Furuya et al. 2010). Instead of using it directly as an input to our inversion, we use it as a constraint to simplify the fault ruptures (Funning et al. 2007).

GPS coseismic displacements
A quantity of GPS observation sites was deployed on both sides of the LMSF, primarily for campaign observations, for which the duration was 3-4 d. The last term measured before the Wenchuan earthquake was conducted between April and July of 2007 and it was resurveyed immediately after the earthquake. The data has been processed using the GAMIT/GLOBK GPS processing software (King & Bock 2004). We refer to the original article for details relating to the processing of the GPS data . The coseismic GPS data illustrates the state of crustal deformation caused by the Wenchuan earthquake, which appears as an eastward motion on the hanging wall and a westward motion on the footwall of the fault. In general, the GPS coseismic deformation demonstrates a relative motion between the LMSF zone and Sichuan basin (Fig. 2). 122 GPS horizontal and 33 vertical displacements are adopted to conduct our inversion.

Inversion method
A constrained least-square method to determine the slip vector is introduced here. The rms is defined as a measure of the difference between the observed and the modelled, while the data-model correlation measures the linear correlation between the observation and the prediction by model. In order to avoid unphysical fluctuation of the slip distribution, we include a smoothing term in the misfit function. A smoothing factor (arbitrarily normalized) of 0.24 4 G. Zhang et al. is determined by analysing the trade-off between the squared data misfit and the squared slip roughness (Fig. 4a). In addition, a static offset in each InSAR data set is inverted as an unknown and used to correct the InSAR. The mathematical formulation of the joint inversion is expressed by where i stands for GPS or InSAR measurements, U i and W i.ratio are the data uncertainties and the relative weight between GPS and InSAR, respectively. The uncertainty of InSAR data is set to be 50 mm and the uncertainty of GPS data is set to be their original uncertainties. s is slip vector, d is the matrix of observation data, d 0 is the static offset for InSAR measurements, M is the Green's function for an elastic half-space or layered elastic halp-space, which describes the relation between the model prediction and the observation. We calculate the Green's functions of the homogeneous elastic halfspace model using Okada (1992), assuming a Poisson ratio of 0.25 and construct the Green's functions of the layered elastic model using EDGRN program (Wang et al. 2003). For the last term, β 2 is defined as the smoothing factor, H is the Laplacian operator, and Hs 2 represents a measure of the slip roughness. It is still subject to determine the relative weighting between GPS and InSAR data (Shen et al. 2009). In order to achieve a desired weighting, we use the data misfit as a measurement to select the weight (Fig. 4b). The basic idea is to give more weight to the data  set that could minimize the misfit function. Our favoured weight of InSAR relative to GPS is 0.1, which means one GPS station equals ten InSAR data points. No further weighting within the GPS data is performed, but the down-sampling method to pick InSAR data can be deemed as a relative weighting within InSAR data (Simons et al. 2002).

Fault segmentation and dip angles
We assume nine parameters to describe the fault dislocation model, including strike, dip, slip, rake, starting longitude, starting latitude, depth, length and width. The solution of these nine parameters simultaneously is a non-linear inversion problem (Yukitoshi & Wright 2008). To solve this high non-linear problem, we make the segmentation of the fault model through analysing offset tracking of ALOS PALSAR images at first (Fig. 3). The offset map gives a strong constraint on fault segments (Simons et al. 2002;Funning et al. 2007;Shen et al. 2009). The location, length and strike of each segment are estimated from the offset maps and adjusted by running a few preliminary inversions. A simplified fault model with five segments is chosen to simulate the earthquake ruptures (Fig. 5). Four segments are used to delineate the main rupture, that is, the Yingxiu-Beichuan fault (YBF) and one segment for the secondary rupture, the Guanxian-Jiangyou fault (GJF). The southwest end of the fault is set to start at ∼20 km southwest of Yingxiu where the epicentre is located. The northeast end of the fault is inferred from the InSAR interferograms at around southwest of Qingchuan, where the aftershocks are abundant but the rupture may not have broken the surface ( Fig. 5; Shen et al. 2009). For the base of the segments, we set the YBF to be at 32 km depth and the GJF to be at 18 km.
The fault dip is difficult to determine. Previous studies of aftershock relocation could not give the dip angle in a strict sense (Huang et al. 2009;Chen et al. 2009). There is also still a lack of enough geologic profiles, which could delineate the fault geometry in 3-D view (Burchfiel 2008). We then test the sensitivity of the weighted normalized RMS (WNRMS) to the fault dip of each segment with the help of a grid searching method, similar to that of done by Tong et al. (2010).
As shown in Fig. 6, the fault dip is in a wide range of 40 • -80 • . The reason for the difference between our findings and those of Tong is that they invert the ascending and descending InSAR measurements simultaneously, which could have a better constraint than just using the ascending data. However in any of the two cases, the dip angle could not be determined by the current geodetic measurements with much confidence. Before settling down on our preferred dip angle, we run a few inversions to test the listric-shaped dip angles within each segment and don't see an obvious decrease in the WNRMS and the spatial slip distribution doesn't show much difference either, in contrast to the constant dip angle within each segment. This seems to be a clue that the current InSAR data and the sparsely distributed GPS measurements are insensitive to the variation of the dip angle along depth, as also found by Furuya et al (2010). The fault dip of each segment in our final fault model has a moderate dip angle (47 • ) at southwest end of the rupture, becomes gradually steeper northeastward and is nearly vertical at the northernmost segment (80 • ). This is in a general consensus between other studies Shen et al. 2009;Tong et al. 2010). Detailed fault parameters can be found in Table 3. After doing so, we fix all of the nine parameters except slip and rake.

With or without boundary constraint on rake
We firstly compare the inversion results, considering whether applying a constraint on rake or not. We make the boundary constraint on rake based upon seismological studies, which has shown that the Wenchuan event is dominantly thrusting with some dextral strike slip ). Therefore the initial value of the rake is confined within the range of 90-180 • . Afterwards, we relax the constraint on rake by step of 90 • . In order to examine its effect to inversion, we keep all other inversion settings as the same during these schemes of inversion, such as the dimension of the subfault 6 G. Zhang et al. element (approximately 4 km along both strike and dip directions), smoothing coefficient (0.24) and the Poisson ratio (0.25, i.e. a single layered homogeneous earth model is adopted here). Fig. 7 shows the inversion result from making the strong constraint on rake within the range of 90-180 • . We notice that the rake reaches their boundary values on some of the subfaults (180 • on YBF at depth greater than 15 km and 90 • on GJF). In order to verify this boundary constraint, we relax it by step of 90 • . Here we only show the inversion results without any boundary constraint on rake (Fig. 8). When we relax the upper bound value of the rake from 180 • to 360 • by step of 90 • , the slip motion exhibits normal faulting or normal plus small left lateral striking on some parts of the YBF at depth greater than 15 km and on the whole GJF, which seems physically wrong (Fig. 8). When we relax the lower bound of the rake from 90 • to 0 • , there is almost no change on the entire YBF, whereas the GJF shows dominant thrusting and small left lateral striking at shallower depth and almost pure left lateral motion at greater depth, which also seems physically improper (Fig. 8).
There are at least two explanations for the unphysical slips. First of all, the normal faulting artifacts on the YBF at greater depth are caused by the so-called InSAR measurements ambiguity. In-SAR measures only the range changes along LOS and one orbit of InSAR measurements is not sufficient to determine the real 3-D deformation, which leads to relative limited depth resolution in terms of slip distribution inversion. We conduct a series of synthetic resolution tests to confirm our conclusion. Using the input source model from the real data inversion with rake constraint from 90 • to 180 • (Fig. 7), we compute the geodetic synthetics at their original observation position, but along both ascending and descending orbits. Due to the seriously incoherence in Longmenshan mountain area, we do not process the descending SAR images. However by forward modelling, we could synthesize the descending measure-ments and use it in our resolution test. Before running the synthetic inversions, we also add 10 cm noise to account for the possible atmospheric effect in InSAR measurements. All the parameters for synthetic inversion are kept as the same as for the real data inversion, except that we remove the boundary constraint on rake. We first invert the GPS and the ascending InSAR synthetics (which we have their counterpart observations). The resolution inversion results show that with sparsely distributed GPS stations and the ascending InSAR measurements only, the input source model cannot be successfully recovered. The output slip motion is still physically wrong along the YBF at depth and the entire GJF, which is very similar to the real data inversion without constraint on rake (Fig. 8). Then we invert all the geodetic synthetics, including GPS, ascending and descending InSAR. This improved synthetic inversion successfully recovered almost all the input slip along the YBF, but the input slip on the GJF remains poorly resolved at depth greater than 10 km (Fig. 9). These schemes of synthetics prove that the depth resolution of slip is mainly related to the InSAR LOS ambiguity. With a proper constraint on rake we can obtain an optimal slip distribution. We also argue that the resolution of one orbit of InSAR measurements is most probably confined within the upper 15 km crust and two orbits of InSAR measurements can efficiently reduce the LOS ambiguity and thus have better depth resolution in slip source inversion.
Second of all, the interference of the two subparallel ruptures should be responsible for the unphysical left strike slip along GJF. In case of the Wenchuan earthquake, even two orbits of InSAR measurements plus a certain number of GPS data could not have a satisfying resolving ability along the secondary rupture of the GJF, which is nearly parallel to the main rupture of the YBF. Also there is no observation data available between the YBF and the GJF, which makes the slip partitioning even more difficult to determine. However, no matter what kind of constraint on rake applied, some common characteristics of the slip distribution can be established. First of all, the slip direction always shows dominantly thrusting and some dextral strike slip at southwest segments of the YBF, then changes into nearly equally thrusting and dextral striking at middle segments of the YBF, and finally exhibits almost pure dextral strike slip at northeast segments of the YBF at depth range of 0-15 km. Three main asperities are inverted, located at near Xingyiu county, Yuejiashan county and Beichuan county, respectively. We also find that the data fitting among all the inversions is almost in the same level and the relaxing of constraints on rake does not improve the data fitting level dramatically (Table 4 and Fig. 10).

Uniform-based discretization and resolution-based discretization
We use a uniform-based discretization in our former schemes of inversion, which means that the subfault patches remain as the same throughout the fault model. We choose a commonly-used value of approximately 4 km along both strike and dip direction (Shen et al. 2009;Xu et al. 2010). In this section we will show the results from different discretizations, including a uniform-based discretization but with a smaller subfault element of approximately 2 km (Figs 11c and d) and a resolution-based discretization (Figs 11e and f), for which the subfault element dimension increases progressively with  depth by a factor of 1.5 (Tong et al. 2010). We apply the boundary constraint on rake from 90 • to 180 • in all the inversions hereafter. The results from uniform-based discretizations with different patches (2 or 4 km) have a similar pattern of slip distribution, which show two major thrusting concentrated slip areas and three major right lateral strike slip areas at shallower depth (Figs 11a-d). However two big discrepancies do occur between these two inversions of uniform subfault patches. One is the maximum slip. The 4 km discretization has a 7 m maximum value both for thrusting slip and dextral strike slip, but the 2 km discretization reaches 11 and 9 m for thrusting slip and dextral strike slip, respectively. Other dimensions of the subfault patch between 4 and 1 km are tested and we find that, the more tightened of regularization, the larger of the maximum slip is required for the geodetic data to converge to a satisfying solution. The other noticeable difference is that the 2 km discretization scheme of inversion has more concentrated slip distribution, which makes the slip at shallower depth separated from the slip at greater depth, whereas the slip distribution of the 4 km discretization inversion is connected between the shallower depth and the greater depth (Fig. 7).
In the resolution-based discretization inversions, we set the first row and column of the subfault patch to be 2 km and others to increase progressively with depth by a factor of 1.5 (Tong et al. 2010). Other values of the factor (from 1 to 2) and the first row and column dimension (from 1 to 2 km) are also tested and they show only minor differences among them. The results from resolutionbased discretization inversion have a more scattered pattern of slip distribution, in contrast to the uniform-based discretizations (Fig. 11). However, we think the slip should be more precisely located at shallower depth in this inversion because of the consistency between the data resolution and the discretization strategy of the subfault patches. Similar to those of the results from all the former inversions, this scheme of inversion has better resolution at shallow depth. The slip distribution forms two major asperities of thrusting motion at southwest segments of the YBF at depth range of 0-15 km (Fig. 11). Significant dextral strike slip occurs at northeast segment of the YBF at near surface (Fig. 11). In the mean time, the slip on the GJF is found to be mainly thrusting and almost no dextral strike slip.
We choose to ignore the slip at depth greater than 15 km, as they might be artefacts caused by the geodetic data. This could be seen from the synthetic inversions, which show some unphysical slip at depth using the ascending orbit synthetics only and almost no unphysical slip at depth using both the ascending and descending orbits of InSAR synthetics (Fig. 9). Another thing worthy to notice is that greater depth is required for dextral strike slip than for thrusting slip. We think this is also related to the LOS ambiguity of the InSAR measurements, which has high precision in vertical but poor precision in horizontal.

Layered elastic model
The final inversion we applied is to introduce the effect of a layered elastic model. The layered crust model is compiled from profiles of 10 G. Zhang et al. Figure 9. Slip distribution inverted from synthetic data without any constraint on rake. Geodetic synthetics include two orbits of InSAR and GPS. Note that the colourbar is the same as in Fig. 7. See details in text. crust 2.0 by Bassin et al (2000). Detailed parameters are listed in Table 5.
We compare inversion results from assuming an elastic halfspace model and a layered elastic model, using a resolution-based discretization. Slip distribution from these two elastic models has a similar scattered pattern, although we find the layered-space model produces more slip than does the half-space model. Same conclusion is obtained by Simons et al (2002) studying the fault slip of the 1999 M w 7.1 Hector Mine Earthquake. They claim that the main effect of the rigidity stratification is to increase the inferred slip at depth by about 20-30 per cent. Here in our example, the increase seems to basically distribute at depth range of 0-15 km and it increases by 15-20 per cent.
The seismic moment estimated from different inversions varies from 6.6 × 10 20 to 9.4 × 10 20 N m ( Table 2). The largest estimation of the scalar moment is from inversion with no boundary constraint on rake, which is also the best-fitting model of the data (Table 4). However as we have already seen before, the fault model is actually not well resolved during that scheme of inversion. We think the average value of 8.0 × 10 20 N m is the most robust estimation, which is also the value from resolution-based discretization of inversion in a layered elastic earth model. Assuming a shear modulus of 3.0 × 10 10 Pa, the moment magnitude of the Wenchuan earthquake is equivalent to M w 7.9.

Comparison with other fault slip models
Numerous studies have been carried out in terms of slip distribution and rupture process of the Wenchuan earthquake. Basically two types of data sets are used to achieve the source model, including teleseismic data and geodetic data. The first documented source model is from Ji & Hayes (2008), which is released on the USGS website only 7 hr after the earthquake. Based on a single fault model, they invert 17 teleseismic broad-band P waveforms and 10 broad-band SH waveforms for the rupture history of the earthquake. The final static slip distribution have two asperities with maximum slip of 10 m, one located at 50 km northeast of the hypocentre at depth of 10 km, the other at 180 km northeast of the hypocentre at depth of 14 km. One major difference between this source model and most of others is that the slip distribution gives the asperities at depth greater than 10 km. Using P teleseismic waveforms only and under a 500-km-long single fault model, Zhang et al (2009) also does a similar job and finds two concentrated slip areas, one distributed in a larger space near the epicentre and the other at near  Beichuan county. They locate the maximum slip area at near surface, which is different from Ji & Hayes (2008) but similar to those of geodetic findings, including ours. These source models are quick but very preliminary. Most of other studies have proved that the fault geometry of the Wenchuan earthquake may vary significantly in space and one single fault segment cannot represent the complexity of the earthquake. Nakamura et al. (2009) does a more delicate study also using the teleseismic waveform. By separating the fault into southwestern and northeastern segments, they conclude that at least two types of striking mechanisms should be considered in the Wenchuan earthquake: one is the thrusting motion mechanism in a low dip angle at southwestern segment; the other is the dextral strike slip mechanism in a high dip angle at northeastern segment. In all of our inversion results, we have also successfully pointed out that the possible mechanism transition between southwest segments and northeast segments of the YBF. Teleseismic inversion has poorer depth resolution, especially in contrast to the spatially dense geodetic data (Tong et al. 2010). There is a rich geodetic data set of the Wenchuan earthquake, including a certain number of GPS and plenty of InSAR measurements. Using the ALOS ascending InSAR measurements, similar to those of used in our study, Hao et al. (2009) presents a source model along two fault planes and yields that thrust fault-slips are dominant in the southwest segment, while dextral fault-slips are dominant along the northeast segment. For comparison, our results also show large dextral strike slip along the southwest segments of the YBF at near surface, in addition to the significant thrust slip. This is most probably caused by the inclusion of GPS data in our inversion. Based upon two orbits of InSAR data (the ALOS PALSAR ascending data and the Envisat ASAR descending data) and GPS data, Shen et al (2009) achieves the fault geometry and the coseismic slip distribution simultaneously within a nonlinear inversion scheme. The fault geometry they obtained has a moderate dip angle of ∼43 • at the southwest end and gradually becomes steeper northeastward. This is in a general agreement with our fault model. Due to the similarity of the adopted data sets and the fault geometry as well, similar slip distribution characteristics can be established between the results from our resolution-based discretization and those of Shen et al (2009). Two major asperities of thrusting motion are located in near epicentre area and Beichuan area at southwest segments of the YBF and most of the dextral strike slip only occurs at northeast segments of the YBF at depth range of 0-15 km. Moderate thrust slip occurs on the GJF. However the results from theirs do not have large slip at depth under 15 km as our inversion results do. Considering that they use both ascending and descending orbits of InSAR measurements, we think the ambiguity of the slip at depth could be successfully reduced along the YBF. Feng et al (2010) processes the ALOS PALSAR data that is unlikely contaminated by the possible ionospheric effect before the Wenchuan event, except for track 475 that no uncontaminated data is available. They use four segments to simulate the earthquake ruptures, three segments for the YBF and one segment for the GJF. We find enormous similarities of slip distribution between their results ( fig. 2 in their article) and the results from our uniform-based discretization with constraint on rake (Figs 7, 11a and b). As they also apply a strong constraint on rake, both of the slip directions show dominantly thrusting with small dextral strike slip at southwest segments of the YBF, yet change into nearly equal thrusting and dextral striking at middle segments of the YBF, and finally exhibit almost pure dextral strike slip at northeast segments of the YBF. We both locate three main asperities at near Yingyiu county, Yuejiashan county, and Beichuan county, respectively. The main difference seems to occur at the middle asperity near Yuejiashan county, where our results show much larger slip than theirs. This is most probably caused by the possible ionospheric effect of the SAR images. Tong et al (2010) uses the most comprehensive of data sets so far, including ascending and descending InSAR measurements, GPS and surface scarp data as well. However they do not unwrap the InSAR interferograms as they think it is too difficult to do so (Tong et al. 2010). To avoid the unwrapping problems they digitize the fringes manually and convert them to displacement in the LOS direction. Though we all adopt a resolution-based discretization of the subfault, the slip distribution shows some discrepancies at greater depth than 15 km, due to the big difference in the data sets we used. They obtain a much smoother slip 12 G. Zhang et al. Figure 11. Slip distribution inverted from different discretizations. Left (a, c and e) are thrust slip from the 4 km uniform-based discretizaiton, the 4 km uniform-based discretizaiton and the resolution-based discretization inversion, respectively. Right (b, d and f) are dextral strike slip from the 4 km uniform-based discretizaiton, the 4 km uniform-based discretizaiton and the resolution-based discretization inversion, respectively. distribution even at depth, owing to the combination of two InSAR orbits.

C O N C L U S I O N S
Using geodetic data, we conduct a series of inversions, each with different parameter settings. Our main purpose is to evaluate the uncertainties of the slip distribution caused by some subjective decision-makings during inversion, which is well known in seismic data inversion but is not so clear in geodetic data inversion. Our conclusions are as follows.
(1) A proper constraint on the slip rake is recommended to acquire a more physically reasonable slip model, due to the resolving ability at depth and the LOS ambiguity of InSAR measurements. We find that the inversion results without an appropriate constraint on rake exhibit a physically wrong slip model, especially at depth under 15 km and the whole secondary rupture of the GJF. We thus claim the importance of a proper constraint on rake.
(2) Through synthetic inversion test we conclude that the slip of the YBF (the main rupture of the Wenchuan earthquake) could be well resolved, especially at depth range of 0-15 km. However the slip on the secondary rupture of the GJF is difficult to determine, due to the slip partitioning between the two subparalleled ruptures.
(3) Different discretizations of the subfault patches will result in some discrepancies. With a more tightened and uniform-based discretization, smoother slip distribution and larger maximum slip value are expected. With a resolution-based discretization, the slip distribution might be scattered and the slip should be better located at shallow depth.
(4) Assumed a layered elastic model would predict more slip at depth than that of an elastic half-space model. Data fitting of a layered elastic model is somewhat higher than that of an elastic half-space model. The computed magnitude also seems to be more consistent than those of the elastic half-space model. This gives us a clue that a layered elastic model may be needed in terms of slip distribution inversion, though we are unable to consider in our model any lateral heterogeneities existing in the Wenchuan earthquake, as revealed by seismic studies Wu et al. 2009).
(5) Neither constraint on rake, nor different discretizaiton makes the data fitting improving dramatically and the residuals among all of our inversions seems to be at the same level and it is one order higher than a simple rupture earthquake such as the 1997 M w 7.6 Manyi earthquake in 1997 (Funning et al. 2007). This is partially because of the complicated rupture processes of the Wenchuan earthquake. Error factors may include atmosphere vapour, topographic error, ionospheric effect, aftershocks and possible post-seismic deformation.
(6) Through comparisons with other studies, we find some common features of slip distribution among all of the schemes of inversion. Three peak-slip areas along the YBF are inverted, at Yingxiu, Yuejiashan and Beichuan areas, respectively; and the inverted rake shows dominantly thrusting at southwest segments and changes progressively to dominantly dextral striking or even pure dextral striking northeastward.
(7) The slip distribution is heavily depended on data used, fault geometry and the inversion method. In this work, we used geodetic measurements, including GPS and InSAR to simulate the spatial slip distribution of the Wenchuan earthquake with a linear inversion method. Due to the complex faulting and rupture process of the Wenchuan earthquake, the slip distribution cannot be taken for granted for all of the details, but only allow gross characteristic analysis. Further study on the slip distribution inversion should be carried on to obtain a more robust slip distribution by joint inversion of all available data sets, including seismic data (both teleseismic data and near-field strong motion data), geodetic data and surface scarps as well.

A C K N O W L E D G M E N T S
Authors thank the editor and two anonymous reviewers for comments and suggestions. This work is co-supported by grants of State Key Laboratory of Earthquake Dynamics, China Earthquake Administration (LED2008A06, LED2010A02), Basic Scientific Funding of Institute of Geology, China Earthquake Administration (IGCEA0908) and National Science Foundation (40874006). It is also supported by grant of the National Science Council 100-2116-M-002-004-. JAXA owes the copyright of the PALSAR level 1.0 data. The topographic data is downloaded from: http://srtm.csi.cgiar.org. The Global Centroid Moment Tensor Project database is searched using www.globalcmt.org/CMTsearch.html (last accessed on 2008 September 12). Generic mapping tools (GMT4.3) are used to generate figures (Wessel & Smith 1998).