Source Model and Stress Disturbance of the 2017 Jiuzhaigou Mw 6.5 Earthquake Constrained by InSAR and GPS Measurements

Seismogenic fault geometry, especially for a blind fault, is usually difficult to derive, based only on the distribution of aftershocks and interference fringes of Interferometric Synthetic Aperture Radar (InSAR). To better constrain the fault geometry of the 2017 Jiuzhaigou Mw 6.5 earthquake, we first carried out a nonlinear inversion for a single fault source using multi-peak particle swarm optimization (MPSO), Monte Carlo (MC), and Markov Chain Monte Carlo (MCMC) algorithms, respectively, with constraints of InSAR data in multiple SAR viewing geometries. The fault geometry models retrieved with different methods were highly consistent and mutually verifiable, showing that a blind faulting with a strike of ~154◦ and a dip angle of ~77◦ was responsible for the Jiuzhaigou earthquake. Based on the optimal fault geometry model, the fault slip distribution jointly inverted from the InSAR and Global Positioning System (GPS) data by the steepest descent method (SDM) and the MC method showed that the slip was mainly concentrated at the depth of 1–15 km, and only one slip center appeared at the depth of 5–9 km with a maximum slip of about 1.06 m, some different from previous studies. Taking the shear modulus of μ = 32 GPa, the seismic moment derived from the distributed slip model was about 7.85 × 1018 Nm, equivalent to Mw 6.54, which was slightly larger than that from the focal mechanism solutions. The fault spatial geometry and slip distribution could be further validated with the spatial patterns of the immediate aftershocks. Most of the off-fault aftershocks with the magnitude > M2 within one year after the mainshock occurred in the stress positive stress change area, which coincided with the stress triggering theory. The static Coulomb stress, triggered by the mainshock, significantly increased at the Tazang fault (northwest to the epicenter), and at the hidden North Huya fault, and partial segments of the Minjiang fault (west of the epicenter).


InSAR Coseismic Deformation Measurements
The SAR data used in this study included the Single Look Complex (SLC) data in ascending and descending orbits from the Sentinel-1A/B and RADARSAT-2 satellite platforms (Table 1), operating at the C band with a radar wavelength of about 5.6 cm. The coseismic interferometry pairs of the Sentinel-1 ascending, descending, and RADARSAT-2 ascending orbits are referred to the S1A, S1D, and R2A in this study, respectively. A and D represent the ascending orbit and descending orbit, respectively. The superscript A and B represent the Sentinel-1A and Sentinel-1B satellites, respectively. ΔT is the time baseline with a unit of day. B⊥ is the vertical baseline with a unit in meters. ΔH2π is the height of ambiguity with a unit in meters. Abbr. is the abbreviation of corresponding interferogram pairs.
We used the InSAR Scientific Computing Environment (ISCE) [17] software developed by the Jet Propulsion Laboratory (JPL) and Stanford University to generate interferograms of the three The black rectangle is the InSAR measurement region in this study. The focal mechanism solutions are from the USGS, GCMT and CSI.

InSAR Coseismic Deformation Measurements
The SAR data used in this study included the Single Look Complex (SLC) data in ascending and descending orbits from the Sentinel-1A/B and RADARSAT-2 satellite platforms (Table 1), operating at the C band with a radar wavelength of about 5.6 cm. The coseismic interferometry pairs of the Sentinel-1 ascending, descending, and RADARSAT-2 ascending orbits are referred to the S1A, S1D, and R2A in this study, respectively. A and D represent the ascending orbit and descending orbit, respectively. The superscript A and B represent the Sentinel-1A and Sentinel-1B satellites, respectively. ∆T is the time baseline with a unit of day. B ⊥ is the vertical baseline with a unit in meters. ∆H 2π is the height of ambiguity with a unit in meters. Abbr. is the abbreviation of corresponding interferogram pairs.
Remote Sens. 2018, 10, 1400 4 of 19 We used the InSAR Scientific Computing Environment (ISCE) [17] software developed by the Jet Propulsion Laboratory (JPL) and Stanford University to generate interferograms of the three image pairs ( Table 1). The precise orbit data for S1A and S1D were used to reduce potential orbit errors. A Shuttle Radar Topography Mission (SRTM) digital elevation model with a spatial resolution of about 30 m [18] was used to remove the topographic contribution of the interferograms. The interferograms were then unwrapped with the SNAPHU method [19]. For the Sentinel-1A/1B data, we mosaicked the multiple bursts in the same swath and the three swaths in a track. A strict, high-precision data alignment is required for the Sentinel-1 SAR data. In addition, a spectral method has been proposed to remove the phase jumps due to steep frequency variations along the azimuth [20]. The ISCE can process the Sentinel-1 interferometry efficiently and allows for the extraction of the Region of Interest (ROI) from the entire mosaic image, which can greatly reduce the computation time.
The interferogram pairs have short vertical baselines smaller than 50 m (Table 1), which are beneficial for decreasing the contribution of potential topographic phase errors and maintaining high coherence. The heights of ambiguity of the S1A, S1D, and R2A interferogram pairs for a 1-m error in the terrain can introduce a~6.5 × 10 −3 cm,~2.2 × 10 −3 cm, and~1.1 × 10 −2 cm deformation deviation, respectively. The vertical accuracy of SRTM 1 is generally better than~10 m [21]; therefore, the influence of topography errors is very small and was ignored in this study. In addition, timely satellite observations after the mainshock are beneficial for excluding postseismic deformation contamination. The ionospheric errors were not negligible for the InSAR measurement with a long wavelength of the L or P band; however, they were greatly weakened for the C band with a relatively shorter wavelength [22]. Significant topography correlated atmospheric screens (APS) could be observed in all three interferograms. To reduce the effects of APS, independent APS simulations with a GACOS online system [23,24] were used to correct the interferograms for APS. We found that the tropospheric disturbance of the S1D interferograms was effectively reduced, and the variance of the far field noise in the interferogram decreased from 52 mm 2 to 41 mm 2 .
The coseismic deformation fields of the InSAR measurements are shown in Figure 2. The coseismic deformation fields can be divided into east and west zones by the nearly NW-trending boundary. The west zone has a larger deformation magnitude and coverage than the east zone. The western deformation fields of S1A and R2A with ascending orbits demonstrated subsidence in the line-of-sight (LOS) direction, opposite to those from the S1D with descending orbit, which indicated that the predominant deformation was a horizontal component as the major difference between the ascending and descending InSAR measurements was due to the contributions of horizontal motions. Therefore, the InSAR results in multiple viewing geometries, particularly in the descending and ascending directions, can better serve the determination of the focal mechanism of the Jiuzhaigou earthquake. image pairs (Table 1). The precise orbit data for S1A and S1D were used to reduce potential orbit errors. A Shuttle Radar Topography Mission (SRTM) digital elevation model with a spatial resolution of about 30 m [18] was used to remove the topographic contribution of the interferograms. The interferograms were then unwrapped with the SNAPHU method [19]. For the Sentinel-1A/1B data, we mosaicked the multiple bursts in the same swath and the three swaths in a track. A strict, high-precision data alignment is required for the Sentinel-1 SAR data. In addition, a spectral method has been proposed to remove the phase jumps due to steep frequency variations along the azimuth [20]. The ISCE can process the Sentinel-1 interferometry efficiently and allows for the extraction of the Region of Interest (ROI) from the entire mosaic image, which can greatly reduce the computation time.
The interferogram pairs have short vertical baselines smaller than 50 m (Table 1), which are beneficial for decreasing the contribution of potential topographic phase errors and maintaining high coherence. The heights of ambiguity of the S1A, S1D, and R2A interferogram pairs for a 1-m error in the terrain can introduce a ~6.5 × 10 −3 cm, ~2.2 × 10 −3 cm, and ~1.1 × 10 −2 cm deformation deviation, respectively. The vertical accuracy of SRTM 1 is generally better than ~10 m [21]; therefore, the influence of topography errors is very small and was ignored in this study. In addition, timely satellite observations after the mainshock are beneficial for excluding postseismic deformation contamination. The ionospheric errors were not negligible for the InSAR measurement with a long wavelength of the L or P band; however, they were greatly weakened for the C band with a relatively shorter wavelength [22]. Significant topography correlated atmospheric screens (APS) could be observed in all three interferograms. To reduce the effects of APS, independent APS simulations with a GACOS online system [23,24] were used to correct the interferograms for APS. We found that the tropospheric disturbance of the S1D interferograms was effectively reduced, and the variance of the far field noise in the interferogram decreased from 52 mm 2 to 41 mm 2 .
The coseismic deformation fields of the InSAR measurements are shown in Figure 2. The coseismic deformation fields can be divided into east and west zones by the nearly NW-trending boundary. The west zone has a larger deformation magnitude and coverage than the east zone. The western deformation fields of S1A and R2A with ascending orbits demonstrated subsidence in the line-of-sight (LOS) direction, opposite to those from the S1D with descending orbit, which indicated that the predominant deformation was a horizontal component as the major difference between the ascending and descending InSAR measurements was due to the contributions of horizontal motions. Therefore, the InSAR results in multiple viewing geometries, particularly in the descending and ascending directions, can better serve the determination of the focal mechanism of the Jiuzhaigou earthquake. Figure 2. The InSAR coseismic deformation fields. (a-c) Corresponding coseismic deformation fields of the interferogram pairs of S1A, S1D and R2A, respectively. Figure 2. The InSAR coseismic deformation fields. (a-c) Corresponding coseismic deformation fields of the interferogram pairs of S1A, S1D and R2A, respectively.

Geodetic Modeling
First, we estimated the fault geometry parameters from down-sampling the InSAR data. Then, we jointly inverted the InSAR and GPS data to obtain the fault slip distribution as well as its uncertainties.

Data Down-Sampling
The conventional quadtree down-sampling method may produce high densities of data far from the deformation center in the presence of non-deformation noise, decorrelation, or phase-unwrapping errors, which may possibly affect the solution accuracy of fault geometric parameters and slip distributions. Here, we used the Resolution Based method (R_Based) [25] to down-sample the coseismic deformation fields. The R_Based method considers the resolution between the deformation data and fault slip model. It uses a priori fault plane source to determine the data resolution matrix and down-sample the data with spatial variables in which it becomes possible to down-sample the data densely in the near-field of the seismogenic fault and sparsely in the far-field (Figure 3), which can effectively decrease the far-field noise influence in the following inversion. Through R_Based resampling, we obtained 1953, 1581, and 1517 deformation points from the interferograms of the S1A, S1D, and R2A, respectively, with variable incidence angles and azimuth angles to the coordinates of InSAR data. First, we estimated the fault geometry parameters from down-sampling the InSAR data. Then, we jointly inverted the InSAR and GPS data to obtain the fault slip distribution as well as its uncertainties.

Data Down-Sampling
The conventional quadtree down-sampling method may produce high densities of data far from the deformation center in the presence of non-deformation noise, decorrelation, or phase-unwrapping errors, which may possibly affect the solution accuracy of fault geometric parameters and slip distributions. Here, we used the Resolution Based method (R_Based) [25] to down-sample the coseismic deformation fields. The R_Based method considers the resolution between the deformation data and fault slip model. It uses a priori fault plane source to determine the data resolution matrix and down-sample the data with spatial variables in which it becomes possible to down-sample the data densely in the near-field of the seismogenic fault and sparsely in the far-field (Figure 3), which can effectively decrease the far-field noise influence in the following inversion. Through R_Based resampling, we obtained 1953, 1581, and 1517 deformation points from the interferograms of the S1A, S1D, and R2A, respectively, with variable incidence angles and azimuth angles to the coordinates of InSAR data.

Nonlinear Inversion for Fault Geometry
According to the focal mechanism of the mainshock and aftershocks of the Jiuzhaigou earthquake [5], the seismic sequence was dominantly controlled by a vertical left-lateral strike-slip fault with a SSE strike direction. The long axis orientation of isoseismic, spatially linear distributions of the aftershocks, and earthquake-triggered bedrock collapses or landslides in the epicenter all indicate that the seismogenic faulting was oriented NNW [4]. The focal mechanism solutions of the GCMT, USGS, and CSI also showed that the SSE-oriented rupture controlled this event. The long axis azimuth of the InSAR deformation fringes coincided with the SSE-oriented rupture. Therefore, we assumed that the seismogenic fault was an NNW-SSE striking and left-lateral strike-slip fault, which has also been adopted by most of the preceding researches [6,[8][9][10]12].
To build the fault model more reasonably, we adopted two methods to invert the fault geometry parameters. If the geometry parameters inverted by different methods were similar, we believed that the results were reliable and could be cross-validated. However, if the opposite was true, the results were unreliable.

Nonlinear Inversion for Fault Geometry
According to the focal mechanism of the mainshock and aftershocks of the Jiuzhaigou earthquake [5], the seismic sequence was dominantly controlled by a vertical left-lateral strike-slip fault with a SSE strike direction. The long axis orientation of isoseismic, spatially linear distributions of the aftershocks, and earthquake-triggered bedrock collapses or landslides in the epicenter all indicate that the seismogenic faulting was oriented NNW [4]. The focal mechanism solutions of the GCMT, USGS, and CSI also showed that the SSE-oriented rupture controlled this event. The long axis azimuth of the InSAR deformation fringes coincided with the SSE-oriented rupture. Therefore, we assumed that the seismogenic fault was an NNW-SSE striking and left-lateral strike-slip fault, which has also been adopted by most of the preceding researches [6,[8][9][10]12].
To build the fault model more reasonably, we adopted two methods to invert the fault geometry parameters. If the geometry parameters inverted by different methods were similar, we believed that the results were reliable and could be cross-validated. However, if the opposite was true, the results were unreliable.

MPSO+MC Method
Based on the hypothesis of uniform slip, the MPSO algorithm [13] was used to invert for the geometric parameters of the fault including the fault location, length, width, strike, buried depth, dip angle, slip magnitude, and rake angle. The MPSO nonlinear inversion algorithm has been successfully applied to multiple earthquakes [26,27]. The relative weights of the coseismic deformation fields of S1A, S1D, and R2A were separately determined at 0.17, 0.43, and 0.40 by the far-field noise levels, which were evaluated by the best-fit 1D covariance function [14].
By minimizing the root mean square (RMS) deviation between the measurements and model values, an optimal solution of the geometric parameters of the fault can be obtained. However, the MPSO method alone cannot estimate the inversion errors of the fault parameters. Referring to the MC method [14], we further generated 100 groups of InSAR measurement data with random noise through the full Variance-Covariance Matrix (VCM) in the far field, assuming that the noise was isotropic in the image. With regard to the constraints imposed by the noisy data, we again performed global search 100 times to obtain 100 sets of solutions of fault parameters. Then, the mean and standard deviation of the fault parameters were solved by the Gaussian distribution statistics of the 100 sets of solutions ( Figure 4). The standard deviation was considered as the error of the fault geometry parameter.

MPSO+MC Method
Based on the hypothesis of uniform slip, the MPSO algorithm [13] was used to invert for the geometric parameters of the fault including the fault location, length, width, strike, buried depth, dip angle, slip magnitude, and rake angle. The MPSO nonlinear inversion algorithm has been successfully applied to multiple earthquakes [26,27]. The relative weights of the coseismic deformation fields of S1A, S1D, and R2A were separately determined at 0.17, 0.43, and 0.40 by the far-field noise levels, which were evaluated by the best-fit 1D covariance function [14].
By minimizing the root mean square (RMS) deviation between the measurements and model values, an optimal solution of the geometric parameters of the fault can be obtained. However, the MPSO method alone cannot estimate the inversion errors of the fault parameters. Referring to the MC method [14], we further generated 100 groups of InSAR measurement data with random noise through the full Variance-Covariance Matrix (VCM) in the far field, assuming that the noise was isotropic in the image. With regard to the constraints imposed by the noisy data, we again performed global search 100 times to obtain 100 sets of solutions of fault parameters. Then, the mean and standard deviation of the fault parameters were solved by the Gaussian distribution statistics of the 100 sets of solutions ( Figure 4). The standard deviation was considered as the error of the fault geometry parameter.  The geometrical parameters of the seismogenic fault retrieved by the MPSO+MC algorithm are shown in Table 2, and represented by F1. According to the MPSO+MC inversion, the strike, dip angle and rake angle of the fault were closer to the GCMT and USGS solutions than the CSI solutions. Assuming the shear modulus was µ = 32 GPa, the seismic moment of the Jiuzhaigou earthquake derived by the MPSO+MC algorithm was about (5.54 ± 0.71) × 10 18 Nm, equivalent to a moment magnitude of Mw 6.44 ± 0.03, slightly smaller than those of the focal mechanism solutions.

MCMC Method
One problem of the Akaike's Bayesian Information Criterion (ABIC) inversion algorithm [28] is that the model space integration cannot be correctly resolved under non-negative constraints. Therefore, one full Bayesian inversion method has been proposed, which adopts the more reasonable sampling method of the Markov Chain Monte Carlo [29]. However, if the priori constraints are too strong, the solution of the MCMC method will be not unique, which may cause the multimodal posterior distribution to possibly fall into a local minimum solution. The maximum value of complex functions with multiple peaks can be found by the simulated annealing (SA) algorithm, which can effectively avoid falling into the local minimum solution [15]. Furthermore, the MCMC method improved by the SA algorithm can make the acceptance probability cool and converge faster.
Although the above MCMC method was suitable for the optimization problem with non-negative constraints, we adopted the method to invert the geometric parameters of the fault with the boundary constraints of variables. The Markov chain can be divided into the burn-in and resampling stage. A relatively long Markov chain is needed to ensure the solution convergence. The total length of the Markov chain was composed of 1 × 10 6 seeds, the length of the burning and resampling stage were separately composed of 0.5 × 10 6 seeds, and the resampling interval was set as 100 seeds. The advantage of the MCMC algorithm is that not only can the mean value of geometric parameters be retrieved, but also the error range and marginal probability density function (MPDF) of each parameter can be estimated. After resampling the second stage of the Markov chain, the mean model can be obtained by averaging the sampling space. The results showed that the MPDF of the fault parameters were consistent with the normal Gauss distribution ( Figure 5). The geometrical parameters of the seismogenic fault retrieved by the MPSO+MC algorithm are shown in Table 2, and represented by F1. According to the MPSO+MC inversion, the strike, dip angle and rake angle of the fault were closer to the GCMT and USGS solutions than the CSI solutions. Assuming the shear modulus was μ = 32 GPa, the seismic moment of the Jiuzhaigou earthquake derived by the MPSO+MC algorithm was about (5.54 ± 0.71) × 10 18 Nm, equivalent to a moment magnitude of Mw 6.44 ± 0.03, slightly smaller than those of the focal mechanism solutions.

MCMC Method
One problem of the Akaike's Bayesian Information Criterion (ABIC) inversion algorithm [28] is that the model space integration cannot be correctly resolved under non-negative constraints. Therefore, one full Bayesian inversion method has been proposed, which adopts the more reasonable sampling method of the Markov Chain Monte Carlo [29]. However, if the priori constraints are too strong, the solution of the MCMC method will be not unique, which may cause the multimodal posterior distribution to possibly fall into a local minimum solution. The maximum value of complex functions with multiple peaks can be found by the simulated annealing (SA) algorithm, which can effectively avoid falling into the local minimum solution [15]. Furthermore, the MCMC method improved by the SA algorithm can make the acceptance probability cool and converge faster.
Although the above MCMC method was suitable for the optimization problem with non-negative constraints, we adopted the method to invert the geometric parameters of the fault with the boundary constraints of variables. The Markov chain can be divided into the burn-in and resampling stage. A relatively long Markov chain is needed to ensure the solution convergence. The total length of the Markov chain was composed of 1 × 10 6 seeds, the length of the burning and resampling stage were separately composed of 0.5 × 10 6 seeds, and the resampling interval was set as 100 seeds. The advantage of the MCMC algorithm is that not only can the mean value of geometric parameters be retrieved, but also the error range and marginal probability density function (MPDF) of each parameter can be estimated. After resampling the second stage of the Markov chain, the mean model can be obtained by averaging the sampling space. The results showed that the MPDF of the fault parameters were consistent with the normal Gauss distribution ( Figure 5). The geometrical parameters of the seismogenic fault retrieved by the MCMC algorithm are shown in Table 2, and represented by F2. Taking the shear modulus as µ = 32 GPa, the seismic moment of the Jiuzhaigou earthquake was about (5.43 ± 0.10) × 10 18 Nm, equivalent to the moment magnitude Mw 6.43 ± 0.01, which was very close to the MPSO+MC algorithm.

Inversion for Fault Slip Distribution
Except for about a 3 • deviation in the dip angle of the F1 and F2 models, the other parameters were highly consistent with each other. Due to the high similarity of the fault geometry parameters determined by the different nonlinear inversion methods, we believed that the fault geometry models were reliable. As there was only a slight difference between the F1 and F2 models, the F1 model was considered for inversion for the fault slip distribution. Besides the InSAR data, we also adopted the coseismic displacements of 10 GPS stations [8] to constrain the slip distribution.
The fault model of F1 extended along the strike and dip directions to a length of 40 km and 22 km, respectively. In order to discretize the fault plane into a reasonable patch size, we first performed checkerboard tests. Synthetic slip distributions with patch sizes of 2 km, 4 km and 6 km were used as input models to calculate LOS displacements at the multiple view geometries of the InSAR and GPS locations, and then the synthetic displacements were subsequently used to invert for the distributed slip basing on the same regularization that was applied in the forward-model. The checkerboard tests showed that the slips down to a depth of 18 km were well recovered, with a patch size of 6 km ( Figure S1). According to previous studies [30,31], a uniform fault patch size may be experimentally about 20% of the recoverable spatial resolution of the checkerboard at the depth, above which most of slip occurs. Therefore, the extended fault plane was divided into 880 patches, each with a size of 1 km × 1 km in this study.
Based on the result from the field investigation [4], we assumed there was no sliding at the top boundary of the fault plane. The Crust1.0 model [32] and EDGRN software, a layered half-space model [33] were employed to compute Green's functions. Considering that the accuracy of GPS is better than that of InSAR, we set its initial relative weight to 1 and set the relative weight ratio of InSAR to (0, 1) with the same value used in the nonlinear inversion of fault geometry. Then, by using the mean RMS residuals ratio of GPS and the InSAR measurement in the initial inversion of the slip distribution, an optimal factor of 10:1 was chosen as the relative weight ratio of the GPS and InSAR data after multiple calculations.
The reasonable smoothing factor was set to 0.35, which was determined by the trade-off curve of the relative fitting residuals of the Misfit and Roughness (Figure 6a). In order to further optimize the dip angle, the slip distributions with a variable dip angle from 70 • to 80 • with an interval of 0.5 • were inverted, separately; and the optimized dip angle was determined to be 77 • by the trade-off curve between the Misfit and dip angle (Figure 6b). The geometrical parameters of the seismogenic fault retrieved by the MCMC algorithm are shown in Table 2, and represented by F2. Taking the shear modulus as μ = 32 GPa, the seismic moment of the Jiuzhaigou earthquake was about (5.43 ± 0.10) × 10 18 Nm, equivalent to the moment magnitude Mw 6.43±0.01, which was very close to the MPSO+MC algorithm.

Inversion for Fault Slip Distribution
Except for about a 3° deviation in the dip angle of the F1 and F2 models, the other parameters were highly consistent with each other. Due to the high similarity of the fault geometry parameters determined by the different nonlinear inversion methods, we believed that the fault geometry models were reliable. As there was only a slight difference between the F1 and F2 models, the F1 model was considered for inversion for the fault slip distribution. Besides the InSAR data, we also adopted the coseismic displacements of 10 GPS stations [8] to constrain the slip distribution.
The fault model of F1 extended along the strike and dip directions to a length of 40 km and 22 km, respectively. In order to discretize the fault plane into a reasonable patch size, we first performed checkerboard tests. Synthetic slip distributions with patch sizes of 2 km, 4 km and 6 km were used as input models to calculate LOS displacements at the multiple view geometries of the InSAR and GPS locations, and then the synthetic displacements were subsequently used to invert for the distributed slip basing on the same regularization that was applied in the forward-model. The checkerboard tests showed that the slips down to a depth of 18 km were well recovered, with a patch size of 6 km ( Figure S1). According to previous studies [30,31], a uniform fault patch size may be experimentally about 20% of the recoverable spatial resolution of the checkerboard at the depth, above which most of slip occurs. Therefore, the extended fault plane was divided into 880 patches, each with a size of 1 km × 1 km in this study.
Based on the result from the field investigation [4], we assumed there was no sliding at the top boundary of the fault plane. The Crust1.0 model [32] and EDGRN software, a layered half-space model [33] were employed to compute Green's functions. Considering that the accuracy of GPS is better than that of InSAR, we set its initial relative weight to 1 and set the relative weight ratio of InSAR to (0, 1) with the same value used in the nonlinear inversion of fault geometry. Then, by using the mean RMS residuals ratio of GPS and the InSAR measurement in the initial inversion of the slip distribution, an optimal factor of 10:1 was chosen as the relative weight ratio of the GPS and InSAR data after multiple calculations.
The reasonable smoothing factor was set to 0.35, which was determined by the trade-off curve of the relative fitting residuals of the Misfit and Roughness (Figure 6a). In order to further optimize the dip angle, the slip distributions with a variable dip angle from 70° to 80° with an interval of 0.5° were inverted, separately; and the optimized dip angle was determined to be 77° by the trade-off curve between the Misfit and dip angle (Figure 6b). The uncertainty of the slip distributions could not be directly obtained by the single SDM inversion. Therefore, we integrated the MC method with the SDM method. First, the SDM inversion was performed 100 times, each time with the input data generated by observational LOS displacement plus random noise using the VCM method in the far field. Then the 100 sets of slip distribution solutions were analyzed statistically to calculate the average and standard deviation of the slip distributions ( Figure 7 and Table 2). The uncertainty of the slip distributions could not be directly obtained by the single SDM inversion. Therefore, we integrated the MC method with the SDM method. First, the SDM inversion was performed 100 times, each time with the input data generated by observational LOS displacement plus random noise using the VCM method in the far field. Then the 100 sets of slip distribution solutions were analyzed statistically to calculate the average and standard deviation of the slip distributions ( Figure 7 and Table 2). Slip distribution inverted based on the combination of the Sentinel-1A/1B and Radarsat-2 data, which the slip distribution in the lower left is more reasonable than (a); (c) Slip distribution inverted based on the combination of the Sentinel-1A/1B, Radarsat-2, and GPS data, which was the optimum solution with the S1 model in Table 2; (d) Uncertainties with the slip distribution of the S1 model. The maximum uncertainty of the fault sliding was about 3.2 cm. The maximal fault sliding was about 1.06 m. The red pentagram is the mainshock by Fang et al. [34].
Using the distributing slip model of S1 to forward calculate the deformation fields and to compare with the InSAR and GPS data, the residuals obtained are shown in Figure 8. The mean residuals of the retrieved S1A, S1D, and R2A interferogram pairs were about 1.9 cm, 1.4 cm, and 0.8 cm, respectively. The mean residual of the GPS horizonal deformation was about 2.0 mm. Overall, the residuals were not systematic spatially, indicating that the fault slip distributions model was a good match for the geodetic measurement data. (b) Slip distribution inverted based on the combination of the Sentinel-1A/1B and Radarsat-2 data, which the slip distribution in the lower left is more reasonable than (a); (c) Slip distribution inverted based on the combination of the Sentinel-1A/1B, Radarsat-2, and GPS data, which was the optimum solution with the S1 model in Table 2; (d) Uncertainties with the slip distribution of the S1 model. The maximum uncertainty of the fault sliding was about 3.2 cm. The maximal fault sliding was about 1.06 m. The red pentagram is the mainshock by Fang et al. [34].
Using the distributing slip model of S1 to forward calculate the deformation fields and to compare with the InSAR and GPS data, the residuals obtained are shown in Figure 8. The mean residuals of the retrieved S1A, S1D, and R2A interferogram pairs were about 1.9 cm, 1.4 cm, and 0.8 cm, respectively. The mean residual of the GPS horizonal deformation was about 2.0 mm. Overall, the residuals were not systematic spatially, indicating that the fault slip distributions model was a good match for the geodetic measurement data. The InSAR residuals of S1A, S1D, and R2A, respectively; (d) Red and blue arrows are the observation and simulation deformation of GPS, separately. The red solid line is the fault model of S1 in the surface. The maximal residuals were smaller than 0.1 m and the mean residuals were separately 1.9 cm, 1.4 cm, and 0.8 cm for the S1A, S1D, and R2A. The simulation deformations of GPS were generally consistent with the observations with a mean residual of about 2.0 mm.

Characters of the Slip Distribution Model
The fault slip was mainly concentrated in a rectangular area 30 km long along the fault strike direction and 16 km wide along the dip direction, with a major slip at a depth of 1-15 km within the recoverable depth range of the checkerboard test, and one slip center at a depth of about 5-9 km (Figure 7c). The slip within the depth of 0-1 km was smaller than 0.25 m. The InSAR residuals of S1A, S1D, and R2A, respectively; (d) Red and blue arrows are the observation and simulation deformation of GPS, separately. The red solid line is the fault model of S1 in the surface. The maximal residuals were smaller than 0.1 m and the mean residuals were separately 1.9 cm, 1.4 cm, and 0.8 cm for the S1A, S1D, and R2A. The simulation deformations of GPS were generally consistent with the observations with a mean residual of about 2.0 mm.

Characters of the Slip Distribution Model
The fault slip was mainly concentrated in a rectangular area 30 km long along the fault strike direction and 16 km wide along the dip direction, with a major slip at a depth of 1-15 km within the recoverable depth range of the checkerboard test, and one slip center at a depth of about 5-9 km (Figure 7c). The slip within the depth of 0-1 km was smaller than 0.25 m. The F1 and F2 are separately the uniform slip fault model retrieved by the MPSO+MC and the MCMC method; the S1 are the distributed slip models derived by the SDM+MC method based on the extension of F1 with an optimize dip angle. Models from References [6,7] are with the inversion of the InSAR data; the model of Reference [5] is with the inversion of seismic waves; the model in Reference [8] is with the inversion of the GPS data; the model from Reference [12] is derived from the joint inversion of the InSAR and seismic wave data; the models in References [9,10] are with the joint inversion of the InSAR and GPS data. The superscript C represents the center point of the fault plane, the superscript A represents the average value, the superscript M represents the maximum, and the superscript MS represents the maximal slip.
The maximum uncertainty of the slip model was about 3.2 cm (Figure 7d), equivalent to about 3.0% of its maximum slip (about 1.06 m). The average slip estimated by the MC method was about 0.28 ± 0.01 m, with an average rake angle about −7.86 ± 0.31 • ; the maximum slip was about 1.06 ± 0.02 m, located at a depth of 6.84 ± 0.33 km, with a rake angle about −10.37 ± 0.34 • . When the shear modulus was µ = 32 GPa, the seismic moment was about (7.85 ± 0.13) × 10 18 Nm, which was equivalent to a moment magnitude of about Mw 6.54 ± 0.01. The slip models showed that the Jiuzhaigou earthquake was a typical sinistral strike-slip rupture event.

Relationship between Fault Slip and Aftershocks
According to the precise location of aftershocks [34] during the month after the mainshock, the aftershocks were distributed over a zone of~42 km with NWW-SEE trending about 150 • , predominantly at a depth of 4-20 km. In addition, the dip angle of the seismogenic fault was steep with an average value of about 84 • . Finally, the average positioning error of the aftershocks in the three directions of E-W, N-S, and U-D was 0.16, 0.15, and 0.18 km, respectively.
To reveal the relationship between the slip distribution and the aftershocks, we projected the positions of the aftershocks to the parallel and vertical direction of the fault trending, respectively ( Figure 9). We found that there was a high consistency between the spatial distribution of the aftershocks and fault slip in the section parallel to the fault trending from NW to SE (Figure 9a). However, it is worth noting that the aftershocks were relatively scarce in the sliding peak center within the depth of 5-10 km, which may have been due to a full rupture of the mainshock in this region.
In addition, the fault geometric model essentially passed through the concentration zone of aftershocks in the section perpendicular to the fault trending from southwest (SW) to northeast (NE) (Figure 9b). There is also a need to note that the optimal dip angle of 77 • was slightly smaller than that (~84 • ) directly suggested by the aftershocks [34]. From the multiple view geometries of the InSAR deformation fields (Figure 2), the deformation range and amplitude west of the fault were obviously larger than that in the east. In the case of a purely sinistral strike-slip rupture with vertical dip angle, the elastic deformation on both sides of the fault should be completely symmetrical. Apparently, the deformation of the Jiuzhaigou earthquake with a typical sinistral strike-slip rupture did not show symmetrical characteristics. Therefore, we suggest that the 77 • dip angle was possibly more reasonable than the nearly vertical dip angle, considering the coincidence with the geodetic observation. The maximum uncertainty of the slip model was about 3.2 cm (Figure 7d), equivalent to about 3.0% of its maximum slip (about 1.06 m). The average slip estimated by the MC method was about 0.28 ± 0.01 m, with an average rake angle about −7.86 ± 0.31°; the maximum slip was about 1.06 ± 0.02 m, located at a depth of 6.84 ± 0.33 km, with a rake angle about −10.37 ± 0.34°. When the shear modulus was μ = 32 GPa, the seismic moment was about (7.85 ± 0.13) × 10 18 Nm, which was equivalent to a moment magnitude of about Mw 6.54 ± 0.01. The slip models showed that the Jiuzhaigou earthquake was a typical sinistral strike-slip rupture event.

Relationship between Fault Slip and Aftershocks
According to the precise location of aftershocks [34] during the month after the mainshock, the aftershocks were distributed over a zone of ~42 km with NWW-SEE trending about 150°, predominantly at a depth of 4-20 km. In addition, the dip angle of the seismogenic fault was steep with an average value of about 84°. Finally, the average positioning error of the aftershocks in the three directions of E-W, N-S, and U-D was 0.16, 0.15, and 0.18 km, respectively.
To reveal the relationship between the slip distribution and the aftershocks, we projected the positions of the aftershocks to the parallel and vertical direction of the fault trending, respectively ( Figure 9). We found that there was a high consistency between the spatial distribution of the aftershocks and fault slip in the section parallel to the fault trending from NW to SE (Figure 9a). However, it is worth noting that the aftershocks were relatively scarce in the sliding peak center within the depth of 5-10 km, which may have been due to a full rupture of the mainshock in this region.
In addition, the fault geometric model essentially passed through the concentration zone of aftershocks in the section perpendicular to the fault trending from southwest (SW) to northeast (NE) (Figure 9b). There is also a need to note that the optimal dip angle of 77° was slightly smaller than that (~84°) directly suggested by the aftershocks [34]. From the multiple view geometries of the InSAR deformation fields (Figure 2), the deformation range and amplitude west of the fault were obviously larger than that in the east. In the case of a purely sinistral strike-slip rupture with vertical dip angle, the elastic deformation on both sides of the fault should be completely symmetrical. Apparently, the deformation of the Jiuzhaigou earthquake with a typical sinistral strike-slip rupture did not show symmetrical characteristics. Therefore, we suggest that the 77° dip angle was possibly more reasonable than the nearly vertical dip angle, considering the coincidence with the geodetic observation.  [34]; the slip distribution model is S1; the red line is the fault location in the section. The fault geometry and slip distribution had a high consistency with the aftershocks.  [34]; the slip distribution model is S1; the red line is the fault location in the section. The fault geometry and slip distribution had a high consistency with the aftershocks.

Coulomb Stress Change by the Mainshock
According to the seismic stress trigger theory, the fault rupture of an earthquake will always lead to stress disturbance at the epicenter and its surrounding zones. According to previous research [35][36][37], the stress change is perhaps highly correlated with the distribution of immediate aftershocks and has an effect on the seismic potentiality of the surrounding faults. According to the Coulomb failure criterion [38], the Coulomb stress is written as: where σ f is the Coulomb stress change on a specific receiving fault; τ β (positive in the direction of fault slip) and σ β (positive in compression) are the shear stress and normal stress change on the receiving fault, respectively; µ is the effective friction coefficient with a range 0-1. Using the fault slip distribution derived from the InSAR and GPS data, we calculated the coseismic Coulomb stress changes using the software package PSCMP/PSGRN [39]. The coseismic Coulomb stress changes with different depths of 0 to 20 km with an interval of 5 km were calculated using the optimal rupture plane of F1 with empirical friction of 0.4 for the receiving fault to reveal the triggering relationship between the mainshock and aftershocks ( Figure 10).  [34]. The far-field aftershocks denoted by rose-red solid circles are from CSI. The magnitude of near-field aftershocks is referred to in Figure 9b. The gray lines are faults referred to in Figure 1. The depth of the aftershocks is <5 km with a, ≥5 km and <10 km with b, ≥10 km and <15 km with c, ≥15 km and <20 km with d, and ≥20 km with e, respectively.  [34]. The far-field aftershocks denoted by rose-red solid circles are from CSI. The magnitude of near-field aftershocks is referred to in Figure 9b. The gray lines are faults referred to in Figure 1. The depth of the aftershocks is <5 km with a, ≥5 km and <10 km with b, ≥10 km and <15 km with c, ≥15 km and <20 km with d, and ≥20 km with e, respectively.
The dislocation model is based on linear elasticity where only infinitesimal strain is assumed. In the vicinity of the rupturing surface, especially on their edges, an extremely large stress concentration occurs, and linear elasticity does not work. Due to the limitation of the dislocation model, the Coulomb stress changes in the near-field of the fault are unreliable; therefore, we do not discuss the spatial correlation between the near-field aftershocks and Coulomb stress changes.
The twenty-three aftershocks (M > 2) in the far-field of the fault within one year after the mainshock from 8 August 2017 to 7 July 2018 were compared with the Coulomb stress changes ( Figure 10). We found that only five shocks occurred in the stress shadow area (Figure 10a-d), about 22% of the total number shocks in the far-field. Additionally, about 78% of the shocks happened in the positive stress change area, which were obviously consistent with the stress trigger theory. A notable phenomenon was that a total of eight shocks happened in the vicinity of the NHY fault ( Figure 10).
To evaluate the effect of effective friction on the Coulomb stress changes, we also calculated the Coulomb stress changes with the frictions of 0.1 and 0.7 (Figures S2 and S3). The comparison showed that the different friction had a few influences on the Coulomb stress changes. However, the spatial correlation between the distribution of the off-fault aftershocks and Coulomb stress changes was not altered by the apparently different friction. Most of the aftershocks in the far-field happened in the positive stress change area with the different frictions.

Stress Disturbance at the Surrounding Faults
The Coulomb stress disturbances of the surrounding faults have been calculated in previous studies [7,40,41]. As our slip distribution model was somewhat different from previous models (Table 2), therefore, the stress disturbance to the surrounding faults due to the coseismic rupture needed to be reanalyzed for the receiving faults (Table S1). We not only considered the receiving faults [42], but also the Tazang fault and North Huya faults near the mainshock. The Tazang fault is considered as a southeastern extension of the East Kunlun fault. The North Huya fault is a hidden northern segment of the Huya fault. The strike direction, dip angle, and rake angle of the North Huya fault were assumed to be those of the mainshock fault model in the Section 3.2.
The stress distributions were calculated by the software PSCMP/PSGRN [39] with a stratified crustal model of Crust1.0. The stress distributions at the depth of 5 km showed that the remarkable stress increase (>0.1 MPa) at the Tazang fault was primarily located in the intersection with the Minjiang fault (Figure 11a). The stress increase of the Minjiang fault was localized west of the epicenter. The overall stress of the Huya fault was increased, especially in the northern hidden section, with a stress increase larger than 0.01 MPa. The stress disturbances to other faults were within ±0.01 MPa, approximate to the previous results [7,40,41]. The uncertainties of the stress distribution were also calculated based on the errors of the slip distribution model. Except for two small fragments in the middle-segment of the Tazang fault and the northern-segment of the north Huya fault (Figure 11b), most of the uncertainties were smaller than 0.003 MPa, which were within a reasonable range. section, with a stress increase larger than 0.01 MPa. The stress disturbances to other faults were within ±0.01 MPa, approximate to the previous results [7,40,41]. The uncertainties of the stress distribution were also calculated based on the errors of the slip distribution model. Except for two small fragments in the middle-segment of the Tazang fault and the northern-segment of the north Huya fault (Figure 11b), most of the uncertainties were smaller than 0.003 MPa, which were within a reasonable range.  [34]. The fault names are referred to in Figure 1 and the fault parameters are referred to in Table S1.

Seismic Potentiality of the Near-Field Fault
The static Coulomb stress was significantly enhanced (>0.01 MPa) in the northwestern Tazang fault that intersects with the Minjiang fault, with a static Coulomb stress increase at a segment more than 0.03 MPa, three times larger than the threshold to trigger future earthquakes. The static Coulomb stress also significantly increased in partial sections of the Minjiang fault, west of the epicenter with a maximum value of more than 0.03 MPa. However, the static Coulomb stress reduced in other sections of the Minjiang fault ( Figure 11a). Therefore, this suggests that the Jiuzhaigou earthquake may increase the seismic potentiality of the northwestern section of the Tazang fault and partial sections of the Minjiang fault.
Some researchers have claimed that if the collective effect of the 2017 Jiuzhaigou Mw 6.5 earthquake, 1973 Huanglong MW 6.5 earthquake, and the 1976 Songpan earthquake clusters break the Huya fault entirely from south to north, the seismic potentiality of the Huya fault may possibly be reduced [6,9]. However, other studies have found that the static Coulomb stress in the northern part of the Huya fault increased [40,43], and the latest activity of the Huya fault could possibly cause further earthquakes with a relatively small magnitude [11]. Our result confirmed a significant increase in the static Coulomb stress in the North Huya fault due to the Jiuzhaigou earthquake (Figure 11a).  [34]. The fault names are referred to in Figure 1 and the fault parameters are referred to in Table S1.

Seismic Potentiality of the Near-Field Fault
The static Coulomb stress was significantly enhanced (>0.01 MPa) in the northwestern Tazang fault that intersects with the Minjiang fault, with a static Coulomb stress increase at a segment more than 0.03 MPa, three times larger than the threshold to trigger future earthquakes. The static Coulomb stress also significantly increased in partial sections of the Minjiang fault, west of the epicenter with a maximum value of more than 0.03 MPa. However, the static Coulomb stress reduced in other sections of the Minjiang fault ( Figure 11a). Therefore, this suggests that the Jiuzhaigou earthquake may increase the seismic potentiality of the northwestern section of the Tazang fault and partial sections of the Minjiang fault.
Some researchers have claimed that if the collective effect of the 2017 Jiuzhaigou Mw 6.5 earthquake, 1973 Huanglong M W 6.5 earthquake, and the 1976 Songpan earthquake clusters break the Huya fault entirely from south to north, the seismic potentiality of the Huya fault may possibly be reduced [6,9]. However, other studies have found that the static Coulomb stress in the northern part of the Huya fault increased [40,43], and the latest activity of the Huya fault could possibly cause further earthquakes with a relatively small magnitude [11]. Our result confirmed a significant increase in the static Coulomb stress in the North Huya fault due to the Jiuzhaigou earthquake (Figure 11a).

Comparison with Previous Slip Models
We collected the partly existing models of Ji et al. [7], Zhao et al. [9], Nie et al. [10], and Zhang et al. [12] and redrew them in the same color scales ( Figure S4). A visual comparison of the models showed that there were some differences among them.
Our slip model showed that there was no significant slip in the lower left and right corners of the fault, which was different from previous models [7,9]. The maximum slip of our model was smaller than the model of Zhao et al. [9] and larger than the other models [7,10,12]. Unlike the other models, the seismic moment of our model was slightly larger than the focal mechanism solutions of GCMT and USGS (Table 2). Considering the role of aftershocks during the observation of the InSAR data, it would be reasonable to suggest that the seismic moments were slightly larger than the focal mechanism solution.
With the sole constraint of the geodetic data, the fault slip distribution models still differed from one another (Figure S4a-d). The discrepancies are possibly due to the fault geometric models. Our fault geometric parameters were derived from multiple view geometries of the InSAR data with nonlinear inversion methods, which were different from the previous studies [7,9,10,12]. We found that the fault geometry models retrieved by the nonlinear MPSO+MC and MCMC methods were highly consistent and the standard deviations of the fault geometry parameters were controlled within small values.
Overall, the fault slip distribution derived from the geodetic deformation measurements were principally located at the depth of 1-15 km ( Figure S4a-d); whereas the addition of seismic waveform data caused dominant sliding to be located at a deeper depth ( Figure S4e). Although there were some differences among the models, the dominant fault rake angles were all between -20 • and 20 • , indicating that the mainshock was attributable to a left-lateral strike-slip rupture.

Discrepancies of Coulomb Stress Changes of Mainshock Rupture
The coseismic Coulomb stress changes at the depth of 15 km, from the multiple coseismic slip distribution models, were calculated with the same crustal model, and receiving fault of the mainshock rupture ( Figure S5) to evaluate the influence of the coseismic slip distribution models to the static Coulomb stress changes. We found some discrepancies in the coseismic Coulomb stress changes.
Coulomb stress distributions are obviously different from each other in the near-field of the fault ( Figure S5). If we hypothesize that the nonlinear elastic actions in the vicinity of fault are similar, the discrepancies in the five dislocation models may dominantly influence the Coulomb stress distributions in the near-field of the fault. There were also some diversities in the far-field of the fault; however, these diversities will not influence the spatial correlation comparison between the far-field aftershock and Coulomb stress change ( Figure S5).
It is noteworthy that if we ignored the nonlinear elastic action in the vicinity of the fault, the discrepancy of the existing coseismic slip distribution models may bias the spatial relationship analysis of the near-field aftershocks with static Coulomb stress changes. For example, the near-field aftershocks at the depth range of 15-20 km occurred mostly in the positive Coulomb stress area identified with the models of Zhao et al. [9], Zhang et al. [12], and this work ( Figure S5a,b,e). However, they happened mostly in the shadow area with the models of Ji et al. [7] and Nie et al. [10] (Figure S5c,d).

Discrepancies of Coulomb Stress Disturbance to Nearby Faults
Using the different slip distribution models ( Figure S4), the stress disturbances to nearby faults at the depth of 5 km were calculated by the software PSCMP/PSGRN with the same crustal model and receiving faults (Table S1). We found that the differences of the coseismic slip distribution models may bias the Coulomb stress disturbances to nearby faults ( Figure S6). The magnitude of positive stress changes at the MJ fault west of the epicenter was larger than 0.03 MP using the slip models of Nie et al. [10] and this work ( Figure S6a,d). Whereas, it was smaller than 0.03 MPa with the slip models of Ji et al. [7], Zhao et al. [9], and Zhang et al. [12] (Figure S6b,c,e). The positive stress perturbations to the Tazang fault east of the epicenter were not significant when calculated with the slip models of Ji et al. [7], Zhao et al. [9], and this study-however, they were significant based on the estimate obtained from the slip model of Nie et al. [10] and Zhang et al. [12].

Conclusions
We inverted for the seismogenic fault geometry models of the Jiuzhaigou Mw 6.5 earthquake from InSAR data by the MPSO+MC and MCMC nonlinear algorithms, and the slip distribution by the SDM+MC method. The fault geometry and slip model were compared with the distribution of the relocated aftershocks. Then, the static Coulomb stress disturbances to the mainshock rupture plane and nearby faults were analyzed. Finally, the discrepancies of several existing slip distribution models and Coulomb stress changes due to the slip models were discussed. The main findings were as follows: (1) A blind faulting with a length of about 23 km, a width of about 11 km, a strike angle of about 154 • , an optimized dip angle of 77 • , was causative of the Jiuzhaigou earthquake. (2) The fault slip was majorly concentrated at the depth of 1-15 km, and only one slip center appeared at the depth of 5-9 km with a maximum slip of~1.06 m. The average rake angle was about 7.84 • , indicating the pure left-lateral strike-slip fracture of the mainshock. (3) The seismic moment derived from the slip distribution was about 7.85 × 10 18 Nm, equivalent to a moment magnitude about Mw 6.54, which was slightly larger than the previous results and the focal mechanism solutions. (4) Most of the off-fault aftershocks with the magnitude >M 2 within one year after the mainshock happened in the stress positive stress change area, which coincided with the stress triggering theory.  (Table S1) at the depth of 5 km; Table S1: Receiving fault parameters for the calculation of stress disturbance.