InSAR Constrained Downdip and Updip Afterslip Following the 2015 Nepal Earthquake: New Insights into Moment Budget of the Main Himalayan Thrust

: We use ALOS-2 and Sentinel-1 data spanning 2015–2020 to obtain the post-seismic deformation of the 2015 Mw 7.8 Nepal earthquake. ALOS-2 observations reveal that the post-seismic deformation was mainly distributed in four areas. A large-scale uplift deformation occurred in the northern subsidence area of the co-seismic deformation ﬁeld, with a maximum uplift of ~80 mm within 4.5 yr after the mainshock. While in the southern coseismic uplift area, the direction of the post-seismic deformation is generally opposite to the co-seismic deformation. Additionally, two notable deformation areas are located in the region around 29 ◦ N, and near the MFT, respectively. Sentinel-1 observations reveal post-seismic uplift deformation on the north side of the co-seismic deformation ﬁeld with an average rate of ~20 mm/yr in line-of-stght. The kinematic afterslip constrained by InSAR data shows that the frictional slip is distributed in both updip and downdip areas. The maximum cumulative afterslip is 0.35 m in downdip areas, and 0.2 m in the updip areas, constrained by the ALOS measurements. The stress-driven afterslip model shows that the afterslip is distributed in the downdip area with a maximum slip of 0.3 m during the ﬁrst year after the earthquake. Within the 4.5 yr after the mainshock, the estimated moment released by afterslip is ~1.5174 × 10 20 Nm,about 21.2% of that released by the main earthquake.


Introduction
On 25 April 2015, a Mw 7.8 earthquake occurred in Nepal, which is located in the middle of the Himalayan orogenic belt. According to the USGS, the epicenter was at 28.2 • N, 84.7 • E, about 82 km away from Kathmandu. On 12 May, 17 days after the Mw 7.8 earthquake, a Mw 7.3 (27.8 • E, 86.1 • N) earthquake occurred about 150 km southeast of the mainshock epicenter. The Indian plate was subducted northwards toward the Eurasian plate at a speed of 40 mm·a −1 , and a 50% convergence rate was absorbed by the Himalayan orogenic belt, causing uplift, thus shortening and producing a lateral extrusion on the Tibetan Plateau [1]. The Himalayas have become the most active and complex orogenic belt in the world, and Nepal is located in the middle of the collision zone between the Indian and Eurasian plates. The southward thrust nappe tectonic system has developed several thrust fault belts from south to north, i.e., the main front faults, the main boundary faults, and the main central thrust faults (Figure 1a). These faults converge on a low-angle main thrust fault or detachment surface (Figure 1b). The numerical experiments that have been performed demonstrate that a weaker, intermediate decollement is needed to form the accretionary and sharp topographic gradients found between the Lesser and High Himalaya, respectively [2]. Two types of co-seismic slip models have been proposed according to previous studies, one is a low-angle fault plane [3,4] the other is a fault plane composed of a flat mid-crustal ramp [5], which induced the permanent uplift on the mid-crustal ramp [6]. The geodetic data shows that the co-seismic rupture depth range of the Nepal earthquake was 10-15 km, and that the shallow part remained locked. A 10-20 km rupture gap between the Mw 7.8 and Mw 7.3 earthquakes could cause a Mw 6.9-7.0 earthquake in the future [7]. Many studies have been carried out on the post-seismic deformation and afterslip that followed the Nepal earthquake, based on GPS and InSAR data, some studies have shown that the afterslip filled the slip-gap between the Mw 7.8 and Mw 7.3 earthquakes [8][9][10]. The afterslip mainly occurred in the downdip region, while in the updip zone afterslip is small and is still in a locked state [4,[8][9][10]. However, whether the afterslip exists in the updip zone and is characteristic of the spatiotemporal distribution remains controversial. Some studies show that the afterslip is mainly distributed in the downdip zone [4,8,10], however, other studies have shown that afterslip exists in the updip area of co-seismic rupture [8,11].
In addition to afterslip, in terms of the viscoelastic relaxation and poroelastic rebound, previous studies have suggested that the post-seismic deformation was mainly dominated by afterslip in the early period (~2 yr), and that viscoelastic relaxation has taken a leading role in later years [9,10,12]. Tian et al. [9] show that the post-seismic deformation was mainly dominated by downdip afterslip in the first 2 years, with viscoelastic relaxation playing a leading role in the following years. Diao et al. [12] refer that the joint contributions of afterslip and viscoelastic relaxation can better explain the observed post-seismic deformation. Hong & Liu jointly inverted the~3 years of InSAR and GPS post-seismic deformation data and suggested that afterslip may last for a longer time [13]. Zhang et al. used the 5 years of available GPS data and found that in the intermediate-field post-seismic deformation was caused by equal parts of afterslip and viscoelastic relaxation [14].
Therefore, the previous studies remain controversial as to whether afterslip exists in the updip zone. Most of the previous studies are based on GPS data from~2 yr after the earthquake. However, the spatial resolution of the GPS data is insufficient, and few GPS stations are distributed in the southern part of the co-seismic deformation zone, where the Sentinel-1 SAR data also cannot obtain good results because of their low coherence. Therefore, the previous studies hardly observed the post-seismic deformation in the south of the co-seismic deformation zone. It is believed that the updip zone is still locked after the Nepal earthquake, and in this study we use the 4.5 yr ALOS-2 data for after the mainshock, with better coherence, to observe the characteristics of the post-seismic deformation in this area. Meanwhile, we processed the 4.5 yr Sentinel-1 data based on the PS-InSAR (Permanent Scatterer InSAR) technique to obtain the average deformation-rate field, and the short-term post-seismic deformation based on the DInSAR (Differential InSAR) and Stacking technologies, to observe the post-seismic deformation and analyze the presence of afterslip in the updip zone.

InSAR Data and Processing Methods
We used the InSAR data for the 4.5 yr after the Mw 7.8 earthquake, that is, the ALOS-2 PALSAR data of the descending orbit T48 from May 2015 to December 2019, and the Sentinel-1 data of the ascending orbit T85 and descending orbit T19 (Figure 1a). Both types of InSAR data cover the co-seismic deformation field well. The ALOS-2 data-processing method we adopted was: (1) Interferometric pair construction and differential interferogram generation. First, all the images were co-registered to a single first image (20170108). Then, interferometric processing was carried out by the GAMMA software [15]. The SRTM-90 (Shuttle Radar Topography Mission) digital elevation model was used to eliminate the topographic phase, and the adaptive filter, and the minimum cost flow algorithm was used to obtain the unwrapped interferogram based on the same reference pixel [15,16]. We retained the pixel points with coherences higher than 0.1 and the interferograms with perpendicular baselines less than 350 m ( Figure 2); (2) Atmospheric phase removal. The influence of the atmospheric delay phase was removed by using the GACOS (Generic Atmospheric Correction Online Service for InSAR) model generated by the fusion of GPS and ECMWF data by Li Zhenhong's team from Newcastle University [17]; (3) Orbital phase removal. We masked the post-seismic deformation area, and used the quadratic function to simulate the orbit phase. For the expression as shown in (1), (x,y) represents the longitude and latitude coordinates of the pixel points, and a, b, c, d, e, f are constants. Figure 3c shows the simulated orbit phase. We subtracted the atmospheric phase and orbital phase to obtain the post-seismic deformation phase ( Figure 3d); (4) Post-seismic time-series de-formation field. After removing the orbit and topographically related atmospheric delay phase, we calculated the time series LOS deformation, and average velocity fields of the Nepal earthquake within 4.5 yr. We used GIAnT [18] software based on the SBAS-InSAR (Small Baseline Subset InSAR) approach, which is mainly based on the utilization of the least squares method for revealing the linear deformation in the reference and secondary images. The basic formula is (2), where Φ ij represents the phase of the two images, δ ϕ n represents the phase increment of the n and n+1 images, G represents Green's function, d represents the data vector, and m represents the model vector.
W orbit = ax 2 + by 2 + cx + dy + exy + f (1)  The method we used to process the Sentinel-1 data was as follows: The Sentinel-1 InSAR data was from May 2015 to December 2019, i.e., the time span was the same as for the ALOS-2 data, the ascending orbit (T85) and descending orbit (T19) observations both covered the co-seismic deformation field well ( Figure 1a). The SLC file, baseline information and geocoded first-image were obtained through the co-registration to a single first image, then we used Stamps V4.1 [19] for further processing. We set the amplitude dispersion to 0.4 and removed the spatially uncorrelated error. The SRTM3 digital elevation model was imported to correct the topographic phase, and the pixel with a DEM error higher than 20 m was eliminated to improve the accuracy of the PS point. The Goldstein adaptive filter was used to improve the signal-to-noise ratio of the data, and the atmospheric model GACOS was used to remove the atmospheric delay phase and obtain the average rate deformation field.
We used the C-band Sentinel data, which is sensitive to minor deformation, to observe the post-seismic deformation in the short-term baseline to ensure the coherence and credibility of the deformation. For the C-band Sentinel data, the phase delay caused by the ionosphere can be ignored [20], only the influence of the phase delay on the troposphere is considered here. The tropospheric phase delay includes vertical stratification and turbulence. The turbulent portion can be removed by interferogram stacking [21,22], while the vertical stratified portion is related to the terrain elevation. We obtained the unwrapped phase by differential interference of five interferograms. At first, we stacked the unwrapped phase of the five interferograms to obtain the average phase rate (Figure 4a), then we selected the topographically related deformation zone from outside of the post-seismic deformation field. Figure S1 shows the fitting curve for the topography and displacement before and after topographic delay correction, with the deformation being linearly related to the topographic height before correction. We subtracted the simulated topographic phase (Figure 4b) from the average phase rate. Finally, we calculated the post-seismic deformation field based on the corrected average phase rate. The post-seismic deformation was mainly distributed in the north of the co-seismic deformation field, as found with the Stacking method ( Figure 4c). The maximum uplift displacement was~0.13 m and the maximum subsidence displacement was~0.04 m in one and a half months. The characteristics of the deformation field are the same as those observed by Sreejith [23].  Figure 5a. There are four main post-seismic deformation areas, as shown by P1, P2, P3 and P4 in Figure 5, the regions of P2 and P3 are located in the range of the co-seismic deformation field. The large-scale east-west uplift zone shown in P3 occurred in the subsidence area of the co-seismic deformation field (north of the co-seismic deformation field), with a maximum uplift of~80mm in 4.5yr (Figure 7a). The small-scale subsidence region of P2 is distributed in the Kathmandu basin with a magnitude of over 100 mm, where there was an uplift area in the co-seismic deformation field. It may be related to the groundwater exploitation in Kathmandu due to its small size and large displacement. The other two significant deformation areas are distributed in the far-field of the co-seismic deformation field, the uplift deformation area P1 is located near the MFT updip from the co-seismic fracture with a maximum deformation of~40 mm in 4.5 yr. The subsidence deformation area P4 is distributed in the north about 200 km away from the co-seismic deformation field, with a maximum deformation of~40 mm in 4.5 yr. Figures 6 and 7 reflect the dynamic evolution process of the post-seismic deformation fields, while the deformation gradient (velocity) is gradually attenuating, which is similar to the post-seismic deformation process of other large earthquakes [24][25][26]. However, significant differences are shown among the four deformation areas regarding the postseismic deformation process, the deformation attenuation amplitudes of the areas P3 and P4 (AA' CC') are much larger than those of areas P1 and P2 (BB' DD'). The deformation rate was quick for the first 2.5 yr after the earthquake and stable in the later period. This may be related to the frictional and rheological properties of the lower crust and upper mantle, and the degree of the disturbance caused by the co-seismic rupture. In order to highlight the deformation characteristics of different time periods after the mainshock, we calculated the average post-seismic deformation rate in two short periods, i.e., the time intervals of 2.5 yr and 1.5-4.5 yr after the Nepal Mw 7.8 earthquake mainshock, and also compared these with the results of the whole observation time span of 4.5 yr (Figure 8). Figure 8 indicates that although the post-seismic deformation patterns in the three periods are generally consistent, there are obvious differences in deformation amplitude. The average rates for the 2.5 yr (Figure 8b) and 4.5 yr (Figure 8a) after the earthquake are basically the same, while the deformation rate of about 1.5-4.5 yr is very small. This suggests that the observed post-seismic deformation mainly occurred in the first~2.5 years, then it decreases and tends to be stable.   In addition, we use GPS post-seismic deformation observation results [27] in December of the first year after the earthquake from 8 stations distributed in the InSAR deformation field for comparison with our results. The locations of GPS stations for comparison are shown in Figure 5a. The eight GPS sites are distributed in four deformation areas, BTNI in area P1, KKN4 and SLBL in area P2, CHLM and JILZ in area P3, and CANA, YANI and GUCU in area P4 (Figure 5a). The ALOS post-seismic displacement sequences in four deformation areas are all well fitted to the GPS data, which used Formula (3) to project the LOS direction with the corresponding incidence (θ inc ) and heading angles (α azi ) of each site of the descending orbit of ALOS data (Figure 9), indicating the reliability of our post-seismic deformation observations from the ALOS-2 data. The trend of post-seismic deformation obeys the logarithmic function, which can be seen in other earthquake events [28][29][30]. The blue circle represents the GPS time series data; the orange circle represents the ALOS-2 data; the error bar indicates the standard deviation of 64 adjacent pixels; the GPS data is projected onto the LOS direction. Figure 10 shows the average deformation rate field of the Sentinel-1 data in 4.5 yr after the mainshock based on the PS-InSAR method. The ascending and descending data shows roughly the same deformation characteristics similar to the ALOS-2 post-seismic deformation. It shows the uplift deformation in the north of the co-seismic deformation field with a maximum deformation rate of 20mm/yr, while small-scale subsidence deformation was observed near Kathmandu in the south. However, the uplift deformation located on the MFT (27 • N) is much weaker than that observed via the ALOS data ( Figure 5).

Kinematic Afterslip Model
Kinematic afterslip is widely used in the interpretation of early post-seismic deformation. The dynamic source of the afterslip is the Coulomb disturbance caused by the co-seismic rupture on the fault plane, and the range of the afterslip is generally distributed around the co-seismic slip zone and does not extend far away [31]. We use the SDM (steepest descent method, [32]) to invert the kinematic afterslip based on the GPS and InSAR data to explore the temporal and spatial distributions of the afterslip. We adopt the co-seismic fault model [33], and extend the strike and dip directions to ensure coverage of the afterslip area. Moreover, we resampled the deformation field using the quadtree decomposition method and laid-off the southern deformation before inversion. According to our available post-seismic deformation observation results from GPS, ALOS-2 and Sentinel-1, we have carried out kinematic afterslip inversion of seven groups of dataset ( Figure 11). If only GPS data are set as constraints, the afterslip is mainly distributed in the downdip direction of the fault model and the maximum slip is~0.3 m in about 1yr (Figure 11a), the range of the afterslip extends northwards to 29 • N, and minor slip is observed on the updip zone. We jointly used GPS and ALOS-2 data for 1yr after the main earthquake, we analyzed the weight ratio between the two datasets based on the method of Wang & Fialko [34], which looks for the optimal value of α such that misfits for each data set are within 30% of the respective minimum values. We chose the weight value of 0.7 for InSAR/GPS, the kinematic afterslip inverted by the two datasets is shown in Figure 11b. The characteristics of the joint afterslip distribution are approximately consistent with the GPS results, only the afterslip magnitude in the north is slightly larger. The afterslips constrained by ALOS-2 data at different periods are distributed in both updip and downdip areas (Figure 11c-e). However, the magnitude of afterslip varies greatly in different time periods after the earthquake. The afterslips of 4.5 yr (Figure 11c) and 2.5 yr (Figure 11d) are similar, for example, the maximum afterslip is 0.21 m in the downdip area located at a depth of 19 km and 0.1 m in the updip area located at a depth of 10 km (Figure 11c), while the afterslip of 2.5-4.5 yr after the earthquake was much weaker (Figure 11e), the maximum afterslip in this stage was only 0.07 m in the downdip area. The characteristics of kinematic afterslip constrained by both ascending and descending Sentinel-1 post-seismic deformation are similar to those inverted by ALOS-2 only (Figure 11f,g). Furthermore, the kinematic afterslip inverted by the two kinds of InSAR data for the downdip area is significantly greater than that of updip angle direction, this may indicate that the afterslip mainly occurred in the deep brittle-ductile transition zone below the co-seismic rupture area. Figure 11e shows the afterslip inverted from Sentinel-1 post-seismic deformation obtained by the PS-InSAR method in 4.5 yr, the maximum slip in the downdip area is~0.3 m. Figure 11f shows the kinematic afterslip inverted from Sentinel-1 post-seismic deformation data obtained by Stacking method, with a maximum slip of 0.35 m in the downdip and 0.2 m in the updip areas, respectively. The range of the distribution of afterslip is roughly as same as for the ALOS-2 results, only in the north is the scope of afterslip slightly smaller than that of the ALOS-2 constrained values. Therefore, the results of the inversion suggest that the afterslip does exist in the updip area, and that the updip afterslip attenuation was fast in the first 2.5 yr and slow in the later period, while downdip afterslip continued during the observation period of 4.5 yr. The above data simulation results are shown in Figure S4, all simulation data match the observation data well.

Stress-Driven Afterslip
Although the kinematic afterslip results can fit the post-seismic deformation results well, they cannot truly reflect the physical mechanism of slip. We used a stress-driven afterslip model to predict the future post-seismic deformation and the afterslip on the fault plane. The open-source program RELAX was used to predict the post-seismic deformation [35]. The velocity-strengthening part of the fault can be expressed as Formula (4): where ∆τ is the co-seismic shear stress change, v 0 is the reference slip rate of the timescale afterslip, a represents the friction parameter of the rock, and σ is the effective normal stress on the fault. The slip rate is controlled by the local stress drop ∆τ. We treat the receiver fault the same as the kinematic afterslip fault model. Figure S5 shows the simulated GPS velocity field, when v 0 = 1 m/yr, a σ= 8 MPa. The simulated results overfit the observation data in the updip zone, while it is much smaller than that observed in the downdip zone, and the updip slip momentum is much larger than the downdip slip momentum. The stress-driven afterslip model allows creep to occur in the updip zone and produces large post-seismic deformation, while both GPS and InSAR data observe only small displacements in the updip zone. Therefore, we assume the updip zone locked after the Mw 7.8 earthquake.
Assuming that y 0 represents the distance from the surface trace of the MFT to the projection of the surface trace of the locked position, we searched for the optimal values of y 0 based on the RMS between the forward modeling data and the GPS time series data in Formula (5): where (D) ij 0bs and (D) ij mod represent the observed and simulated values, respectively, N and E represent the north and east deformation components, i represents the GPS station, j represents the time interval, m represents the number of time separations and n represents the total number of GPS stations. We obtained the optimal value of y 0 as 80 km ( Figure S6). We searched for the optimal values of v 0 and aσ in the ranges of 4-12 m/yr and 4-10 MPa, respectively. Figure 12 shows the optimal value of v 0 = 7 m/yr and aσ = 6 MPa when RMS is smallest. The forward model GPS velocity for the first year post-earthquake is shown in Figure 13, and the simulated data show a good agreement with the observations. The stress-driven afterslip is mainly distributed within the depth range of 15-25 km with a maximum value of~0.35 m.

Discussion
In order to further explore the temporal characteristics of afterslip and its released moment, based on the kinematic afterslip evolution model constrained by ALOS-2 postseismic deformation, we calculated the seismic moment released by the afterslip in three different time intervals, these are the 4.5 yr, 2.5 yr and 2.5-4.5 yr periods after the Nepal Mw 7.8 earthquake. The moment released in 4.5yr based on the kinematic afterslip is 1.5174 × 10 20 Nm, which is about~21.2% of co-seismic moment released by the mainshock (~7.14×10 20 Nm). The moment released in 2.5 yr is~1.3134 × 10 20 Nm, accounting for 18.4% of the co-seismic moment release, while the moment released from 2.5-4.5 yr only accounts for~2.8% of the co-seismic moment release. Results reveal that the afterslip mainly occurred in the early period, especially in the first 2.5 yr after the mainshock. Previous studies have also shown that energy released during the first 24 hrs of postseismic afterslip was comparable to 30% of the corresponding energy released during co-seismic slip events [36], and similar examples exist in other studies [37][38][39]. The first observation date of the ALOS-2 data used in this study is 17 May 2015 which was 22 days after the mainshock. While the GPS dataset observation was following the mainshock, so the afterslip magnitude constrained by GPS data in 1yr is approximately equal to that from ALOS-2 data in 4.5 yr. The afterslip takes effect mainly within 2.5 yr after the mainshock. The above observation and inversion based on ALOS-2 and Sentinel-1 data show that the post-seismic deformation exists in the updip area, although it is smaller than that in the downdip region due to the different frictional properties of the updip and downdip areas. The maximum cumulative displacement in the updip area was~25 mm and mainly occurred within 2.5 yr. The kinematic afterslip was~0.2 m in 2.5 yr and only 0.04 m in 2.5-4.5 yr, and the energy released by afterslip & aftershocks in updip area was small, and was possibly relocked soon after. Co-seismic slips of up to 8-10 m can occur with average recurrence times of~600 yr [40], partial rupture accounts for only 4-6 m of the slip for the Nepal earthquake, which is suggestive of heterogeneous frictional properties, so the risk of a future earthquake in the updip area is still high. The updip area is locked based on the results of the stress-driven afterslip model, while the kinematic afterslip results show that afterslip exists in the updip area, which suggests there are some patches with different properties in the updip region. In fact, the post-seismic deformation we observed in the updip region may include the contribution of aftershocks, not just the aseismic creep. From the distribution of the aftershocks, we can see that most aftershocks, including several Mw > 5 events exist in the co-seismic fracture zone and its upwardly dipping side within 45 days after the Mw 7.8 earthquake (Figure 11b), while there are few aftershocks in the deep part of the co-seismic rupture zone, which may indicate that the deep region is close to the ductile shear zone where aseismic slip is dominant. The aftershock distribution area is consistent with the kinematic afterslip area on the updip side. The Coulomb stress change we calculated with the co-seismic slip distribution shows a stress loading effect in both the updip and downdip areas ( Figure S7). Therefore, we speculate from all above analysis that the afterslip in the downdip area is caused by the stress-driven afterslip, while the afterslip in the updip area is mainly related to brittle fracture caused by aftershocks.
The long-term inter-seismic stress accumulation and release in the co-seismic period, causing permanent deformation from deep structures to the surface is the long-term pattern of action in the Himalayan orogenic belt [6,41,42]. The uplift of the Tibet Plateau is located on the mid-crustal ramp [5], and the hypothesis of Elliot has shown that the Nepal earthquake will promote the building of topography at the front of the mountain chain. Post-seismic deformation has leads to uplift deformation in the north of the co-seismic deformation field, which conforms to the hypothesis. Laboratory experiments have shown that an average recurrence cycle of~600 yr will produce a slip of~8-10 m, the Nepal earthquake was only~4-6 m, indicating the heterogeneous frictional properties of the region [40]. The width of the Himalayan arc determines the recurrence cycle [43]. Based on the paleo-earthquake and geodetically inferred moment deficit rate, the possibility of a millenary Mw > 9 earthquake is high [44]. When accounting for the physical fault slip and the frictional barriers of fault dynamic rupture, the MHT may induce an Mw > 8 earthquake every 200 years [45]. For velocity-strengthening bodies, the post-seismic effects manifested in the form of aftershocks and afterslip. Several previous studies have shown that the shallow part of the MHT is fully locked [29,41]. The relocation of aftershocks have shown that an aftershock swarm existed in the shallow part of the region, which is indicative of brittle fractures in the shallow part that are concentrated in one area, where no GPS data exists for the post-seismic period. Our observational results and the relocation of aftershocks show partial rupture in the shallow part of the region, which may have been fueled by a large earthquake.

Conclusions
In this study, we obtain 4.5 yr post-seismic deformation from ALOS-2 and Sentinel-1 data after the 2015 mainshock. The post-seismic deformation of the ALOS-2 data is mainly distributed in four areas. The direction of the post-seismic deformation is opposite to the co-seismic deformation near the co-seismic deformation field, and another two deformation areas are located in the region of 29 • N and near the MFT. The Sentinel-1 post-seismic deformation is mainly distributed in the downdip area with an average deformation rate of~20 mm/yr. The kinematic afterslip shows a slip distribution in both the updip and downdip areas, with the maximum slip being~0. 35  The afterslip in the updip area is mainly caused by brittle fracture due to the aftershocks.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/rs14020306/s1, Figure S1. The relationship of elevation and deformation; The blue dot is before linear correction; The red dot is after linear correction, The orange line is the trend of linear fitting before correction, Figure S2. The contribution of the GPS sites; the color background is the InSAR deformation field in 4.5 yr; the red triangular arrows is the GPS velocity in horizontal (a) and vertical dimension (b), Figure S3. The optinal weight of GPS/InSAR. The red star represents the optimal value, Figure S4. Simulated and residual maps of ALOS-2 and Sentinel-1A data.  Figure S5. Stress-driven afterslip mdel, the updip area not locked, Figure S6. The optinal value of y 0 , Figure S7. The coulomb stress change of the mainshock; the purple contour line represents the stress-driven afterslip; the black contour line represents the coseismic afterslip; the red contour line represents the coseismic deformation; the red star represents the epicenter.