Joint Inversion of GPS, Leveling, and InSAR Data for The 2013 Lushan (China) Earthquake and Its Seismic Hazard Implications

: On 20 April 2013, a moment magnitude (Mw) 6.6 earthquake occurred in the Lushan region of southwestern China and caused more than 190 fatalities. In this study, we use geodetic data from nearly 30 continuously operating global positioning system (GPS) stations, two periods of leveling data, and interferometric synthetic aperture radar (InSAR) observations to image the coseismic deformation of the Lushan earthquake. By using the Helmert variance component estimation method, a joint inversion is performed to estimate source parameters by using these GPS, leveling, and InSAR data sets. The results indicate that the 2013 Lushan earthquake occurred on a blind thrust fault. The event was dominated by thrust faulting with a minor left-lateral strike–slip component. The dip angle of the seismogenic fault was approximately 45.0 ◦ , and the fault strike was 208 ◦ , which is similar to the strike of the southern Longmenshan fault. Our ﬁnite fault model reveals that the peak slip of 0.71 m occurred at a depth of ~12 km, with substantial slip at depths of 6–20 km. The estimated magnitude was approximately Mw 6.6, consistent with seismological results. Furthermore, the calculated static Coulomb stress changes indicate that the 2013 Lushan earthquake may have been statically triggered by the 2008 Wenchuan earthquake. rake angle of 90 ◦ . The checkerboard tests at the shallowest and middle subfaults were conducted to conﬁrm high resolution at the surface to middle depth (Figure 8a,b,e,f). The tests also showed that the most of the slip patches were smeared in the retrieved slip distribution mode for the complex conﬁgurations (Figure 8i-p) while large asperity patch could be recovered easily in the resulting slip models (Figure 8a–h), which is similar to the result of Huang et al. [20]. Furthermore, the two new CGPS stations of QLA1 and MINS at the nearﬁeld of the seismogenic fault were also tested to ﬁnd whether they will a ﬀ ect the result of the slip distribution. We found that the QLA1 and MINS stations increased ~10% maximum slip (Figure 8d,h,l,p), which proved that two CGPS stations could provide more constraint to the slip distribution. These results suggest that the slip asperity of four or less patches can be well determined from the joint inversion of GPS, leveling, and InSAR observations.


Introduction
On 20 April 2013, an earthquake with a moment magnitude (Mw) of 6.6 struck Lushan County, Ya'an city, Sichuan Province, southwestern China, generating considerable economic losses in the local region [1,2]. The event was the result of a thrust earthquake that occurred in this area following the 2008 Mw 7.9 Wenchuan earthquake. The focal mechanism solution released by the United States Geological Survey (USGS) revealed a focal depth of approximately 12.3 km and a fault with a dip of 33 • [3]. After this event, both near-field and far-field broadband digital seismic data were used to constrain the seismic source rupture process and focal mechanism in combination with other relevant information [4][5][6][7][8][9]. For example, Zhao et al. [5] used waveform data from the Global Seismographic Network (GSN)/Integrated Research Institutions for Seismology (IRIS), China Earthquake Networks Center (CENC) and Sichuan regional seismic network to perform a focal mechanism inversion; their findings indicated that the mainshock was a thrust earthquake that occurred on a fault with a dip of 45 • . Liu et al. [6] reported that the main rupture asperity was located around the hypocenter with a range of~28 km, and the peak slip was approximately 1.5 m. Considering the rupture source process, Wang et al. [9] reported that the dip was 38.5 • , and the maximum slip occurred at depths of 10-12 km.
However, no obvious surface rupture was observed in a field geological survey, and thus the Lushan earthquake was classified as a typical blind reverse-fault earthquake [10].
With an increasing amount of data, a variety of researchers have carried out a large number of inversions of the coseismic deformation associated with the 2013 event using global positioning system (GPS), leveling, interferometric synthetic aperture radar (InSAR), and strong motion and gravity data either separately or jointly [11][12][13][14][15][16][17][18][19][20]. For example, Jiang et al. [11] used observations from 33 continuously operating GPS (CGPS) stations to invert the coseismic slip distribution; their results suggested that the rupture was dominated by thrust faulting with a slight but still statistically significant left-lateral strike-slip component. Liu et al. [13] used InSAR data to image the coseismic deformation and calculated a further distributed slip model. Hao et al. [12] employed two campaigns of precise leveling observations across the southern section of the Longmenshan fault (LMSF) zone to quantify the coseismic vertical deformation and suggested that the seismogenic fault was dominated by a thrust component. Tan et al. [16] combined near-field coseismic offsets acquired both by continuously operating and by campaign GPS stations as well as seismic strong motion stations to reveal the rupture model and analyze its tectonic implications. Zhang et al. [18] jointly inverted local strong motion data and geodetic (including GPS and InSAR) measurements to obtain a robust model of the rupture process. Most recently, Huang et al. [20] performed a joint inversion with leveling, GPS, InSAR, and strong motion data to constrain the fault geometry and slip distribution.
Nevertheless, due to differences among the data sets, the findings of the previous studies all contrasted considerably with regard to the details of the fault geometry and coseismic slip distribution. The main difference among these results is the dip angle of the rupture plane, which varies by about 8 • . This discrepancy may be attributed to the low spatial resolution of far-field seismic waveforms and an insufficient amount of near-field geodetic data. As a result, previous studies generally assumed that the rupture plane was rectangular with a single dip angle; thus, the geometry of the rupture surface at depth could not be finely constrained.
Strain accumulation and stress change due to the mainshock have also been rigorously discussed in previous researches [21][22][23][24]. For instance, Shan et al. [21] calculated the co-and post-seismic stress changes caused by the 2008 Mw 7.9 Wenchuan and 2013 Ms 7.0 Lushan earthquakes to investigate their relationship and the stress changes on the major faults in this region. Zhang et al. [18] found a Coulomb stress change of approximately 0.5 MPa that occurred on the secondary fault segments simultaneously with the initiation of coseismic slip on the main fault of the 2013 event, indicating that the slip was likely triggered by coseismic slip on the main blind thrust fault. In addition, Lin et al. [24] investigated the Coulomb stress changes during the 2017 Mw 6.5 Jiuzhaigou and 2013 Mw 6.6 Lushan earthquakes and discussed their relationships to the 2008 Wenchuan earthquake.
In this paper, in addition to the CGPS data as used in Jiang et al. [11], we collected measurement data from six CGPS stations (two of them situated atop the footwall of the seismogenic fault responsible for the 2013 Lushan earthquake); these GPS data will provide better constraints for the slip distribution inversion along the rupture plane. We first acquired CGPS observed coseismic deformation around the seismogenic fault area of the Lushan earthquake, as well as InSAR data and leveling data. Then, we adopted a two-step method to retrieve the source parameters associated with this event to determine the geometry of the seismogenic fault and to estimate the slip distribution along the fault plane using these GPS, leveling, and InSAR data. Finally, we calculated the Coulomb stress change based on our preferred slip model and discussed the relationship to the Wenchuan earthquake.

Tectonic Background
The Lushan earthquake occurred in a region characterized by complex tectonism to the immediate north of an area that can be referred to as a triple junction, with the Songpan-Ganzi block to the northwest, the Chuandian block to the southwest and the Sichuan Basin to the east [18] (Figure 1). The epicenter of the Lushan earthquake is located in the southern segment of the LMSF zone, which also hosted the 2008 Wenchuan Mw 7.9 earthquake to the northeast. The LMSF constitutes the boundary Remote Sens. 2020, 12, 715 3 of 18 between the Songpan-Ganzi block and the Sichuan Basin; this boundary is characterized by crustal shortening caused by the southeastward movement of the Songpan-Ganzi block and collision with the South China Block [10]. The Songpan-Ganzi block is bounded by the Xianshuihe fault (XSHF) to the west, while to the south, the Anninghe fault (ANHF) separates the Chuandian block from the Sichuan Basin ( Figure 1). Both the XSHF and the ANHF are large and active strike-slip faults with dextral strike-slip rates of >10 mm/yr and 3-8 mm/yr, respectively [25]. In contrast, the LMSF experiences movement of less than 3 mm/yr with both thrust and dextral strike-slip components, as determined from the surface geology and GPS measurements [26,27].
Through a field geological survey, the Lushan earthquake was classified as a typical blind reverse-fault earthquake [10]. Li et al. [28] found that the Lushan event was initiated along a detachment fault underneath the front range of the LMSF between the Shuangshi-Dachuan fault (SDF) and the Xinkaidian fault (XF) based on observations of surface cracks. Xu et al. [10] found several northeast-striking imbricate reverse faults in the southern segment of the LMSF, but no apparent earthquake surface rupture zones were located along these active faults or their adjacent areas. Based on an analysis of precisely relocated aftershocks, most aftershocks were confined to depths from 5 to 20 km, while few earthquakes occurred in the shallow crust [29].
Remote Sens. 2020, 12,715 3 of 20 boundary between the Songpan-Ganzi block and the Sichuan Basin; this boundary is characterized by crustal shortening caused by the southeastward movement of the Songpan-Ganzi block and collision with the South China Block [10]. The Songpan-Ganzi block is bounded by the Xianshuihe fault (XSHF) to the west, while to the south, the Anninghe fault (ANHF) separates the Chuandian block from the Sichuan Basin ( Figure 1). Both the XSHF and the ANHF are large and active strikeslip faults with dextral strike-slip rates of >10 mm/yr and 3-8 mm/yr, respectively [25]. In contrast, the LMSF experiences movement of less than 3 mm/yr with both thrust and dextral strike-slip components, as determined from the surface geology and GPS measurements [26,27]. Through a field geological survey, the Lushan earthquake was classified as a typical blind reverse-fault earthquake [10]. Li et al. [28] found that the Lushan event was initiated along a detachment fault underneath the front range of the LMSF between the Shuangshi-Dachuan fault (SDF) and the Xinkaidian fault (XF) based on observations of surface cracks. Xu et al. [10] found several northeast-striking imbricate reverse faults in the southern segment of the LMSF, but no apparent earthquake surface rupture zones were located along these active faults or their adjacent areas. Based on an analysis of precisely relocated aftershocks, most aftershocks were confined to depths from 5 to 20 km, while few earthquakes occurred in the shallow crust [29].

GPS Data
We collected nearly two years (2013-2014) raw GPS data from near-field GPS stations surrounding the epicenter of the Lushan earthquake managed by different agencies, including 10 permanent GPS sites from the Crustal Movement Observation Network of China, 10 regional CGPS stations deployed by Dr. Jiang from the China Earthquake Administration, two CGPS stations (YAAN and SCQL) from the Sichuan Provincial Seismological Bureau, and six CGPS stations constructed by the Sichuan Bureau of Surveying, Mapping and Geoinformation (SBSMG) (Figure 1, blue squares). Two sites (QLA1 and MINS) from the SBSMG were atop the footwall of the seismogenic fault, which we believe will provide better constraints on the fault distribution than those offered in previous researches [11,18,20].
We used BERNESE 5.2 software [32] to process the raw data, including the implementation of daily solutions and adjustments. The data processing of a daily solution mainly includes data preprocessing, baseline combination, eliminating integer ambiguities, and calculating free network solutions. Pseudorange and phase processing modes were adopted for the integer ambiguity elimination strategy with a quasi-ionosphere-free (QIF) wide-lane (WL) ambiguity resolution strategy based on the pseudorange and an L5 strategy based on the phase. To maintain the regional framework throughout China and establish a connection with the global framework during the daily solutions, six International Global Navigation Satellite System (GNSS) Service (IGS) stations throughout China and its adjacent areas, namely, BJFS, DAEJ, IISC, POL2, TCMS, and ULAB, were used to acquire joint daily solutions. During the daily network adjustment process, Solution Independent Exchange (SINEX) files of global daily solutions from the CODE Data Processing and Analysis Center were downloaded and converted into normal equation files, which were then superimposed onto the aforementioned normal equation files of daily solutions. By fixing the aforementioned six IGS stations and global framework stations, the regional daily solutions and global framework daily solutions were integrated. With the baseline constraint, all stations under the IGB08 framework were used as baseline framework points. Additionally, a priori unit weight mean error of all the reference points were set as 0.001 m for the baseline constraint [33][34][35]. For the least-square constraint, the global seven-parameter transformation method was adopted instead of the block adjustment three-parameter transformation method [36][37][38]. Finally, single-day net solution time series with high stability and high reliability were obtained. Figure 2 illustrates examples of the high-precision time series of LS05 and SCTQ, from which coseismic deformation can be clearly observed.
After acquiring high-precision single-day measurements, we calculated the coseismic deformation produced by the Lushan earthquake with a five-day daily position time series before and after the earthquake. By averaging the preseismic and postseismic positions and subtracting the latter from the former, we obtained the coseismic deformation in the region. Here, we subtracted the coseismic deformation and preseismic deformation from the postseismic deformation to define a new time series, and we propose using the standard deviation of this new time series to evaluate the coseismic deformation error. The weighted root mean squared (wrms) is as follows: where y i stands for the site position before or after the earthquake from the new time series, y describes its average, and n is the number of GPS observations in the new time series. The corresponding coseismic deformation error, rms, is calculated on the error propagation law as, rms = wrms 2 pre + wrms 2 post . The coseismic deformation from these CGPS stations are documented in Table 1 and Figure 3a. The largest coseismic displacement was approximately 68.0 mm horizontally and 82.6 mm vertically at site LS05, which are in good agreement with the results from Jiang et al. [11]. After acquiring high-precision single-day measurements, we calculated the coseismic deformation produced by the Lushan earthquake with a five-day daily position time series before and after the earthquake. By averaging the preseismic and postseismic positions and subtracting the latter from the former, we obtained the coseismic deformation in the region. Here, we subtracted the coseismic deformation and preseismic deformation from the postseismic deformation to define a new time series, and we propose using the standard deviation of this new time series to evaluate the coseismic deformation error. The weighted root mean squared (wrms) is as follows: where stands for the site position before or after the earthquake from the new time series, describes its average, and n is the number of GPS observations in the new time series. The corresponding coseismic deformation error, rms, is calculated on the error propagation law as, = + . The coseismic deformation from these CGPS stations are documented in Table   1 and Figure 3a. The largest coseismic displacement was approximately 68.0 mm horizontally and 82.6 mm vertically at site LS05, which are in good agreement with the results from Jiang et al. [11].   Note: E, N, and U represent the coseismic deformation of east, north, and up direction, respectively. Se, Sn, and Su represent the coseismic deformation error of east, north, and up direction, respectively.

Leveling Data
The leveling data is from the work of Hao et al. [12]. Here we just give a simple illustration of the leveling observation. A historical leveling line is located in the southwestern part of the study area across the LMSF starting at Jiajiang, after which it passes through Ya'an and Lushan and ends at Xiaojin in the southeast (Figure 1). With a length of approximately 294 km, this leveling line traverses the Xinkaidian fault (XF) and the Shuangshi-Dachuan (SDF) ruptures. The China State Bureau of Surveying and Mapping took first-order measurements from Jiajiang to Ya'an in October 1985 and from Ya'an to Xiaojin in July 1986. After the 2008 Wenchuan earthquake, the SBSMG carried out second-order leveling remeasurements along this route from April to November 2010 and embedded the benchmarks that ruptured during the earthquake. Based on these two periods of leveling measurement data, the vertical crustal deformation across the southern segment of the LMSF between 1985/1986 and 2008 was obtained. In 2013, the Second Crustal Deformation Monitoring Center, China Earthquake Administration, acquired first-order precise leveling measurements along the Ya'an-Baoxing route within three months after the Lushan earthquake struck Sichuan Province. Similarly, the vertical deformation was calculated along this route based on the two leveling data sets from 2010 and 2013. The calculated deformation also included the interseismic crustal deformation from 2010 to 2013 (prior to the occurrence of the Lushan earthquake) and the effects of afterslip and viscoelastic relaxation following the 2008 Wenchuan earthquake [12].
The following methods were employed by Hao et al. [12] to remove the impacts of crustal deformation before the Lushan earthquake first, they computed the interseismic vertical deformation rate by subdividing the vertical deformation during 1985-2010 by the time interval; next, they determined the corresponding vertical deformation from 2010 to 2013 and subtracted it from the two leveling data sets described above. Additionally, they adopted a cubic spline function to interpolate the interseismic vertical deformation rate of the newly added benchmarks during 2010 and 2013. By these means, they finally obtained the coseismic deformation of the Lushan earthquake according to leveling observations (Table 2 and Figure 4a). Based on theoretical estimates and analyses of the vertical deformation results in the cases of some earthquakes, including the 2008 Wenchuan

Leveling Data
The leveling data is from the work of Hao et al. [12]. Here we just give a simple illustration of the leveling observation. A historical leveling line is located in the southwestern part of the study area across the LMSF starting at Jiajiang, after which it passes through Ya'an and Lushan and ends at Xiaojin in the southeast (Figure 1). With a length of approximately 294 km, this leveling line traverses the Xinkaidian fault (XF) and the Shuangshi-Dachuan (SDF) ruptures. The China State Bureau of Surveying and Mapping took first-order measurements from Jiajiang to Ya'an in October 1985 and from Ya'an to Xiaojin in July 1986. After the 2008 Wenchuan earthquake, the SBSMG carried out second-order leveling remeasurements along this route from April to November 2010 and embedded the benchmarks that ruptured during the earthquake. Based on these two periods of leveling measurement data, the vertical crustal deformation across the southern segment of the LMSF between 1985/1986 and 2008 was obtained. In 2013, the Second Crustal Deformation Monitoring Center, China Earthquake Administration, acquired first-order precise leveling measurements along the Ya'an-Baoxing route within three months after the Lushan earthquake struck Sichuan Province. Similarly, the vertical deformation was calculated along this route based on the two leveling data sets from 2010 and 2013. The calculated deformation also included the interseismic crustal deformation from 2010 to 2013 (prior to the occurrence of the Lushan earthquake) and the effects of afterslip and viscoelastic relaxation following the 2008 Wenchuan earthquake [12].
The following methods were employed by Hao et al. [12] to remove the impacts of crustal deformation before the Lushan earthquake first, they computed the interseismic vertical deformation rate by subdividing the vertical deformation during 1985-2010 by the time interval; next, they determined the corresponding vertical deformation from 2010 to 2013 and subtracted it from the two leveling data sets described above. Additionally, they adopted a cubic spline function to interpolate the interseismic vertical deformation rate of the newly added benchmarks during 2010 and 2013. By these means, they finally obtained the coseismic deformation of the Lushan earthquake according to leveling observations ( Table 2 and Figure 4a). Based on theoretical estimates and analyses of the vertical deformation results in the cases of some earthquakes, including the 2008 Wenchuan earthquake, the postseismic vertical deformation caused by afterslip and viscoelastic relaxation in the lower crust was determined to range from approximately several millimeters to centimeters. Thus, these effects were disregarded when the coseismic deformation of the Lushan earthquake was calculated. For more details about instruments and data process strategy for the leveling observations, please refer to the work of Hao et al. [12]. calculated. For more details about instruments and data process strategy for the leveling observations, please refer to the work of Hao et al. [12]. Note: U and represent the vertical coseismic deformation and its error, respectively.

InSAR Data
We employed Canada RADARSAT-2 single look complex (SLC) data to obtain the coseismic deformation of the Lushan earthquake along the satellite line of sight (LOS). Considering the

InSAR Data
We employed Canada RADARSAT-2 single look complex (SLC) data to obtain the coseismic deformation of the Lushan earthquake along the satellite line of sight (LOS). Considering the coverage area and temporal-spatial baseline, we selected a pair of ascending SAR images, acquired on 10 January 2013 and 10 May 2013 with a perpendicular baseline of 186.3 m, to acquire the interferogram. The SAR data were processed by Switzerland GAMMA software [39] in two-pass differential interferometry mode [40]. The interferogram was corrected for the effects of topography using a 90-m-resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) and MDA-provided orbit data. The interferograms were filtered using a power spectrum filter to reduce the effects of phase noise, unwrapped with the branch-cut method [41,42], and finally geocoded into the World Geodetic System 1984 (WGS84) coordinate system.
The interferograms resolve coseismic displacements only to the east of the Lushan epicenter, primarily in the Sichuan Basin and on the footwall of the thrust fault (Figure 5a). Large LOS signals were observed far from the fault in the Sichuan Basin, where an earthquake of this magnitude would not produce such significant coseismic ground displacements. These signals were likely caused by factors such as atmospheric turbulence. Excluding these errors, the InSAR-derived displacements qualitatively agree well with those observed by the GPS measurements. The standard deviation and decay distance estimated from a 1D covariance function [43] for the interferogram were 1.2 cm and 25.6 km, respectively.
Remote Sens. 2020, 12, 715 9 of 20 coverage area and temporal-spatial baseline, we selected a pair of ascending SAR images, acquired on 10 January 2013 and 10 May 2013 with a perpendicular baseline of 186.3 m, to acquire the interferogram. The SAR data were processed by Switzerland GAMMA software [39] in two-pass differential interferometry mode [40]. The interferogram was corrected for the effects of topography using a 90-m-resolution Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) and MDA-provided orbit data. The interferograms were filtered using a power spectrum filter to reduce the effects of phase noise, unwrapped with the branch-cut method [41,42], and finally geocoded into the World Geodetic System 1984 (WGS84) coordinate system. The interferograms resolve coseismic displacements only to the east of the Lushan epicenter, primarily in the Sichuan Basin and on the footwall of the thrust fault (Figure 5a). Large LOS signals were observed far from the fault in the Sichuan Basin, where an earthquake of this magnitude would not produce such significant coseismic ground displacements. These signals were likely caused by factors such as atmospheric turbulence. Excluding these errors, the InSAR-derived displacements qualitatively agree well with those observed by the GPS measurements. The standard deviation and decay distance estimated from a 1D covariance function [43] for the interferogram were 1.2 cm and 25.6 km, respectively.

Joint Inversion
Based on the coseismic deformation determined from the GPS, leveling, and InSAR measurements, we used an elastic half-space dislocation model [44] to retrieve the fault geometry and slip distribution of the Lushan earthquake. Here, a two-step method was adopted to estimate the source parameters associated with this event. The first step was to determine the geometry of the seismogenic fault by a nonlinear inversion method while assuming a uniform slip model composed of a rectangular fault plane. The second step was to estimate the slip distribution along the fault plane based on a linear inversion method.
The horizontal and vertical GPS coseismic displacements from 28 sites were adopted in the inversion. Their mean square errors in the east, north and vertical directions were 1.5 mm, 1.6 mm and 5 mm, respectively, as detailed in Table 1. Taking into account the computational efficiency and the feasibility of the inversion, a resolution-based downsampling method [45] was used to deduce the observations of the InSAR surface deformation field. Finally, 762 points with corresponding incidence and azimuth angles were derived from millions of observations.
The relative weighting ratio among the observed values in the joint inversion process was determined by the Helmert variance component estimation method [46] according to the corresponding mean square error . The relative weighting ratio of the GPS horizontal deformation: GPS vertical deformation: leveling: InSAR observation was finally determined to be 1.000: 0.052: 0.021: 0.010.

Joint Inversion
Based on the coseismic deformation determined from the GPS, leveling, and InSAR measurements, we used an elastic half-space dislocation model [44] to retrieve the fault geometry and slip distribution of the Lushan earthquake. Here, a two-step method was adopted to estimate the source parameters associated with this event. The first step was to determine the geometry of the seismogenic fault by a nonlinear inversion method while assuming a uniform slip model composed of a rectangular fault plane. The second step was to estimate the slip distribution along the fault plane based on a linear inversion method.
The horizontal and vertical GPS coseismic displacements from 28 sites were adopted in the inversion. Their mean square errors in the east, north and vertical directions were 1.5 mm, 1.6 mm and 5 mm, respectively, as detailed in Table 1. Taking into account the computational efficiency and the feasibility of the inversion, a resolution-based downsampling method [45] was used to deduce the observations of the InSAR surface deformation field. Finally, 762 points with corresponding incidence and azimuth angles were derived from millions of observations. The relative weighting ratio P i among the observed values in the joint inversion process was determined by the Helmert variance component estimation method [46] according to the corresponding mean square error σ i . The relative weighting ratio of the GPS horizontal deformation: GPS vertical deformation: leveling: InSAR observation was finally determined to be 1.000: 0.052: 0.021: 0.010.

Uniform Slip Inversion
During the nonlinear inversion, we employed a hybrid optimization algorithm based on multi-peak particle swarm optimization (M-PSO) [47] to retrieve the seismogenic fault parameters, including the fault position, strike, dip angle, length, depth, width, and slip. A Monte Carlo bootstrap method [48] was used to estimate the uncertainty of each fault parameter. Based on the statistical characteristics of the one-dimensional covariance function, the uncertainty of each fault parameter was estimated by using 100 combined observation data points with errors. Table 3 lists all the estimated fault parameters and their errors. The inversion results demonstrate that the errors were relatively small and that the seismogenic fault could be characterized by a blind fault with a dip of approximately 45 • that generated an Mw 6.56 earthquake.

Distributed Slip Inversion
With the fault geometry determined, the slip distribution along the fault plane can be inverted based on classic elastic dislocation theory. The optimal fault was extended to 36 km along the strike and 40 km along the dip, and the fault plane was divided into 1440 patches, each with a size of 1 km × 1 km. The constrained least-square method was used to solve the objective function: where G is Green's matrix, d represents the observation data, s is the slip of each patch, and L is the second-order finite difference approximation of the Laplace operator used to avoid a slip distribution characterized by nonphysical oscillations. κ 2 is a smoothing factor. In this study, we first analyzed the changing trends of the model roughness (ψ) and residuals (ξ) with the variation of smoothing weight [47]. To obtain the trade-off between the roughness (ψ) and residual (ξ), we built a Log-function f (κ) = log(ψ + ξ), which indicated a single minimum at κ 2 = 2.5 (grey solid line in Figure 6). Figure 7 shows the slip distributions and uncertainty based on the joint inversion of GPS, leveling, and InSAR observations. The results show that the seismogenic fault is mainly a thrust fault with a left-lateral strike-slip component. The total energy released was approximately 6.1 × 10 18 Nm (Mw 6.58), which is comparable to that estimated in the GCMT catalog but greater than that predicted by the USGS. The maximum fault slip was 0.71 m at a depth of approximately 13 km (21-25 km wide), and the average slip error was approximately 2-3 cm.            5 show the observed and simulated deformation and residuals for the GPS, leveling, and InSAR measurements, respectively, from the best-fitting slip model. Clearly, the coseismic deformation observed by the GPS, leveling, and InSAR observations can be effectively explained by the slip distribution model without any significant residual in the near field. The significant residuals of the interferogram are in the Sichuan Basin (Figure 5c), far from the epicenter, which is mainly attributed to the atmospheric effect. The standard deviations between the observed and simulated GPS, leveling, and InSAR values were 5.1 mm, 2.4 mm and 4.2 mm, respectively, from our joint inversion model.

Checkerboard Tests
Additionally, to assess the reliability of the obtained slip model, we performed checkerboard tests, for which we first construct four kinds of input models: a model with single-patch at the shallowest subfaults, a model with single-patch at the middle subfaults, a model comprising four patches, and a 10-patch classic model (Figure 8a,e,i,m). For these models, the slip of each patch was set to 1 m with a rake angle of 90 • . The checkerboard tests at the shallowest and middle subfaults were conducted to confirm high resolution at the surface to middle depth (Figure 8a,b,e,f). The tests also showed that the most of the slip patches were smeared in the retrieved slip distribution mode for the complex configurations (Figure 8i-p) while large asperity patch could be recovered easily in the resulting slip models (Figure 8a-h), which is similar to the result of Huang et al. [20]. Furthermore, the two new CGPS stations of QLA1 and MINS at the nearfield of the seismogenic fault were also tested to find whether they will affect the result of the slip distribution. We found that the QLA1 and MINS stations increased~10% maximum slip (Figure 8d,h,l,p), which proved that two CGPS stations could provide more constraint to the slip distribution. These results suggest that the slip asperity of four or less patches can be well determined from the joint inversion of GPS, leveling, and InSAR observations. tests, for which we first construct four kinds of input models: a model with single-patch at the shallowest subfaults, a model with single-patch at the middle subfaults, a model comprising four patches, and a 10-patch classic model (Figure 8a,e,i,m). For these models, the slip of each patch was set to 1 m with a rake angle of 90°. The checkerboard tests at the shallowest and middle subfaults were conducted to confirm high resolution at the surface to middle depth (Figure 8a,b,e,f). The tests also showed that the most of the slip patches were smeared in the retrieved slip distribution mode for the complex configurations (Figure 8i-p) while large asperity patch could be recovered easily in the resulting slip models (Figure 8a-h), which is similar to the result of Huang et al. [20]. Furthermore, the two new CGPS stations of QLA1 and MINS at the nearfield of the seismogenic fault were also tested to find whether they will affect the result of the slip distribution. We found that the QLA1 and MINS stations increased ~10% maximum slip (Figure 8d,h,l,p), which proved that two CGPS stations could provide more constraint to the slip distribution. These results suggest that the slip asperity of four or less patches can be well determined from the joint inversion of GPS, leveling, and InSAR observations.

Coseismic Stress Changes
Variations in regional stresses due to coseismic displacements are likely affected by the surrounding faults and seismic activity. To estimate the stress field, we used the static Coulomb stress change method by King et al. [49]. According to Coulomb's failure law, the static Coulomb failure stress ∆CFS on a given fault surface should obey the following formula: where ∆τ is the variation in the shear stress parallel to the fault surface, ∆σ is the variation in the stress vertically perpendicular to the fault plane (release is represented as a positive value), and µ is the friction coefficient. When the Coulomb stress increases, it will accelerate the rupture process on the fault plane; in contrast, when the Coulomb stress decreases, stress relief will occur, thereby suppressing fault rupture. It is worth noting that the Coulomb stresses that could trigger an earthquake are often smaller than tectonic stresses, and a Coulomb stress as low as 0.1 bar is sufficient to trigger an earthquake [49,50].
In this study, we used Coulomb v3.4 software [51], which was developed based on MATLAB, to calculate the variations in the static Coulomb stresses generated by the Lushan earthquake. The effective friction coefficient and shear elasticity modulus were set to be 0.4 and 3.32 × 10 10 N/m 2 , respectively [52]. The specific fault parameters in our region of investigation, that is, the strike, dip, and rake, were 208 • , 43 • , and 83 • (the same as the rupture fault), respectively, as determined from our inversion model ( Table 3). The calculated Coulomb stress distribution at a depth of 13 km (which is nearly the depth of maximum slip) is shown in Figure 9, revealing that the Coulomb stress increased in the epicentral area by more than 0.45 bar.  Next, we calculated the stress changes on nearby active faults (Figure 1) caused by the Wenchuan and Lushan earthquakes. The parameters of these active faults were adopted from Table 2 of Shan et al. [53]. The Coulomb stress changes along the surrounding faults caused by the Wenchuan and Lushan earthquakes (both separately and jointly) are shown in Figure 10. As a result of the Wenchuan earthquake, the Coulomb stress changes increased on the following faults: the middle part of the XSHF, the Wenchuan-Qingchuan fault (WQF), the Longriba fault (LRBF), the Beichuan-Yingxiu fault (BYF), and the southern segment of the Pengxian-Guanxian fault (PGF) (Figure 10a). As a result of the Lushan earthquake, the Coulomb stress changes increased on the following faults: the PGF and the southern segment of the BYF (Figure 10b). From these increases in the Coulomb stress changes, we can conclude that the Coulomb stress changes accumulated mainly along the PGF due to the Wenchuan earthquake, although we further deduce that the Lushan earthquake ultimately ruptured the PGF (Figure 10c).

Discussion
Using GPS, leveling, and InSAR data, the fault geometry and slip characteristics of the 2013 Lushan earthquake (Table 3) were retrieved based on an elastic half-space dislocation model. This investigation presents the full application of high-quality geodetic observations in the study area.

Discussion
Using GPS, leveling, and InSAR data, the fault geometry and slip characteristics of the 2013 Lushan earthquake (Table 3) were retrieved based on an elastic half-space dislocation model. This investigation presents the full application of high-quality geodetic observations in the study area. From the joint inversion of geodetic data, our model indicates that a one-segment fault is good for fitting the observed near-field deformation (Figures 3-5). The observed and modeled displacements from the GPS, InSAR, and leveling data were well fit, especially in the footwall of the seismogenic fault, due to the addition of observations from six CGPS stations. The findings agree well with those of Jiang et al. [11] and Huang et al. [20] but differ from those of Zhang et al. [18]. Jiang et al. [11] found that a one-segment fault geometry based on a layered half-space model [54] can be fit well with CGPS observations, while Zhang et al. [18] reported that a four-segment fault geometry based on a homogeneous elastic half-space model can be fit well when constrained by both CGPS and InSAR data. Huang et al. [20] used GPS, InSAR, leveling, and strong motion data to study different fault geometries based on a homogeneous elastic half-space model and a layered half-space model separately, based on the data residuals, and found that a one-segment fault model constrained by GPS and leveling data can provide a better fit in the near field than other fault models.
Because of the availability of more near-field CGPS observations from the footwall of the seismogenic fault, our slip distribution model precludes the presence of obvious deformation in the upper crust ( Figure 7) and demonstrates that the seismogenic fault is a blind thrust fault, which is in agreement both with the geological survey by Xu et al. [10] and with the study of Huang et al. [20] but different from the results from Jiang et al. [11] and Zhang et al. [18]. Huang et al. [20] verified that their slip distribution model yields no slip on the subfault patches near the surface, while the models developed by Jiang et al. [11] and Zhang et al. [18] presented coseismic slip on the near-surface subfaults. Our checkerboard tests also confirmed that there was no significant slip near the surface (Figure 8). Wen et al. [55] studied the asperities on the deep part of the seismogenic fault responsible for the 2015 Pishan earthquake in Xinjiang; they reported the presence of partial ruptures on the ramp fault at depths of 14-18 km that truncate the strata from the Paleozoic to the bottom of the Paleocene and transfer slip from the Proterozoic decollement beneath the area to the shallow Eocene layer beneath the Tarim Basin. Our research proves that the 2013 Lushan Mw 6.6 earthquake occurred along a blind thrust fault in a boundary zone between the Longmenshan fault and Sichuan Basin. These two moderate earthquakes, namely, the Pishan and Lushan events, illustrate that crustal material migrates outward from the deformed front within the Tibetan Plateau, causing lateral compression.
The centroid depth based on the uniform slip model was 7.4 km, and the maximum slip occurred at a depth of~12 km. These depths are shallower than those reported by the GCMT and USGS but are basically consistent with the GPS-derived results by Jiang et al. [11], with the main slip occurring at a depth of 13 km. One reason for this discrepancy is the variability of the local crust throughout a seismic array. In this paper, a homogeneous and elastic half-space model was adopted, but this model represents only a first-order approximation of the real crust. Moreover, in some cases, the focal depth of a layered crustal model is shallower than that of a half-space model [56]. A second possible reason is that the depth of the focal mechanism is the initial rupture depth, while the modeled depth is the centroid depth or peak slip depth of the fault. Based on over 80 finite-hypocenter rupture models of more than 50 large earthquakes, Mai et al. [57] concluded that shallow crustal ruptures tend to propagate upward from deeper layers along the fault plane. Another explanation for the shallow centroid depth of the model in this paper is the distribution of ground measurements. Near-field and far-field deformation are induced by seismic waves; in our model, near-field (0-100 km) observations, including InSAR, GPS, and leveling data, were used to provide strong constraints on the shallow slip distribution along the fault. The optimal slip model shows that the rupture of the 2013 Lushan earthquake was dominated by thrust faulting (with an 11% left-lateral strike-slip component), which is in agreement with two previous geodetic inversions [11,20] but disagrees with the findings of Zhang et al. [18]. The main reason for this may be that our model employs a one-segment fault geometry to perform the inversion, similar to the models of Jiang et al. [11] and Huang et al. [20], while Zhang et al. [18] [11,18,20]. The total energy released was approximately 6.1 × 10 18 Nm, which is close to that of Huang et al. (7.63 × 10 18 N m) but 30% smaller than that of Jiang et al. (9.5 × 10 18 N m), and 67% smaller than that of Zhang et al. (1.8 × 10 19 N m). Again, the differences among these studies may be affected by the different data sets. We utilized 2 CGPS stations on the footwall of the seismogenic fault, and these stations are very close to the epicenter; this may provide better constraints on the inversion than other studies, leading to a lower moment.
The changes in the ∆CFS calculated from the uniform slip model and the slip distribution model show that the Coulomb stress increased in the epicentral area by more than 0.45 bar (shown in Figure 9), which is larger (by more than 0.3 bar) than the increase in the epicentral area reported by Lin et al. [24] using teleseismic broadband waveforms. We depict the stress changes on nearby active faults caused by the Wenchuan and Lushan earthquakes in Figure 10. Following the Lushan earthquake, the CFS increased on the PGF and BYF by nearly 0.1 bar and 0.2 bar, respectively, which agree well with the results of Lin et al. [24]; it is possible that this stress change will result in future earthquakes. Following the Wenchuan earthquake, the CFS increased on the following faults (as shown in Figure 10a): the middle part of the XSHF, the WQF, the LRBF, the BYF and the southern segment of the PGF. Our results demonstrate that the Wenchuan earthquake increased the CFS by up to 0.2 bar at the hypocenter of the Lushan earthquake, which exceeded the threshold value (0.1 bar) of earthquake triggering [58]. Lin et al. [24] found that the Wenchuan earthquake had a significant impact on the surrounding area and verified that the Lushan earthquake was probably promoted (possibly triggered) by the 2008 Wenchuan earthquake. There were few aftershocks or smaller events about Mw 6 that occurred in this area between the 2008 Wenchuan and 2013 Lushan events. The coseismic static Coulomb stress changes induced by these small events would not significantly affect the stress distribution pattern in the studied area, such as the 2017 Mw 6.5 Jiuzhaigou earthquake [24].

Conclusions
In this study, the fault geometry and slip distribution of the 2013 Lushan earthquake were inverted by jointly using GPS, leveling, and InSAR data. We investigated the rupture parameters of the seismogenic fault by employing the Helmert variance component estimation method to determine the relative weighting ratio among the GPS, leveling, and InSAR observations. The inversion demonstrated that the 2013 Mw 6.6 Lushan earthquake occurred along a blind thrust fault, which was dominated by thrust faulting with a minor left-lateral strike-slip component. The dip angle of the seismogenic fault was approximately 45.0 • , and the fault strike was 208 • , which is generally consistent with the strike of the southern LMSF. The slip was concentrated mainly at depths of 6-20 km and the peak slip was 0.71 m, which occurred at a depth of approximately 12 km. The estimated moment magnitude was approximately Mw 6.6, which is consistent with the results from seismological studies but smaller than other geodetic inversions.
Finally, we calculated the static Coulomb stress changes triggered by the Lushan and Wenchuan earthquakes. The results show that the Lushan event increased the static Coulomb stress changes in the epicentral area by more than 0.45 bar. In addition, the Coulomb stress changed by the Wenchuan earthquake due to surrounding faults increased the stress by up to 0.2 bar at the hypocenter of the Lushan earthquake, confirming that the Lushan earthquake was probably promoted by the 2008 Wenchuan earthquake.
Author Contributions: Z.L. and Y.W. conceived and designed the experiments; Z.L. and P.Z. performed and analyzed the experiments; Z.L., Y.W., P.Z., and Y.L. wrote the paper; Y.Z. performed the data processing. All authors have read and agreed to the published version of the manuscript.