Coseismic slip in the 2010 Yushu earthquake (China), constrained by wide-swath and strip-map InSAR

On 14 April 2010, anMw = 6.9 earthquake occurred in the Yushu county of China, which caused ∼3000 people to lose their lives. Integrated with the information from the observed surface ruptures and aftershock locations, the faulting pattern of this earthquake is derived from the descending wide-swath and ascending strip mode PALSAR data collected by ALOS satellite. We used a layered crustal model and stress drop smoothing constraint to infer the coseismic slip distribution. Our model suggests that the earthquake fault can be divided into four segments and the slip mainly occurs within the upper 12 km with a maximum slip of 2.0 m at depth of 3 km on the Jiegu segment. The rupture of the upper 12 km is dominated by left-lateral strike-slip motion. The relatively small slip along the SE region of Yushu segment suggests a slip deficit there. The inverted geodetic moment is approximately Mw = 6.9, consistent with the seismological results. The average stress drop caused by the earthquake is about 2 MPa with a maximum stress drop of 8.3 MPa. Furthermore, the calculated static Coulomb stress changes in surrounding regions show increased Coulomb stress occurred in the SE region along the Yushu segment but with less aftershock, indicating an increased seismic hazard in this region after the earthquake.


Introduction
On 14 April 2010, a devastating (M w = 6.9) earthquake struck the Yushu county, south area of Qinghai Province, China.This event caused 2698 fatalities and more than 270 people missing (Daniell et al., 2011).The Yushu earthquake ruptured on the Ganzi-Yushu fault (Fig. 1), which strikes northwest and runs through Ganzi, Yulong, Dengke, Yushu, Dangjiang and Zhiduo counties from SE to NW, with a length of about 500 km and dips northeast with a steep angle (60 • -75 • ) (Zhou et al., 1996).Together with the Xianshuihe fault, the Ganzi-Yushu fault forms the northern boundary between the Chuandian and Qiangtang active blocks, and the southern boundary of the Bayan Har active block (Zhang et al., 2003).
The Ganzi-Yushu fault is a Holocene left-lateral strikeslip fault, as a part of the principal left-lateral strike-slip fault system accommodating eastward extrusion of crustal blocks of the eastern Tibetan Plateau, and the plateau was formed due to the collision between the Indian and Eurasian plates (Searle et al., 2011, and references therein).Zhou et al. (1996) reported that the average rate on the Dangjiang segment is 7.3 ± 0.6 mm a −1 based on historical surface rupture, 14 C and thermoluminescence dating.Based on the geomorphological studies, Wen et al. (2003) estimated an average left-lateral slip rate of 12 ± 2 mm a −1 for the past 50 ka, which is close to the slip rate of Xianshuihe fault.Gan et al. (2007) inferred a present average slip rate of ∼14.4 mm a −1 along the whole Yushu-Xianshuihe fault using regional GPS data.Using C-band Envisat images collected between 2003 and 2010, Liu et al. (2011) suggested the interseismic slip rate of the Ganzi-Yushu fault is 6.4 mm yr −1 .Since late Quaternary, the fault is seismically active, with many large earthquakes (Sun et al., 2012).
The fast inversions of rupture process from broadband seismic data by Liu et al. (2010) and Zhang et al. (2010) showed that the earthquake was composed of two sub-events,  (Track 487 and 139).The black-white beach balls are from USGS (http://earthquake.usgs.gov/earthquakes/eqinthenews/2010/us2010vacp) and GCMT (http://www.globalcmt.org),respectively.The red circles are aftershocks (M>2.0)from ISC (http://www.isc.ac.uk) till 16 December 2011.Blue vectors represent interseismic GPS velocities relative to a Eurasian reference frame (Gan et al., 2007).located near the hypocenter and a region to the southeast of the epicenter, respectively, with peak slip of 2.1 m.The event was a unilateral rupture mainly propagating southeastward.Field investigations (e.g.Rao et al., 2011;Guo et al., 2012;Li et al., 2012;Sun et al., 2012) indicated that the surface rupture of the Yushu earthquake was subdivided into three geometric sections, each of which was separated by significant steps and no rupture zone, and the total length of surface rupture is about 65 km with a left-lateral offset up to 2.4 m.With the coseismic surface displacement observed by Interferometric synthetic aperture radar (InSAR), models (e.g.Li et al., 2011;Tobita et al., 2011) revealed that the length of surface and subsurface faults is about 73 km with the maximum left-lateral slip of ∼2.6 m.Studies based on seismological, geological field and InSAR observations indicated that the Yushu earthquake created an ∼70 km long surface rupture zone along the Yushu fault, but there were significant differences in the maximum slip magnitude and spatial slip distribution among these studies.
Combined with ascending and descending tracks InSAR observation, the high-density near field geodetic measurements can provide strong constraints on the focal mechanism and slip distribution in depth, which is useful to estimate the Coulomb stress change on the surrounding faults and therefore important for the seismic hazard assessment.For the past InSAR studies of the 2010 Yushu earthquake, the Green function was only calculated in elastic half-space and the slip smoothing between neighboring fault patches was employed to avoid oscillating slip distribution (e.g.Li et al., 2011;Tobita et al., 2011).In this study, we combine Phase Array Lband Synthetic Aperture Radar (PALSAR) InSAR data collected by the Advanced Land Observation Satellite (ALOS) with both ascending strip mode and descending wide-swath mode, geological field observations, and aftershocks distribution, along with the Green function for a layered crustal model and stress drop smoothing constraint, to derive a refined coseismic slip model for the Yushu earthquake and to investigate the static stress changes in the surrounding crust affected by the mainshock.

InSAR data processing and coseismic deformation
In January 2006, the Japanese Aerospace eXploration Agency (JAXA) successfully launched the ALOS satellite, which carried the L-band PALSAR for acquiring strip-mode SAR imagery on ascending tracks and wide-swath/scanning synthetic aperture radar (ScanSAR) imagery on descending tracks.Compared with the C-band SAR, the PALSAR preserves high coherence in interferometry over regions of rugged topography with high rainfall and dense vegetation.Because the interferograms with different imaging geometries can provide more information about surface deformation (Wright et al., 2004;Y. Liu et al., 2012), a more reliable coseismic slip model can be retrieved from interferometry data from both ascending and descending orbits.After the Yushu earthquake, we obtained four ALOS PALSAR images from the JAXA to determine the surface deformation, fault rupture geometry and slip distribution.Our images include two ascending images (Path 487A) and two descending images (Path 139D) (Fig. 1, Table 1).It should be noted that the descending path 139D operated with the ScanSAR mode, which is composed of five sub-swaths and covers an area of ∼350-km wide.In this study, we only use sub-swath 2 and 3 to construct the ScanSAR interferograms.
All interferograms were generated from the PALSAR level 1.0 (raw data) images with two-pass differential InSAR method using the Switzerland GAMMA software (Werner et al., 2000).Effects of topography were removed from the interferograms using JAXA precision orbits together with 90-m resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) (Farr et al., 2007).The interferograms were then filtered using a power spectrum filter (Goldstein and Werner, 1998) to reduce the effects of phase noise, and the unwrapped effects with the branch cut method.The interferograms were finally geocoded to geographic coordinates.
For the unwrapped interferograms, there are obvious longwavelength signals, which may be attributed to the inaccuracy orbits with the long perpendicular baselines (∼1 km).In addition, the Yushu earthquake occurred in the northern Tibetan Plateau, where is usually very dry, but there are still some atmospheric propagation errors in the interferograms.To remove the residual orbital errors and topography-induced atmospheric delay (Elliott et al., 2008;Cavalie et al., 2008), we assumed that far-field displacement (dashed region in Fig. 2) is zero and that the residual phases can be expressed by a bilinear function z = ax + by + ch + d, where (x, y) are the azimuth and range in the radar coordinate system, h is the elevation.The parameters of a, b, c and d from each interferogram can be estimated using the least-square method.Figure 3 shows the interferograms with tropospheric delay and long-wavelength orbital errors (Fig. 2) corrected.
Clear interferometric fringes caused by the Yushu earthquake can be visualized in both ascending and descending interferograms (Fig. 3), and significant earthquake-induced displacements are distributed over an area of 90 km by 70 km.Lower coherence can be attributed to a combination of high displacement gradient caused by earthquake and landsides, temporal decorrelation from vegetation and seasonal snow, and inadequate topographic phase correction in the area of extreme relief.From Fig. 3 we can see the crustal deformations are concentrated in two areas, the Longbao Lake and the region of 20 km WN of Yushu County, where displacement jumps can be detected.The rupture length of Longbao Lake area is about 20 km with relatively small magnitude displacement, whereas the rupture length of the area of 20 km WN of Yushu County is 30 km with relatively large surface displacement.The maximum line-of-sight (LOS) displacements in ascending 487A and descending 139D are 0.37 m and 0.50 m, respectively.
Although the significant errors are removed from the interferograms, there are still some unmodeled errors remaining in the detrended interferograms.Because the residual errors of interferograms are highly spatially correlated, a 1-D covariance function (Hanssen, 2001) is used to describe the uncertainty characteristics of interferograms (Fig. 3), including the magnitude and spatial scale.The magnitude and spatial scale of errors are 2.2 cm and 8.5 km for descending 139D interferograms, respectively, and those of ascending 487A are 1.9 cm and 8.3 km (Table 1).

Inversion for coseismic slip model
The InSAR coseismic deformation map (Fig. 3) consists of millions of observations.Taking computation efficiency and inversion feasibility into consideration, we used a quadtree sampling approach (Jónsson et al., 2002) to produce a manageable data set with appropriate spatial resolution.Meanwhile, the look vectors of each observation are calculated considering the differences of incidence and azimuth angle based on the precise orbit data and local topography.The final InSAR data set consists of 2676 observations, including 1190 points from the ascending strip-mode path and 1486 points from the descending wide-swath path.
The fault geometry is based on the geologically mapped surface rupture (Sun et al., 2012;Guo et al., 2012), aftershock distribution and ascending and descending InSAR coseismic deformation patterns.Four rectangular fault segments are used to approximate the whole rupture fault, whose strike and length are determined with the information of the surface ruptures locations (Table 2).While the surface trace of the fault is defined based on field investigation, aftershock distribution and InSAR data, the dip angles of fault segments are difficult to determine.In our study, a grid search between the root mean square (RMS) of model residuals and the fault dip with wide range is used to infer the dip angle of each fault segment.In searching for best fitting fault dips, the dip angle of each fault segment is varied at 1 • interval while fixing other parameters.The best fitting model (Fig. 4) shows that the dip angles of F3 and F4 segments are 80 • dipping to north and 85 • dipping to south, respectively and well constrained, while the dip angles of F1 and F2 segments are worse constrained and vary in a wide range (±10 • ).This result is in general agreement with the dip angle of 83 • dipping to north from aftershocks relocation using data from both movable and fixed digital seismic monitoring networks (Q.X. Liu et al., 2012).To determine to what width (depth) slip could be allowed to occur, a series of inversions were carried out for a width range between 10 km and 32 km (e.g.Wen et al., 2012).Figure 3a shows the plot of RMS against inversion width.It is clear that when the allowed width is shallower than ∼20 km, the RMS is large, indicating that the inversion is poor at matching the observations.The maximum width of the fault is fixed to 20 km.After several iterations, the fault geometry is fixed for final slip inversion (Table 2).
It is well known that kinematic slip inversion requires regularization to avoid oscillating slip distribution.The most common regularization approach for slip inversion is to minimize the RMS of residuals and smoothing in the slip distribution (e.g.Li et al., 2011;Tobita et al., 2011).Unlike the regularization methods applied in other inversions, a regularization to minimize the RMS and roughness of stress drop distribution is introduced in this study.The static (shear) stress drop related with earthquakes is a fundamental parameter, which linearly depends on slip (Scholz, 2002).A constant stress drop results normally in a slip pattern like a parabolic or elliptic distribution that is from the physical point of view a more realistic approximation than a uniform slip model.In addition, the slip smoothing constraint usually causes stress drop singularities at fault edges or at connections between neighboring fault segments, which can be overcome using the stress drop smoothing constraint.Using the previously determined fault geometry, the fault plane is divided into 440 sub-faults with a size of ∼2 km by 2 km.Then a constrained least square method (Wang et al., 2009) is used to solve the objection function: where G is the Green function relating unit model slip to the predicted displacement, m is the coseismic slip on each patch, d is InSAR LOS displacement, H is the second-order  Laplacian smoothing operator, α is the positive smoothing factor that can be derived from trade-off between model roughness and data misfit (Fig. 5b), and τ is shear stress drop.The stress drop is defined by τ = G τ m, where G τ is the shear stress kernel matrix linearly related to the slip distribution m on the whole fault plane.And in the inversion, a layered earth model (Table 3) from the receiver function inversion (Wang et al., 2012) is employed for the Green function calculation.
Figure 6 shows the best fitting coseismic slip distribution of the Yushu earthquake.The calculated slip distribution images three areas of high slip (asperities) greater than 0.5 m.The inverted rakes show that the fault of upper 12 km is dominated by left-lateral strike-slip motion.The peak slip of 2.0 m occurs on the Jiegu segment (F3) at a depth of 3 km.The second largest slip (up to 1.5 m) patch occurs on the F2 segment, with motion constrained at depth of 3-14 km.The Longbao Lake segment (F1) has over 1.1 m slip with motion   constrained to the upper 12 km.The Yushu segment (F4) is observed with the maximum slip of 1.3 m and significant slip down to 7 km from the surface.
To estimate the uncertainties in the slip distribution from atmospherically correlated noise in the interferograms, we perturb the original interferogram data set with random synthetic correlated noise to perform a Monte Carlo error analysis (Hanssen, 2001;Wright et al., 2003).The inversion results from 100 data sets perturbed with noise of the statistical properties from previous 1-D covariance function are used to estimate the standard deviation in slip, and the resulting uncertainties are shown in Fig. 6d.The large magnitude of the surface deformation makes the impact of atmospheric noise small on the confidence in the slip value (the average errors of 5 cm), except for right-bottom parts of the segment F4 where the error increases to 10 cm.
The seismic moments released on segments F1, F2, F3 and F4 are 7.7 × 10 18 N m, 4.4 × 10 18 N m, 7.6 × 10 18 N m and 3.0 × 10 18 N m, respectively.The total seismic moment is 2.27 × 10 19 N m (M w = 6.9), which is lower than the moment of 3.7 × 10 19 N m derived from seismic wave modeling (Zhang et al., 2010) and slightly higher than the moment of 2.2 × 10 19 N m from a joint inversion of Envisat and ALOS SAR images (Li et al., 2011).Figure 3c-f shows the predicted and residual interferograms from our best fitting slip model.It is clear that the distributed slip model can explain most patterns of both descending and ascending observations with a small RMS misfit: 1.8 cm for 487A, and 2.5 cm for 139D interferogram, and those numbers are similar to the level of atmospheric noise in the InSAR coseismic displacements.The correlation coefficient between the observations and predictions is 98.5 %.However, there are some significant mismatches between measured and modeled deformation in the left-top and right-bottom regions of 139D path, which can be attributed to phase unwrapping errors, DEM errors, atmospheric propagation errors and other unknown errors.Moreover, the observed and modeled LOS displacements for profiles perpendicular to the strike for the four segments are shown in Fig. 7.There is good agreement between the observed LOS displacements and the modeled data, and the residuals are at a few cm level.

Discussion
We compare other available coseismic slip models (Table 4) with our best fitting model.In our model, the length of significant near-surface slip and the maximum of the segment F1 is 15 km and 0.6 m, respectively; and segments F3 and F4 have 35 km length of significant near-surface slip with a maximum slip of 2.0 m.The geological field investigation (Sun et al., 2012) suggests that the Jielong section (F1) ruptured a length of ∼15 km with a maximum left-lateral 0.66 m and the rupture length of Jiegu section (F3 and F4) is about 31 km with a maximum left-lateral slip of 2.4 m.The field observations of surface rupture (Sun et al., 2012) agree well with the near-surface slip from our model (Fig. 6a) both in distribution and magnitude.Teleseismic P-waves inversion (Zhang et al., 2010) indicated that the slip concentrated in the 10-40 km section near the hypocenter with a maximum slip of 2.1 m.Our slip shows a good agreement with the geological (Sun et al., 2012) and the seismic studies (Zhang et al., 2010) in both the length of rupture and the magnitude of maximum slip, but it slightly differs from the 1.25 m maximum near-surface and 1.5 m maximum slip from the ascending and descending ASAR and ascending PALSAR data (Li et al., 2011) and 1.7 m maximum near-surface and 2.6 m maximum slip at a depth of 7 km from the ascending and descending PALSAR data (Tobita et al., 2011).Because there is no obvious surface rupture occurred on the segment F2, the movement characteristic of this segment is very difficult to determined from field work (e.g.Sun et al., 2012), and it is also not well resolved in teleseismic inversion (Zhang et al., 2010).The slip model determined from InSAR observations shows that there are obvious left-lateral motion occurring at depth of 3-14 km with a maximum slip of 1.5 m.Moreover, our best fitting slip model shows that segment F1, F2, F3 and NW part of F4 have significant slips, but the slip on SE part of F4 is relatively small and therefore indicates slip deficit remains.
The Ganzi-Yushu fault is one of the major left-lateral active faults in eastern Tibet.Historic ruptures on the fault in the last 200 yr include the 1738 Yushu and its northwestern M = 6.5 earthquake, the 1886 Cuoa-Ganzi M = 7.3 earthquake and the 1896 Shiqu-Yushu M = 7.5 earthquake (Chen et al., 2010).Assuming that the 2010 event released all the accumulated strain energy, an recurrence interval of about 253-298 yr can be obtained from the interseismic slip rate of 7.3 ± 0.6 mm a −1 (Zhou et al., 1996) and the maximum coseismic slip of 2.0 m, and this means that the 2010 Yushu earthquake is a repeat on the location of the 1738 Yushu and its northwestern M = 6.5 earthquake.
From a statistical study on stress drops of 2000 earthquakes of mb ≥ 5.5 between 1990 and 2007, Allmann and Shearer (2009) indicate that stress drop estimates for individual earthquakes range from about 0.3 to 50 MPa, but the median stress drop does not vary with moment, and intraplate earthquakes stress drops have a median of 5.95 MPa.The average stress drops of segments F1, F2, F3 and F4 are 1.2, 2.7, 2.0 and 0.7 MPa, respectively, which fall within the typical range of earthquake stress drop by Allmann and Shearer (2009) and Shaw (2009).And the maximum stress drops of each segment are 4.2, 7.3, 8.3 and 5.8 MPa, respectively, which are consist with the other studies (e.g.Sun et al., 2011;Elliott et al., 2012).
In order to evaluate the static stress changes in the surrounding regions affected by the Yushu earthquake, we calculate the Coulomb stress changes (King et al., 1994) where τ and σ n are changes in shear and normal stress on a given fault plane, µ = 0.75 is the friction coefficient and B is the Skempton's coefficient, which has been measured to range between 0.47 and 1.0 (e.g.Parson et al., 1999;Serpelloni et al., 2012).In this study, the Coulomb stress changes are calculated with Skempton's coefficient between 0.47 and 1.0 and depth from 5 to 20 km, then projected on receiver fault with 110 • strike, 90 • dip and −6 • rake, which is the geometry of the mainshock rupture (F1 segment).The calculated Coulomb stress changes show that there was sensitivity in the results to parameter choices, but not to the degree that altered the sign of the calculated stress changes at specific target fault zones, or changed the general pattern of stress change (Fig. 8).The overall stress change pattern consists with the stress changes using a single fault model (Shan et al., 2011).Assuming the same receiver fault geometry for the aftershocks as that specified in the stress calculation, the correlation between M > 2.0 aftershocks until 16 December 2011 and static Coulomb stress changes shows that 67-89 % aftershocks occurred in the zones of increased stress.The aftershocks mainly occurred in the NW part of Longbao Lake segment (F1) with increased stress.However, the SE part of Yushu segment (F4) with increased stress only has little aftershocks, together with the slip deficit in our best fitting slip model; this area has higher seismic hazard possibility.

Conclusions
In this study, we map the coseismic deformation of the 2010 Yushu earthquake with ALOS PALSAR ascending stripmode and descending wide-swath data.The significant displacement is distributed over an area of 90 km by 70 km.The maximum observed LOS displacement in the descending interferogram is 0.37 m, and 0.50 m in the ascending interferogram.Inversion from a layered crustal model based on stress drop smoothing constraint suggests that there are three slip asperities identified.The slip occurred mainly in the upper 12 km, with a maximum slip of 2.0 m at depth of 3 km on the Jiegu segment.The inverted geodetic moment is approximately 2.27 × 10 19 N m (M w = 6.9), which is consistent with seismology results.Moreover, the evidences from the slip distribution derived from the InSAR observations, the calculated static Coulomb stress changes and the distribution of aftershocks in the SE region along Yushu segment (F4), indicate that a continuing seismic hazard exists here.

2 8. 5 *
Columns show track numbers with satellite directions (Ascending and Descending), incidence angle, azimuth angle and perpendicular baseline B ⊥ in the center of the scene.The standard deviation (σ ) and e-folding length scale (Distance) are those calculated for 1-D covariance function derived from the interferogram noise and are used in Monte Carlo estimation of fault parameters errors.

Fig. 2 .
Fig. 2. The ascending track 487 strip mode (a) and descending track 139 wide-swath mode (b) raw interferograms, corrections from orbital errors (c and d), and corrections from a linear tropospheric delay related to topography (e and f).The data in dashed region are used for correction parameters estimation.

Fig. 3 .
Fig. 3. Ascending track 487 strip mode interferogram after trend removal and rewrapped to 6 cm fringes (a), simulated interferogram (c) and residual (e) based upon the four fault segments (red lines) slip distribution.Descending track 139 wide-swath mode interferogram after trend removal and rewrapped to 6 cm fringes (b), simulated interferogram (d) and residual (f) based upon the four fault segments slip distribution.The black lines AA , BB , CC and DD are profiles used in Fig. 7.

Fig. 4 .
Fig. 4. Plots of the RMS misfit with various fault dip angles.(a) F1 segment, (b) F2 segment, (c) F3 segment and (d) F4 segment.The dip angle of each fault segment is varied at 1 • interval while fixing other parameters in searching for best fitting fault dips.The text in upper-right corner indicates the dip direction of each fault segment.

Fig. 5 .
Fig. 5. Plots of data misfit and inversion width and model roughness for a range of smoothing factor.(a) Trade-off curve between the data misfit and the maximum slip width.(b) Trade-off curve between the misfit and model roughness with a width extent of 20 km.The roughness is the normalized value.The pluses represent the values chosen for the inversions presented in Fig. 6.

Fig. 7 .
Fig. 7. Profiles of LOS displacements (blue dots), simulated LOS displacements (red dots) and topography for cross sections taken perpendicular to the four fault segments.(a, c, e and g) are profiles on ascending track 487.(b, d, f, and h) are profiles on descending track 139.

Table 1 .
Details of SAR data used in the study * .

Table 2 .
Modeled fault parameters of Yushu earthquake.

Table 3 .
Layered crustal structure parameters of Yushu fault region.

Table 4 .
Fault rupture parameters of the 2010 Yushu earthquake from different studies.