Deriving Centimeter-Level Coseismic Deformation and Fault Geometries of Small-To-Moderate Earthquakes From Time-Series Sentinel-1 SAR Images

Small-to-moderate earthquakes (e.g. ≤Mw5.5) occur much more frequently than large ones (e.g. >Mw6.0), yet are difficult to study with InSAR due to their weak surface deformation that are severely contaminated by atmospheric delays. Here we propose a stacking method using time-series SAR images that can effectively suppress atmospheric phase screens and extract weak coseismic deformation in centimeter to sub-centimeter level. Using this method, we successfully derive coseismic surface deformations for three small-to-moderate (Mw∼5) earthquakes in Tibet Plateau and Tienshan region from time-series Sentinel-1 SAR images, with peak line-of-sight deformation ranging from 5–6 mm to 13 mm. We also propose a strategy to downsample interferograms with weak deformation signal based on quadtree mesh obtained from preliminary slip model. With the downsampled datasets, we invert for the centroid locations, fault geometries and slips of these events. Our results demonstrate the potential of using time-series InSAR images to enrich earthquake catalog with geodetic observations for further study of earthquake cycle and active tectonics.


INTRODUCTION
Surface deformations produced by shallow earthquakes have been widely studied by Interferometric Synthetic Aperture Radar (InSAR) since the 1992 Landers earthquake (Massonnet et al., 1993). In complementary with seismological observations, InSAR images provide static, high spatialresolution, near-field surface deformation measurements for earthquakes, allowing for better resolving the absolute location and depth of the slip distribution (e.g., Ferreira et al., 2011;Weston et al., 2011). Thanks to the continuously improving data quality provided by spaceborne SAR satellites, more and more studies have been conducted to earthquakes that produce measurable surface deformations (e.g., Elliott et al., 2016;Merryman Boncori, 2019), which improves our understanding of earthquake cycle and tectonic evolution. However, earthquakes that have been studied by InSAR observations mostly have relatively large magnitudes (e.g. >M w 6.0), associated with decimeters to meters of surface deformation. Detailed source properties of smaller earthquakes (e.g. ≤M w 5.5) with coseismic deformation in the level of centimeter or less are in general less studied by InSAR data, despite that they occur much more frequently than larger ones. Some of these small-to-moderate events can cause substantial damages to the infrastructures if occur at very shallow depths (e.g. Wei et al., 2015;Qian et al., 2019). Small-to-moderate earthquakes can also reveal clues about the local fault systems and stress field as they often rupture on unmapped faults (e.g., Xu et al., 2015).
An important reason for the scarce studies on small-tomoderate earthquakes by InSAR data is that, SAR interferograms are often severely contaminated by atmospheric delays, bringing difficulties to extract reliable deformation signal from a single interferogram. Time-series InSAR analysis methods can effectively suppress the atmospheric delays, allowing for extracting millimeter-level deformation along the radar's lineof-sight (LOS) direction. Such techniques have been developed since 2000s, such as PS-InSAR (e.g., Ferretti et al., 2000), SBAS (e.g., Berardino et al., 2002), and Stacking (e.g., Tymofyeyeva and Fialko, 2015). These techniques, however, are generally used for measuring long-term, continuous ground deformation and less applied to coseismic deformation studies. Some efforts have been made aiming at improving the signal-to-noise ratio (SNR) of coseismic deformation measurements from time-series InSAR images. For instance, Xu et al. (2015) stack the coseismic interferograms to reduce atmospheric errors after estimating and removing linearly elevation-dependent atmospheric signals, which is demonstrated through inverting source parameters for the 2004 M w 5.1 Tabuk earthquake in Saudi Arabia. Grandin et al. (2017) implement a step-function model with the sum of a constant term to extract the subtle signal of coseismic displacement from multi-interferograms spanning the 2018 M w 5.8 Pawnee earthquake in Oklahoma. Although these earthquakes are relatively small, they still caused peak coseismic deformations at the level of a few centimeters.
Here, our effort focuses on pushing the capability of InSAR coseismic signal detection to centimeter and sub-centimeter level. To this end, we select the Tibet Plateau and its surrounding areas as the target region, with two reasons. Firstly, Tibet Plateau is one of the most tectonically active regions that has hosted many moderate to large earthquakes (e.g. >M w 6.0) at shallow depths, many of these recent events have been studied by InSAR data (e.g. Elliott et al., 2010). However, much more smaller earthquakes (e.g. ≤M w 5.5) have been left unstudied. Secondly, Tibet Plateau is relatively dry and less vegetated, therefore coherence can be maintained in a higher level, and the atmospheric impact is much weaker than tropical or sub-tropical regions.
We use SAR data acquired by the C-band Sentinel-1A and 1B satellites launched by the European Space Agency (ESA), which began to provide large-coverage SAR images with 250 km-width frame since April 2014 and April 2016, respectively. The satellites image the Earth surface regularly with short revisit period of 12 days for single satellite and 6 days for twin constellation (Yague-Martinez et al., 2016). The high-temporal sampling rate of Sentinel-1 data produce sufficient images for stacking, which is essential to reduce the local atmospheric turbulence that is difficult to be handled by numerical weather models (Tymofyeyeva and Fialko, 2015). Therefore, the large coverage and short revisit period of Sentinel-1 SAR images provide great potential to extract very weak coseismic deformation for small-tomoderate earthquakes.
In this work, we successfully retrieve coseismic deformations for three earthquakes with magnitudes from M w 5.2 to 5.5 in Tibet region by applying a stacking method using time-series SAR images (Figure 1). These three representative events have thrust, strike-slip and normal faulting mechanisms, respectively. To prepare data for slip inversion by these weak coseismic deformation signals, we propose a new downsampling strategy based on initial model-derived surface deformation. Using these downsampled datasets, we invert for their locations, depths, fault geometries and slip vectors. By studying such small-to-moderate events, we show the potential to enrich earthquake catalog with their source imaged by geodetic observations and discuss its potential use in earthquake cycle studies.

Time-Series InSAR Data Processing
To derive centimeter to sub-centimeter coseismic deformation, we use time-series Sentinel-1 SAR images acquired before and after an earthquake. The basic SAR interferometry is achieved using the burst-based processing chain implemented in the Sentinel-1 Interferometry Processor (http://sarimggeodesy. github.io/software) (Jiang et al., 2017). The Digital Elevation Model (DEM) from Shuttle Radar Topography Mission (SRTM) is used for coregistration and topography correction. The burst interferograms are mosaicked into one ∼250-km wide interferogram after coregistration. The phase discontinuity between adjacent bursts is estimated from burst-overlap interferometry and satellite parameters, and is corrected before mosaicking (Yague-Martinez et al., 2016). Then the mosaicked interferogram is mutilooked with factors of 4 and 16 in azimuth and range directions, respectively, and is then filtered using the Goldstein method to further improve the SNR (Goldstein and Werner, 1998).
For the proposed stacking method, we select the first SAR image acquired just before the earthquake as reference image and co-register all the other images with it using the above mentioned interferometry processing. We only keep data points with high mean coherence values (i.e., >0.3). Then we unwrap the phases on these selected points based on the Statistical-cost, Network-flow Algorithm for Phase Unwrapping (SNAPHU) (Chen and Zebker, 2000) and manually correct unwrapping errors. After unwrapping, a linear ramp and a linear relation with the elevation are estimated and removed from each interferogram to account for possible long-wavelength and topography-related atmospheric delays (Ding et al., 2008;Jiang et al., 2014). We also discard interferograms with extreme atmospheric turbulence that is difficult to eliminate, Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 2 or highly decorrelated interferograms that cannot be successfully unwrapped.
If N + M+1 (N+1 before and M after the earthquake) SAR images were acquired around the time of an earthquake, we can derive N interferograms that do not include coseismic deformation and M interferograms that include coseismic deformation using the same reference image. The phases in these unwrapped interferograms could be written as: where k indicates the reference SAR image acquired just before the earthquake, ∅ k,j and ∅ i,k represent the unwrapped phases with and without coseismic signals, respectively. ∅ atm is the atmospheric contribution, ∅ def ,j is the coseismic deformation contribution and ∅ noise denotes other contributions, such as unwrapping errors, decorrelation errors and etc. Because the long-wavelength and topography-related atmospheric signal have been estimated and removed in each interferogram, the remaining phases should be primarily contributed by temporally unrelated atmospheric turbulence for ∅ i,k or the mixture of atmospheric turbulence and coseismic deformation for ∅ k,j . We stack the unwrapped interferograms with and without the coseismic deformation separately as: Through the stacking, the residual noises ∅ noise and the temporally unrelated atmosphere turbulence ∅ atm,n and ∅ atm,m have been largely suppressed by the temporally averaging. The stacking result of the N interferograms without coseismic signals, i.e., ∅ non mainly exhibits the atmospheric delays introduced by the reference image (with the sign flipped as ∅ atm,k ) (Tymofyeyeva and Fialko, 2015). While the stacking result of the M interferograms with deformation signal from the earthquake, i.e., ∅ non+def , mainly include the atmospheric delays introduced by the reference image and the coseismic  (Xu et al., 2016;Wei et al., 2018). Insets show detailed background seismicity with magnitude larger than M2.5. Within each inset, the colored beachballs represent the locations and mechanisms from USGS (red, star if no mechanism is given) and Global Centroid Moment Tensor (GCMT) (green) (Ekstrom et al., 2012), respectively. Yellow stars are earthquake epicenter locations from China Earthquake Networks Center (CENC).
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 deformation. Then, we can get the estimation of coseismic deformation ∅ def by summing the two stacking results (∅ non + ∅ non+def ) with the atmospheric delays of the reference image being canceled out. For small-to-moderate earthquakes (e.g. M w ≤ 5.5), the possible post-seismic deformation is considered to be negligible due to their small coseismic slip. Therefore, we take this estimation (∅ non + ∅ non+def ) as the final measurement of coseismic deformation.
To demonstrate the capability of our method in extracting small coseismic deformation signal, we show a representative application to the 2017/09/16 M w 5.5 Kuche earthquake ( Figure 2). For this earthquake, we collect 37 and 36 Sentinel-1 SAR images from the ascending track AT12 and descending track DT165, respectively. SAR images acquired on 2017/09/08 (AT12) and 2017/09/07 (DT165) are selected as reference images because they are the images acquired just before the earthquake. We stack the unwrapped phases on the selected data points (Supplementary Figures S1 and S2) for InSAR interferograms with and without deformation signals in both tracks ( Figure 2A and D, B and E). In the stacked InSAR interferograms that cover the earthquake, atmospheric turbulence from the reference images are dominating the data where the coseismic deformation signals are hardly visible due to their weaker amplitudes ( Figures 2B,E). As expected, clear coseismic deformation signals emerge when we sum the stacked results with and without the earthquake signals, as this procedure effectively cancels out the atmospheric turbulence phases ( Figures 2C,F). Similar enhancements to the coseismic signals in both ascending and descending tracks verify the reliability of our method. For this case, the earthquake locations from USGS, GCMT and CENC are all biased to the northwest of the deformation center, yet within 10 km to the peak of deformation.
To quantitatively evaluate the improvement, we calculate the standard deviations of unwrapped phases in the area out of coseismic deformation for both ascending and descending tracks, before and after applying the proposed stacking method. The values decrease from 1.27 radians and 1.31 radians to 0.27 radians and 0.28 radians, equal to ∼0.12 cm, both achieving a noise reduction of ∼80% (Supplementary Figure S7). This quantitative analysis clearly demonstrates the capability of our method in reducing the influence of atmospheric turbulence and improving SNR. Note that considering the complicated temporal decorrelation mechanism of InSAR, the quantitative evaluation may vary for different earthquakes.

Downsampling Strategy
To reduce the number of data points for slip inversions, we need downsample the original InSAR data. Quadtree partitioning is one of the widely used spatial variant downsampling methods (Jonsson et al., 2002), which is based on calculating variance of the measurements to split original interferograms into different levels of quadrants. If the variance, which represents the gradient of deformation within the quadrant, excesses a pre-determined threshold, the corresponding quadrant is split into four equal sub-quadrants. This procedure is repeated until the variance within each quadrant reaches the pre-assigned minimum threshold or the quadtree level reaches the pre-assigned maximum threshold. The average or median value in each quadrant is calculated to represent this quadrant. Note that  ) and (F) are the final measurements of coseismic deformation signals derived by summing the two previous stacked results in each track. One phase cycle (from -π to π) represents ∼2.8 cm of LOS displacement. Blue means moving towards the satellite. Stars represent earthquake locations from USGS (epicenter, red), GCMT (epicenter, green) and CENC (centroid, yellow), respectively.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 this averaging or median-filtering process suppresses the outliers (mostly unwrapping errors) and therefore improves the SNR in the downsampled data. If the SNR is high, the variance of data points within each quadrant represents the level of deformation gradient of that quadrant. As a consequence, the regions with the densest quadrants are associated with the largest deformation gradients in an unwrapped interferogram. However, for small-to-moderate earthquakes, the residual atmospheric phases and decorrelation noises may have comparable or even larger amplitude relative to the earthquake signal (e.g., Figures 3B,E,H). Therefore, regular quadtree partitioning may lead to many small quadrants, which bring no extra or even negative contribution to constrain the earthquake source parameters.
To overcome this issue, we propose a refinement to the quadtree downsampling strategy. Firstly, we obtain an initial slip model using the standard-quadtree downsampled data points and set the initial values of the source parameters as that from available catalogs (e.g., focal mechanism and fault geometry). We then derive a new quadtree downsampling mesh on the noise free synthetic LOS displacement field produced by this initial slip model. Finally, we apply this new mesh to real InSAR data for downsampling, which effectively reduces the number of data points ( Figures 3C,F,I). When the SNR of an interferogram is low, we can also apply the new downsampling process to each interferogram before stacking (Supplementary Figures S3 and S4). Coseismic deformation can be confirmed if their signals continually appear in the time-series downsampled data points after the earthquake. Then we apply the proposed stacking method on these points to obtain the final deformation dataset for inversion.
To quantitatively evaluate the improvement, we compare the downsampled results using standard-quadtree method and the proposed method, and count the numbers of data points inside and outside the deformed area (Supplementary Figure S8). By using standard-quadtree downsampling method, although the number of FIGURE 3 | Single coseismic interferograms and final deformations derived from the proposed method for the M w 5.5 Kuche, M w 5.2 Nierong, and M w 5.2 Zhongba earthquakes from top to bottom. From left to right, column one shows the coseismic interferograms, column two shows the coseismic deformation signals extracted by the proposed stacking method, column three shows the downsampled data points, with red indicating moving towards and blue indicating moving away from the satellites along the radar's LOS direction. Stars represent the earthquake locations from USGS (epicenter, red), GCMT (epicenter, green) and CENC (centroid, yellow), respectively.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 5 data points have been reduced from more than one million to about 1,500, many unwanted data points retained due to the high gradients caused by phase noises. After applying the proposed downsampling method, the numbers of data points outside the deformed area have reduce from 1,385 to 260 for the ascending case and from 1,302 to 254 for the descending case, while the numbers of data point inside the deformed area keep nearly the same (217 vs. 199 and 193 vs. 222). The percentages of data points inside the deformed area increase by about 3.2 and 3.6 times for ascending and descending tracks, respectively. Consequently, the InSAR data downsampled by the proposed method show relatively better constraints on fault parameters (e.g. fault width and depth) compared with the result using standard-quadtree downsampling method (Supplementary Figures S9 vs. S7). The misfit in the inversion is also 33% smaller (Supplementary Table S1). This quantitative analysis clearly demonstrates the advantage of our downsampling strategy for studying small-to-moderate earthquakes with weak coseismic deformation.

Fault Geometry Inversion
With the downsampled datasets, we can invert for the source parameters. The earthquake is modeled as a rectangle fault plane with nine unknown parameters describing fault size (width and length), geometry (strike and dip angles), three-dimension location (latitude, longitude and depth) and uniform slip (amplitude and direction) (Okada, 1985). The static Green's functions are calculated from a homogeneous elastic half-space with Poisson ratio of 0.25. During the inversion of the nine unknown parameters, we use the Geodetic Bayesian Inversion Software (GBIS) to search for the optimal fault geometry and slip vector (Bagnardi and Hooper, 2018).
In the Bayesian framework implemented in GBIS, the posterior probability density function (PDF) is calculated based on the residuals between the data and the model prediction. For each iteration, if the likelihood for the new set of model parameters is greater than previous one, the trial model values are retained. Otherwise the previous set of model parameters is retained. The Markov chain Monte Carlo sampling method is applied during the inversion, incorporating the Metropolis-Hastings algorithm. The final optimal values of model parameters are from the last retained set of model parameters with the maximum posteriori probability after pre-assigned iterations (e.g. 500,000 iterations) (Bagnardi and Hooper, 2018).
For those small-to-moderate earthquakes (e.g. ≤M w 5.5) with centimeter-level surface deformation, the searching result for the optimal parameters, particularly, the strike and dip angles, may not be as robust as that for larger (e.g. >M w 6.0) earthquakes (e.g. Wang et al., 2018). Therefore, we fix the initial strike and/or dip angles based on the focal mechanisms from reputed earthquake catalogs (e.g. USGS or GCMT) and perform the source inversion to obtain the optimal location, depth and slip vector.

RESULTS
We apply the above-mentioned procedure to three earthquakes occurred in Tibet and Tienshan region and present the results in the following sections. These three representative earthquakes have thrust, strike-slip and normal faulting mechanisms, respectively (Figures 1 and 3). According to global catalogs (e.g. USGS, GCMT), they are: the September 16, 2017 M w 5.5 Kuche thrust-faulting earthquake occurred in the south margin of Tienshan; the December 4, 2016 M w 5.2 Nierong strike-slip earthquake occurred in central Tibet Plateau, and the February 1, 2017 M w 5.2 Zhongba normal-faulting earthquake occurred in southern Tibet Plateau. For these three small-to-moderate earthquakes, time-series Sentinel-1 SAR images in both ascending and descending tracks are collected to provide better constraint on the fault geometries and slip vectors (see supplementary figures for time-series results from all tracks).

M w 5.5 Kuche Earthquake
The 2017 M w 5.5 Kuche earthquake (Figure 1, Figures 3A-C and Figure 4A) occurred in the southern Tienshan thrust-andfolding zone, resulted from the converging with the Tarim Basin to the south (Avouac et al., 1993). Seismological solutions from USGS and GCMT all suggest thrust faulting mechanisms, while the earthquake either occurred on a shallowly northward dipping fault (dip 20°, USGS or 28°, GCMT), or a steeply northward dipping fault (dip 70°, USGS or 62°, GCMT).
For the surface deformation, although the coseismic interferogram by single pair of SAR images ( Figure 3A) in descending track DT165 is highly coherent without observable phase noises, the atmospheric turbulence seems to be dominating in the interferogram and the deformation signals are hardly visible. It is clear that single interferogram cannot provide reliable measurement for the coseismic deformation for this event. Nevertheless, coseismic signals are much clearer after applying the proposed stacking method (Figures 2 and 3B). We convert the phases to displacement and apply the proposed downsampling strategy.
Besides the descending track DT165 shown in Figure 3C, we also derive the coseismic deformation from ascending track AT12 using 19 and 18 Sentinel-1 SAR images acquired before and after the earthquake, respectively (Supplementary Figure S2). We obtain clear coseismic signal in ascending track AT12 with maximum displacement of 5-6 mm along the LOS direction ( Figure 4B), similar to that in the descending track. To our best knowledge, this is the smallest coseismic deformation signal that is robustly observed with InSAR technique. Earthquakes with such weak coseismic deformation has been seldom studied by InSAR data.
The InSAR images from both tracks show a single deformation patch with dominant LOS displacement towards the satellite, suggesting a thrust faulting mechanism that is consistent with seismological focal mechanism. Although the slip model fits the InSAR observations very well ( Figures 4B-G), the dipping direction of the fault plane still can't be resolved due to very small fault dimension and relatively deep centroid depth (13.5 km) (Table 1, Supplementary Figure S7). Also because of deep centroid depth, the other source parameters are not as well constrained as the other shallower events (e.g. Supplementary Figures S7 vs. S8).
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 6 Based on the local tectonic setting, the earthquake is close to two active, thrust-fold faults in the south margin of Tienshan region ( Figure 4A). Considering the pervasive, north-dipping thrust faults in the region, we prefer the northward, shallowly dipping fault plane for the Kuche earthquake. We therefore fix the dipping angle as 28°according to the GCMT focal mechanism for the slip inversion (Table 1, Supplementary Figure S7). The best centroid depth of the fault plane is 13.5 km, slightly shallower than that from USGS and GCMT catalogs but deeper than that from CENC catalog. The strike and rake angles derived from our inversion exhibit very similar values to the solutions from GCMT catalog ( Table 1).

M w 5.2 Nierong Earthquake
The 2016 M w 5.2 Nierong earthquake (Figure 1, Figures 3D-F and Figure 5A) occurred in the central Tibet. USGS catalog reports a strike-slip focal mechanism, suggesting either a left lateral event on a southwest-northeast striking fault or a right-lateral event on a southeast-northwest striking fault ( Table 2). The coseismic deformation can be hardly recognized from the single interferogram, yet emerges after applying the proposed stacking method to the ascending track AT143 (Figures 3D,E). The stacking result for the Nierong earthquake has a lower SNR compared with  the Kuche case, mainly due to the larger decorrelation noises ( Figure 3E). After applying our proposed downsampling strategy, the downsampled data shows clearer coseismic signals with a maximum amplitude of 8-9 mm along the LOS direction. In the area that is further away from coseismic deformation, the outliers are sufficiently suppressed, although the amplitude in the data (mostly noise) is comparable with the coseismic deformation ( Figure 3F). Besides the data for the ascending track AT143 ( Figure 3F), timeseries Sentinel-1 SAR images from the descending track DT150 are also processed (Supplementary Figures S3, S4). The maximum LOS displacement from descending track is about 1 cm ( Figure 5E), slightly larger than that in the ascending data ( Figure 3F, Figure 5B).  We use the downsampled dataset for slip model inversion. Similar to what we have done for the Kuche earthquake, we assume uniform slip on a rectangle fault and search for optimal values of all the nine source parameters. While different with the Kuche earthquake, the gradient of observed coseismic deformations for the Nierong earthquake is relative larger to better constrain the coseismic slip model. Therefore, we performed geodetic inversion on both nodal planes and selected the fault plane with smaller misfits ( Table 2). As shown in Figures 5B,E, both ascending and descending tracks show typical deformation pattern for a strike-slip faulting mechanism, which is well reproduced by our preferred slip model ( Figures 5C,F) (Table 2). Our InSAR derived fault geometry is highly consistent with the solution from USGS catalog. But note that the centroid depth derived from InSAR inversion (4.8 km) agrees better with the CENC report (5 km), both are much shallower than the USGS catalog (10 km).
Compared with the Kuche earthquake, the Nierong earthquake is much shallower (4.8 km vs 13.5 km), hence the surface deformations have larger spatial gradient. The inverted source parameters therefore are better constraint with smaller uncertainties (Supplementary Figures S8 vs. S7).
The misfit in the southwest-northeast striking plane (strike angle 57°by InSAR) is slightly smaller than the southeastnorthwest striking plane (strike angle 148°by InSAR) ( Table 2). The fault plane in northeast-southwest striking direction is nearly perpendicular to the nearly west-east striking faults to the north ( Figure 5A). Such striking direction and focal mechanism may be explained as a transition fault connecting the normal fault to the south to the left-lateral fault to the north ( Figure 5A). Considering better data fitting and the tectonic background, we prefer the fault plane in southwest-northeast striking direction for the 2016 M w 5.2 Nierong earthquake.  Figure 6A) occurred in southern Tibet Plateau where a series of major NS-trending rifts accommodate the east-west extension of the plateau (Ryder et al., 2012). Seismological solution from GCMT catalog indicates a normal faulting focal mechanism with nodal planes dipping towards east or west (Figure 1). For the surface deformation, the coseismic signals can be recognized from the single interferogram in the descending track DT92 ( Figure 3G), but the signals caused by the temporally uncorrected atmospheric turbulences are also visible with similar magnitude of the coseismic deformation. After applying the proposed stacking method, the atmospheric phase screens have been largely reduced, resulting in very clear coseismic signals ( Figure 3H). We then convert the phases to displacements and downsample the observations, which have a maximum deformation of 1.3 cm along the LOS direction ( Figure 3I). Besides the observations in descending track DT92 ( Figure 3D), coseismic deformation is also measured in ascending track AT158 with the maximum LOS deformation of 1.2 cm (Supplementary Figures S5 and S6, Figure 6E). Although the coseismic deformation is larger than the Kuche and Nierong events, the deformation gradients are still too small to resolve dipping direction of the ruptured fault plane. As Figure 6A shows, the earthquake is located ∼5 km to the west of a graben, which is bounded by two normal faults. We therefore prefer a fault plane that dips towards east, i.e. the graben, for the Zhongba earthquake and fix the dipping angle as 43°according to the GCMT focal mechanism for the slip inversion ( Table 3). The preferred slip model fits the InSAR observations very well ( Figures 6B-G) and the strike and dip angles are very similar to the solution from GCMT catalog ( Table 3).

DISCUSSION AND CONCLUSION
InSAR technique can provide high spatial-resolution surface deformation measurements in sub-centimeter accuracy and has been widely applied in earthquake studies. Nevertheless, for small-to-moderate earthquakes, InSAR observations are restricted by the weak and smooth surface deformation. Thanks to the Sentinel-1A/B satellites that routinely image the Earth surface with a large scan-width of 250 km and shortrepeating time, more earthquakes with centimeter-level coseismic deformation can be captured by the proposed stacking method using time-series SAR images. In this study, the maximum LOS deformation amplitude is only ∼6 mm for the 2017 M w 5.5 Kuche earthquake, which is well less than previous earthquake studies using InSAR in the literature.
Note that our method requires large amount of SAR images to achieve the best performance in eliminating atmospheric turbulence. It may be difficult for previous SAR missions such as ERS and Envisat. Nevertheless, the Sentinel-1 mission is designed to monitor the active Plate boundary regions. At least 12-days repeating time can be guaranteed for tectonic active region, such as Tibet Plateau, allowing us to study such regions with high seismicity rate. Another limitation of our method is that we cannot avoid the leakage of post-seismic motion in our measurements. However, the postseismic motion for small-to-moderate earthquakes (e.g. ≤M w 5.5) should be small and have very limited influence on the centroid location and fault geometry inversion, as postseismic slips likely occur on or near the same fault that ruptured coseismically.
The exact number of images used in the proposed method depends on many factors, such as coherence, topography, and the level of atmospheric delays. For example, for the 2017 M w 6.5 Jiuzhaigou earthquake in Sichuan, decorrelation noise leads to no reliable coseismic measurement along the rupture for the interferogram with only 12-days time interval (Sun et al., 2018). While for the 2017 M w 5.5 Kuche earthquake in this study, even with eight-month time interval (Supplementary Figures S1 and S2), the interferogram exhibits good coherence. Therefore, situation varies much with earthquakes occurred at different places. It is difficult to directly define the specific number of images required for applying the method. Generally, for earthquake occurred in less vegetated area, we can use 10-30 SAR images to achieve a large noise reduction. In the highly vegetated area, we need carefully select interferograms with relative high coherences to avoid the effects of decorrelation noises, and therefore, resulted in less SAR images involved in our stacking method.
To further explore the detection threshold for small-tomoderate earthquakes by time-series InSAR images, we perform a large number of forward modeling for earthquakes with various magnitudes and depths. For each model with given magnitude, a square fault plane is assumed with its area calculated by a typical stress drop of 2.7 MPa (Kanamori and Anderson, 1975;Abercrombie, 1995) using a shear modulus of 30 GPa (Jonsson et al., 2002). Despite that different faulting mechanism may result in different deformation patterns and amplitudes, we take the focal mechanisms from the three events studied in this work as examples in the forward modeling and calculate the maximum surface deformations along the descending LOS direction. The maximum LOS deformations are plotted as a function of event magnitude and centroid depth for these three representative events with thrust-faulting, strike-slip and normal-faulting mechanisms, respectively ( Figure 7).
As expected, there is a clear trade-off between the earthquake magnitude and centroid depth. If the limitation of the maximum LOS deformation that the proposed stacking method can resolve is 5 mm, the detection threshold can then be determined (the black lines in Figure 7), where events beneath the threshold cannot be detected. Because of the large uncertainty of centroid depth determined by tele-seismic observations, such threshold can be generated based on magnitude and fault geometries from reputed earthquake catalogs before InSAR analysis, providing useful information to roughly determine whether the earthquake is suitable to be studied with InSAR. Given similar magnitudes, the resolvable depths associated with dip-slip earthquakes are relative deeper ( Figures 7A,C) than that from the strike-slip event ( Figure 7B). This is because InSAR is more sensitive to the vertical ground motion. For an earthquake with given magnitude and depth, normal-faulting or thrustfaulting events produce mainly vertical surface deformation, which is more possible to be detected by InSAR. For strikeslip event that produces mainly horizontal surface deformation, the situation is more complicated, as different strike angles will cause different LOS projections too. Nevertheless, referring to the detection threshold by InSAR, the proposed method allows us to study shallow dip-slip earthquakes (shallower than 5 km) with moment magnitudes less than M w 5.0. Therefore, we believe that more and more small-to-moderate earthquakes can be studied with InSAR constraints, significantly enriching earthquake catalogs with geodetically derived fault parameters.
Repeating earthquakes that occur on the identical fault segment are critical for earthquake cycle study, while large characteristically repeating earthquakes (e.g. >M w 7.0) are rare with long recurrence interval. Most well-documented repeaters are detected by seismic data with small magnitudes (Uchida and Bürgmann, 2019). Such earthquakes maybe in the detection range of our method if with proper magnitude and depth (refer to Figure 7). The large amount of archived Sentinel-1 data in global coverage provide the great opportunity to study these repeating earthquakes if their surface deformation is measurable. We believe that along with more and more small-to-moderate events that can be detected by the proposed method, we may find repeating earthquakes with InSAR measurements spanning the full earthquake cycle. For such case, both geodetic and seismic data can be used to reveal more details of their fault geometry and rupture process, casting new lights on earthquake cycle studies.
Our study shows that InSAR technique can provide high spatial-resolution surface deformation measurements at subcentimeter accuracy and is successfully applied to earthquake source studies. However, we realize that InSAR data alone may not be sufficient to resolve the ruptured fault plane for the three earthquakes we presented in this study. As Supplementary Figure  S13 shows, due to the deep depth (13.5 km) of the Kuche event, the coseismic deformations are very small (6-8 mm) and smooth, providing weak constraints on slip model. The Markov chains for the model parameters (except earthquake location) cannot converge well after 500,000 iterations. Consequently, the optimal solution (red lines in Supplementary Figure S10) can be off from the peak of the model-parameter distribution. While for the Zhongba earthquake, the centroid depth is very shallow (<5 km), coseismic observations of ground deformations can provide relatively strict constraints on the slip model, the Markov chains for the source parameters are well converged (Supplementary Figure S14) and the red lines of inverted parameters are near the peak of the model-parameter distributions (Supplementary Figure S12). Moreover, it is more difficult to resolve the rupture details when the centroid becomes deeper. Our preferences of the ruptured fault plane are primarily based on independent constraints (e.g. nearby fault mapping and background tectonics). Therefore, such geodetic data should be used together with seismological and geological observations to better constrain the source parameters, such as rupture processes, slip distribution and stress drop on the fault.
To conclude, in this study, we propose a stacking method using time-series SAR images to extract reliable sub-centimeter coseismic 2 Zhongba earthquake. The surface deformation is projected along the descending LOS direction, and plotted along with magnitudes and centroid depths. Black solid line represents 5 mm of LOS deformation, which is supposed to be the maximum coseismic deformation signals that can be extracted by the proposed stacking method.
Frontiers in Earth Science | www.frontiersin.org February 2021 | Volume 9 | Article 636398 deformation, reaching an unprecedented level than regular single interferometry for studying earthquakes. Our results imply that it is possible to significantly enrich earthquake catalogs with geodetic observations, thus improving our understanding of earthquake physics and active tectonics.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
TW initialed and conceived the study. HL conducted the InSAR data processing and coseismic slip model inversions. ML provided funding and supervising for the study. HL wrote the manuscript with contributions from TW and SW. All authors discuss and interpret the results and the manuscript.