Temporal and spatial variation of fault creep along the Xianshuihe fault from InSAR stacking

SUMMARY

The left-lateral Xianshuihe fault is a seismically active fault system located at the eastern boundary of the Tibetan Plateau. We analysed the Sentinel Interferometric Synthetic Aperture Radar (InSAR) data from 2014 to 2021 to study the temporal and spatial variation of fault creep along the Xianshuihe fault. We applied the InSAR stacking method and the coherence-based Small-BAseline Subset (SBAS) method to derive the Line-Of-Sight (LOS) velocity map and time-series from both the ascending and descending orbits. We studied both the secular and the time-dependent components of surface deformation from InSAR. We compare the InSAR-derived velocity maps with the Global Positioning System (GPS)-derived velocity field and found that these two independent measurements are consistent. A 200 km long creeping section is identified along the central segment of the Xianshuihe fault. The surface creep rate is measured to be ranging from 0 to 6 mm yr⁻¹. We combined the elastic dislocation model and the InSAR velocity maps to invert for the geodetic fault slip rate and the aseismic slip distribution in the upper crust. The secular fault creep model shows that most of the Xianshuihe fault is creeping at depth. The time-dependent fault creep model indicates that the maximum aseismic slip rate from Bamei to Kangding accelerated from 30 to 40 mm yr⁻¹ and then decayed to 5 mm yr⁻¹ from 2014 to 2021. The fully creeping segment of the Xianshuihe fault seems to become a partially locked segment in a short time period (a couple of years). We suspect that the acceleration of fault creep from 2017 to 2019 is linked to dynamic triggering by passing seismic waves or fluid migration. Finally, we compare the temporal variation of fault creep with previous studies and discuss the earthquake hazard implications.

1 INTRODUCTION

The left-lateral strike-slip Xianshuihe fault is a major intracontinental fault system located in southwest China. It plays a significant role in accommodating the post-collisional convergence between the India and the Eurasian plates (Tapponier & Molnar 1977; Deng et al. 1979). The India plate moves northward toward the Eurasian plate which caused extensive deformation of the Tibetan crust. A significant part of this deformation is accommodated by the eastward extrusion of the Tibetan crust along the Xianshuihe fault (Zhang 2013). The total length of the Xianshuihe fault reaches 350 km approximately. In general, the Xianshuihe fault strikes from northwest to southeast. To the northwest of the Xianshuihe fault is the Ganzi–Yushu fault and to the southeast of the Xianshuihe fault is the Daliangshan fault, Xiaojinhe fault and Zemuhe fault (Jiang et al. 2015a). It marks the tectonic boundary between the Bayan Har block to the northeast and the Qiangtang block to the southwest (Fig. 1).

Xianshuihe fault is one of the most seismically active fault systems in China. From historical earthquake record and palaeoseismological data, there are 9 M > 7 and 20 M > 6 earthquakes occurred along the Xianshuihe fault since 1700 AD, including the 1973 M_s 7.4 Luhuo earthquake, the 1981 M_s 6.8 Daofu earthquake and the most recent 2014 M_w 5.8 Kangding earthquake doublet. (Allen et al. 1991; Wen et al. 2008; Jiang et al. 2015b; Bai et al. 2018). According to Wen et al. (2008), the Xianshuihe fault in the eastern Tibetan Plateau has hosted large earthquakes (M > 7.5) in 1327, 1786, 1816, 1955 and 1973. Recently there are many small earthquakes occurred near the Xianshuihe fault (Fig. 1). On 2022 September 5, an M_w 6.6 earthquake occurred near the Xianshuihe fault in Luding county Sichuan Province. It is located south of Moxi and north of Shimian (Fig. 1). Focal mechanism solution from Global Centroid-Moment-Tensor (CMT) (Dziewonski et al. 1981; Ekström et al. 2012) shows that this earthquake is probably a left-lateral strike-slip event along the northwest–southeast trending Xianshuihe fault.

It has been found that surface fault creep is occurring along the Xianshuihe fault. From short-baseline geodetic survey data, Allen et al. (1991) found that the Xianshuihe fault is creeping at a rate of 6 mm yr⁻¹ at the Xialatuo site (Fig. 1). Zhang et al. (2018) reported a creep rate of 1.3–3.5 mm yr⁻¹ along the Xianshuihe fault between latitudes 30.5° and 31.6° using short-baseline and short levelling array data. The creep rate at the Xialatuo site is decreasing exponentially with time from 8 mm yr⁻¹ in 1976–1980 to 2 mm yr⁻¹ or less in 2008–2014 (Zhang et al. 2018). It is possible that this decay of fault creep is a result of stress relaxation following the 1973 M7.6 Luhuo earthquake (Li et al. 2021).

Nowadays GPS and InSAR are two commonly used space-geodetic methods in the studies of intersesismic deformation and fault creep. GPS provides point measurement of crustal deformation so it has limited spatial resolution near active fault strands. On the contrary, InSAR can measure surface deformation at a spatial resolution of tens of metres and with an accuracy of 1 mm yr⁻¹. Thus, InSAR is particularly useful to image subtle interseismic deformation and surface creep near active faults. Wang et al. (2009) first used European Remote sensing Satellite (ERS) and Environmental Satellite (Envisat) data from 1996 to 2008 to reveal the interseismic deformation along the northwest Xianshuihe fault. Due to limited data quality and dense vegetation, they did not find a clear signal of fault creep near the Xianshuihe fault. Ji et al. (2020) processed Sentinel-1 data from 2014 to 2018 and found fault creep in the northwest Xianshuihe fault though the fault creep variation is not well determined. Using InSAR observations from 2015 to 2017, Li et al. (2021) indicate that the Xianshuihe fault is creeping at a rate less than 2 mm yr⁻¹. Li & Bürgmann (2021) studied the spatial and temporal variations of fault creep along the Xianshuihe fault using the Sentinel-1 InSAR data from 2014 to 2019. They reported a surface creep rate ranging from 0 to 6 mm yr⁻¹ along the Xianshuihe fault. They also found that the fault creep near Kangding segment is transient and it is probably controlled by the afterslip triggered by the 2014 Kangding earthquake doublet. The temporal characteristic of the fault creep is complicated due to transient effects caused by nearby earthquakes. Recently, Huang et al. (2023) studied the Kangding segment of the Xianshuihe fault with Sentinel-1 InSAR data from 2014 to 2021 and found that the creep rate decays from 40 to 0 mm yr⁻¹ between 2014 and 2019.

Since GPS observation has limited spatial resolution over the Xianshuihe fault and the fault creep signal is usually a small-scale feature, it is essential to employ InSAR data to investigate the interseismic fault creep near active fault strand. To accurately estimate the earthquake hazard of the Xianshuihe fault and understand the fault creep mechanism, we used space geodetic technique InSAR and Sentinel-1 Terrain Observation with Progressive Scan (TOPS) observations to investigate the interseismic deformation and fault creep of the Xianshuihe fault (Fig. 1). Previous work has studied the crustal deformation near Xianshuihe fault extensively but their results do not necessarily agree with each other. There are mainly two caveats in previous studies. First, they did not distinguish the depth extent of fault creep. When reporting the fault creep rates it is necessary to separate the surface creep from the deep creep. Surface creep only extend to the upper few hundred metres of the crust and the deep creep extend to the entire upper crust (∼10 km depth). Secondly, some studies did not distinguish the temporal characteristics of fault creep. The InSAR data are processed in such a way that the temporal variation of fault creep cannot be easily identify. In this study, we propose a new InSAR stacking method to utilize the abundance of the Sentinel-1 InSAR data. By doing so we can reveal the temporal variation of fault creep in detail.

2 InSAR DATA PROCESSING

The Sentinel-1A satellite was launched by European Space Agency (ESA) in 2014 and Sentinel-1B was launched in 2016. The Sentinel-1 TOPS data processed in this study is from 2014 December to 2021 April. During this time period there is no significant earthquakes occurred along the Xianshuihe fault (Fig. 1). We used GMTSAR, an open source InSAR processing software based on Generic Mapping Tools (GMT), to perform data processing (Sandwell et al. 2011). The C-band Sentinel-1 TOPS data is downloaded from the Alaska Satellite Facility (ASF). The orbit information is downloaded from ESA. We used data from both the ascending and the descending orbits to separate the horizontal movements from vertical movements. There are totally 156 SAR images for the ascending orbit and 101 SAR images for the descending orbit (Fig. 2), covering a large area over the Xianshuihe fault at the southeastern boundary of the Tibetan Plateau. The digital elevation model used in topographic phase correction is from the Shuttle Radar Topography Mission with 90 m resolution. SNAPHU software is used in phase unwrapping (Chen & Zebker 2000). Later the InSAR stacking and the following data processing steps are done with GMT.

3 METHOD

3.1 InSAR stacking

We used the stacking method (Sandwell & Price 1998; Tong et al. 2013) to reduce the phase of the selected interferograms to the mean line-of-sight (LOS) velocity. The stacking method can effectively reduce the atmosphere noise and enhance the tectonic deformation signal. The formula we used in stacking is described in Sandwell & Price (1998) and applied in the study of the interseismic deformation along the San Andreas Fault (Tong et al. 2013). The interferograms are selected based on the following rules: the spatial baselines of an interferogram must be smaller than 100 m. The temporal separation of the interferograms has to be longer than 1 yr. Thus, we obtain many long time span interferograms capturing the tectonic movement of the Xianshuihe fault during the interseismic period (Fig. 2). This rule ensures us to have high-quality interferograms with minimum decorrelation. For the ascending orbit we made 1187 interferograms. For the descending orbit, we made 1027 interferograms from the SAR images. Next, we visually inspected the unwrapped phase of all the interferograms and selected data with good phase recovery. The interferograms with bad unwrapping results or lower coherence are discarded.

We selected 79 high-quality interferograms out of the 1187 interferograms for the ascending orbit and 56 high-quality interferograms out of the 1027 interferograms for the descending orbit. These long-term interferograms are used in InSAR stacking and they are able to recover the small signals of the interseismic deformation and fault creep (Fig. 2). Compared to previous work employing InSAR data (Ji et al. 2020; Li & Bürgmann 2021), this study has used data with a much longer time span to image the interseismic deformation of the Xianshuihe fault. After interferograms stacking we used the remove/restore approach to correct for the orbital error and the large-scale atmospheric error. GPS data are interpolated to assist removing a planar trend. The difference between our approach and the SUm-Remove-Filter-restore (SURF) method used in Tong et al. (2013) is that we did not apply a high-pass filter to mitigate atmospheric noise. To further reduce noise in the mean LOS velocity model, we averaged the LOS velocity using a low-pass spatial Gaussian filter with 2 km wavelength. Thus, the spatial resolution of our LOS velocity maps is about 2 km. The filtering step is done with the grdfilter function in GMT. Finally, we removed the mean value from the entire LOS velocity grid so the reference is the average of the LOS velocity map.

3.2 Multitemporal stacking

The InSAR stacking method can only yield an average rate of deformation and is effective in regions with constant strain rate. Because past researches have suggested that the fault creep along the Xianshuihe fault possibly decays with time (Huang et al. 2023), we need to employ a time-dependent approach to process InSAR data. We proposed a modified InSAR stacking method on different time intervals of the observation, taking advantage of the abundance of the Sentinel-1 SAR data. Six years of the SAR data are divided into three equally spaced groups: 2015–2017, 2017–2019 and 2019–2021. For each group we selected good two-years interferograms and stacked them to obtain an averaged deformation rate within that time period (Fig. 2). By comparing the results from the two-years stacking we will gain insights on the temporal behaviour of fault creep. From our experience, 10–20 interferograms with a two-years’ time span can beat down the atmospheric noise significantly and yield a satisfactory LOS velocity map. We name this stacking method the multitemporal stacking and it has the potential to be applied in other regions where surface deformation is non-secular.

For the ascending orbit we selected 10 interferograms for the 2015–2017 period, 18 interferograms for the 2017–2019 period and 28 interferograms for the 2019–2021 period. For the descending orbit, we selected 15 interferograms for the 2015–2017 period, 14 interferograms for the 2017–2019 period and 22 interferograms for the 2019–2021 period (Fig. 2). The number of interferograms used in multitemporal stacking changes because the atmospheric noise and the degree of vegetation decorrelation keep changing with time. The choice of good interferograms depends on experience and it is difficult to quantify.

3.3 Coherence-based SBAS

To validate the multitemporal stacking method, we attempted to use the coherence-based Small-BAseline Subset (SBAS) method (Berardino et al. 2002; Schmidt & Bürgmann 2003; Tong & Schmidt 2016) to derive InSAR time-series. In the past studies, the SBAS method has been successful in deriving time-dependent deformation near the Xianshuihe fault (Qiao & Zhou 2021). Similar to previous study we only selected SAR acquisitions in the dry season (November to next year's February) to improve the time-series result. The baseline of the interferograms has to be less than 50 m for both the ascending and descending orbits. The temporal interval of the interferograms has to be less than 365 d for ascending orbit and less than 700 d for descending orbit. The time window for descending orbit is longer because there is a gap in SAR data acquisition from 2017 to 2018. The selected interferograms are shown in Fig. 2. After we derived the LOS displacement time-series we divided the temporal interval into three periods (2015–2017, 2017–2019 and 2019–2021) and we calculated an average rate within those periods to compare with the results from multitemporal stacking.

We found that for certain periods there is an apparent correlation between the derived LOS velocity and the topography in this region. We corrected this topography correlated atmospheric noise by fitting and removing a linear function between the interferograms’ phase measurements and the elevation. The correction is performed on individual interferograms that is subsequently used in coherence-based SBAS and multitemporal stacking. After this correction we found the LOS velocity less noisy than before and the fault creep signal can be readily identified.

3.4 Combine LOS to estimate fault-parallel velocity

We combined the LOS velocity map from both the ascending and descending orbits to obtain fault-parallel and vertical velocities. We assumed that the fault-perpendicular velocity in this region is negligible and we calculated the fault-parallel and the vertical velocities using an averaged look angle of the satellite and an averaged fault strike angle. In this calculation, we only used the overlap region of the ascending and descending LOS velocity map. The LOS velocity along the ascending and descending orbits can be written as:

where |${V_{\mathrm{ asc}}}$| is the LOS velocity along the ascending orbit, |${V_{\mathrm{ des}}}$| is the LOS velocity along the descending orbit, |${n_e}$|⁠, |${n_n}$|⁠, |${n_V}$| are the unit look vector of the satellite along the ascending orbit, |$n{^{\prime}_e}$|⁠, |$n{^{\prime}_n}$|⁠, |$n{^{\prime}_V}$| are the unit look vector of the satellite along the descending orbit, |${V_e}$|⁠, |${V_n}$|⁠, |${V_V}$| are the east velocity, north velocity and vertical velocity at a certain point on the ground. Formulae (1) and (2) can be written in the form of the fault-parallel velocity assuming there is no fault-perpendicular movement of the ground:

where |${V_{\textrm{para}}}$| is the fault-parallel velocity and |$\theta $| is the average strike angle of the Xianshuihe fault. By combining eqs (3) and (4), we can estimate the fault-parallel velocity and the vertical velocity from the LOS velocities. We used these formulae when estimating both the secular and the time-dependent deformation in this study.

3.5 Estimate surface creep rate

We investigated the surface fault creep rate variation using the high-resolution interseismic velocity maps from Sentinel-1 InSAR. To estimate the surface creep rate of the Xianshuihe fault, we took 18 cross-sections of the fault-parallel velocity grid along the fault trace. Each cross-section is about 20 km in length and the spacings between cross-sections is about 10 km. We attempted three different cross-section length (1, 5 and 10 km on either side of the fault) when estimating the creep rate. Next, we fit a linear function to the sampled velocity on either side of the fault. The difference of the intercept of the linear function in the middle of the cross-section is determined to be the estimated surface creep rate. The standard deviation of the surface creep rate measurements is based on the goodness of the linear fitting. The details of this method are described in Tong et al. (2013). We estimated the surface creep rate using this method for both the secular and the time-dependent fault creep.

3.6 Construct fault creep model

The surface deformation data can be used to constrain the slip distribution of fault creep in the upper crust (Li & Bürgmann 2021). Here, we describe the method to construct the fault creep model. There are mainly two contribution factors to the surface deformation imaged by InSAR. The first-order contribution is due to the interesismic locking of the Xianshuihe fault. It can be modelled with a fault plane slipping from the locking depth to infinite depth in an elastic half-space (Okada 1985). The second-order contribution is due to fault creep in the upper crust. And it is modelled with a fault plane slipping from surface to the locking depth (Okada 1985). Because the first-order contribution is significantly larger than the second-order contribution we need to remove the effect from interseismic locking before modelling fault creep. Otherwise, the effect from interseismic locking will lead to unrealistically large fault creep estimation. We consider the residuals from the interseismic locking and use that to constrain the fault creep model for both the secular fault creep model and the time-dependent fault creep model.

We investigate the interseismic locking of the Xianshuihe fault with the LOS velocities from Sentinel-1 InSAR. When modelling interseismic deformation we consider the variation of fault strike along the Xianshuihe fault and the variation of look angles across the satellite swath. The uniform fault strike assumption of the Xianshuihe fault is only used in constructing the fault-parallel and the fault-perpendicular velocities. This assumption is not applied to the fault slip modelling here. We used the ascending/descending LOS rates directly and we did not use the fault-parallel rates in modelling interseismic motion. The interseismic fault slip model is set up in an elastic half-space. We modelled the interseismic deformation with free slip on planar fault extending from the locking depth to 1000 km deep. The locking depth is assumed to be 12 km. The fault plane is assumed to be vertically dipping. We divided the fault model into five segments based on the surface trace of the Xianshuihe fault. The northern and southern ends of the fault model has been extended for about 100 km to avoid the edge effect of the elastic dislocation model (Okada 1985). The fault model is allowed to slip only in strike-slip direction so the dip-slip component is set to be zero. We impose a positivity constraint to allow left-lateral strike-slip on the Xianshuihe fault. After setting up the interseismic locking model we invert for the geodetic fault slip rate of the Xianshuihe fault. The slip rate is mainly constrained by the InSAR LOS velocity data. We attempted to vary the locking depth to investigate the effect of varying locking depth on the slip rate estimates and the result is described in the following section.

We developed the fault creep model by inverting the residuals of the interseismic locking model. The residuals are derived by subtracting the prediction of the interseismic locking model from the observed InSAR LOS velocity data. The fault creep model extends from the surface to a depth of 12 km in an elastic half-space. The strike of the fault creep model follows the surface trace of the Xianshuihe fault. The fault plane of the fault creep model is discretized into small rectangular patches. The size of the discretized fault patch is 1 km by 1 km at the top layer and it increases with depth to accommodate lower resolution at deeper part of the model. A smoothing factor is applied to the fault creep model to find the smoothest model that fit the data reasonably well. The smoothing constraints penalize the first derivative of the creep model along both the along-strike and along-dip directions. We have experimented with different values of the smoothing factor (0.2, 0.5 and 1). To view more details of the creep model we chose the preferred model with a smoothing factor of 0.5. In the geodetic inversion problem, we employed the non-negativity constraint to allow for only left-lateral strike-slip. We have added another constraint to minimize slip at the bottom of the creep model to produce more realistic models, because the slip at much deeper depth cannot be constrained tightly by the InSAR data due to large atmospheric noise far away from the Xianshuihe fault.

4 RESULTS

4.1 Secular fault creep model

We used GMTSAR to process the Sentinel-1 TOPS data to derive the LOS velocity maps as shown in Fig. 3. These maps are obtained from stacking long-term interferograms spanning from the end of 2014 to the beginning of 2021. A remove-restore method based on GPS velocities (Wang & Shen 2020) is used to correct for long-wavelength errors of the InSAR data (Tong et al. 2013). The long-wavelength errors could be caused by inaccurate orbit or large-scale atmospheric and ionospheric sources. We assumed that the six-year InSAR stacking result can represent the long-term secular movement of the crust. The mean LOS velocity results show an increase by a few mm yr⁻¹ for the ascending orbit when it goes from the northeast to the southwest side of the Xianshuihe fault and it shows a decrease by a few mm yr⁻¹ for the descending orbit. This sense of motion agrees with the tectonic movement of a left-lateral strike-slip fault. In Fig. 3, the shading of the velocity map indicates the gradient of the LOS velocity field. The shading highlights the places where LOS velocity changes abruptly. It is found that along the Xianshuihe fault there is a sharp change in the LOS velocity, indicating fault creep is occurring along the fault. From the LOS velocity gradient, we identified a long creeping segment along the Xianshuihe fault and the length of the surface creeping reaches 200 km approximately.

We combined the LOS velocity and solved for the fault-parallel velocity and vertical velocity (Fig. 3). To the northeast of the Xianshuihe fault, the fault-parallel velocity reaches 5 mm yr⁻¹ and to the southwest of the Xianshuihe fault it reduces to −5 mm yr⁻¹. As shown by the GPS velocity data (Fig. 3), the main motion of the crust is the left-lateral shear across the Xianshuihe fault (King et al. 1997; Chen et al. 2000; Shen et al. 2005; Zheng et al. 2017). The fault-parallel velocity from InSAR also shows a clear difference across the Xianshuihe fault. Since the major movement is horizontal shear in this region the vertical velocity does not change significantly across the Xianshuihe fault.

Next we took three profiles of the fault-parallel and vertical velocities which are perpendicular to the Xianshuihe fault. The profiles are approximately 200 km long and the coverage of these profiles can be seen in Fig. 3. In Fig. 4, the fault-parallel velocity on the northeast side of the profiles is higher than the southwest side by a few mm yr⁻¹. The horizontal velocity from the GPS data is projected into the LOS direction of the satellite and are plotted together with the InSAR results. The GPS results match well with the LOS velocity profiles and there is a general agreement between GPS and InSAR. This comparison has validated the InSAR results and it shows that the methods we took in data processing are reasonable. Fig. 4 shows the three profiles of vertical velocity and we found that the deformation rate in the vertical direction is subtle, only 1–2 mm yr⁻¹. The northeast side of the Xianshuihe fault is moving relatively upward by 1–2 mm yr⁻¹ compared to the southwest side of the fault.

A clear signal of surface creeping is observed along the Xianshuihe fault in both the LOS velocity and the fault-parallel velocity. As shown in Fig. 5, the gradient of the fault-parallel velocity delineate the surface trace of the Xianshuihe fault. We found that surface creep is occurring to the north of Kangding for a length of about 200 km along the Xianshuihe fault. Fig. 5(b) shows the creep rate as a function of distance along the Xianshuihe fault. The optimal choice of the cross-section length is determined by trial and error. We have experimented with the 1 and 5 km across-fault distance and the result is shown in Fig. 5(b). Because we used 2 km low-pass filter to the LOS velocity map, the minimum resolution of the velocity map is 2 km. The 1 km across-fault distance is lower than the resolution of the velocity map so the estimated creep rate is close to zero. The creep rate estimation from the 5 km across-fault distance is significantly smaller than previous studies (Li & Bürgmann 2021). Thus, we chose the across-fault distance to be 10 km to obtain the optimal estimation of creep rate. We also used the 10 km across-fault distance to determine the temporal variation of the surface creep rate. The fastest creeping region is to the south of Bamei and the surface creep rate is about 6 mm yr⁻¹. The surface creep rate decreases to 4 mm yr⁻¹ near Daofu region and to 2 mm yr⁻¹ near Luhuo region. It decreases to nearly 0 mm yr⁻¹ near Kangding. Apparently, the surface creep rate variation exhibits a wave-like pattern: there are three peaks and three troughs in the creep rate profile.

The interseismic fault slip model is constrained by the surface velocity from InSAR. The surface velocity data have a spatial sampling of hundreds of metres so there is a redundancy in the surface deformation provided by InSAR. We sampled the InSAR LOS velocities using two different schemes. One is for the fault slip rate estimation and the other is for the fault creep rate estimation. To avoid the creep signal affecting the estimation of slip rate we have masked the near-fault data out. Fig. 6 shows the InSAR LOS velocity data points, the predicted LOS velocities from the best-fitting model, and their residuals. For the fault creep rate estimation, we varied the downsample spacing based on the distance to the fault, with denser points near the fault and sparser points in the far field. This approach helps to eliminate the effect of atmospheric noise on the creep rate inversion (Fig. 7). As seen in Fig. 6, the modelled velocities capture the first-order variations of the surface velocity. We found that the best-fitting fault slip rate is 11.5 mm yr⁻¹ by modelling the InSAR LOS velocities. The modelling also indicates that a uniform slip rate of 11.5 mm yr⁻¹ is essential to explain the surface velocity observations. Thus, a variation of the fault slip rate along the Xianshuihe fault is not required by our modelling (Zhang et al. 2022). The geodetic slip rate estimate depends on the assumption on the locking depth. To explore the influence of this assumption, we have tried a different value of the locking depth, that is, 20 km. We found that the fault slip rate increased to 14 mm yr⁻¹ when modelling the interseismic velocity with a 20 km locking depth. The residuals of the interseismic fault slip model exhibit a systematic pattern near the surface trace of the Xianshuihe fault (Fig. 7). This pattern is related to the fault creep of the Xianshuihe fault. The magnitude of the residuals is on the order of 3 mm yr⁻¹.

Fig. 7 shows the data, model prediction and the residual from the fault creep model. The input of the fault creep model is the residual of the interseismic fault slip model. Fig. 8 shows the slip distribution of the secular fault creep models for different smoothing factors. The creep model with a smoothing factor of 1 is overly smoothed and we cannot detect any details of the creep distribution. The prominent feature of this smooth model is the fast-creeping area beneath Bamei. The model with a smoothing factor of 0.2 shows that the creeping patches are more isolated and scattered. This model has a creep rate much larger than the long-term fault slip rate, leading to a much negative interseismic coupling ratio. The preferred model has a smoothing factor of 0.5 and it is between the smooth and the rough models. Although the temporal window is only 6.5 yr and the fault creep might vary through time we deem this fault creep model a secular model. In the preferred model, the maximum creep rate reaches 20 mm yr⁻¹ and it is distributed from south of Bamei to north of Kangding. This rate is slightly larger than the long-term fault slip rate so the interseismic coupling ratio of those sections can be interpreted to be less than zero. This moderate negative interseismic coupling ratio is consistent with previous study by Li & Bürgmann (2021). The other fault segments from Luhuo to Bamei have smaller creep rates around 3–5 mm yr⁻¹. There are two sections of the Xianshuihe fault exhibit close-to-zero creep rates, that is, high interseismic coupling: one is from Zhuwo to Luhuo and the other is from Moxi to Shimian. These deep creeping signals should be distinguished from the surface creeping signal obtained directly from the InSAR velocity map.

4.2 Time-dependent fault creep model

The multitemporal stacking method yields LOS velocities for three time periods: 2015–2017, 2017–2019 and 2019–2021 (Fig. 9). To illustrate the fault creep signal, we simply removed a linear trend in the LOS velocity grid after stacking. We did not carry out the remove-restore method to correct for the long-wavelength errors. We applied the remove-restore method to the InSAR LOS velocities later when modelling the interseismic locking and fault creep, similar to what we have done for the secular fault creep model. We can identify that the Xianshuihe fault is creeping from both the ascending and descending orbits. Because of the looking geometry of the satellite and the fault strike orientation the surface creeping signal is more obvious in the descending orbit than it is in the ascending orbit. In Fig. 9, the near-fault deformation is clearest during the 2015–2017 period. The 2017–2019 period show a large near-fault deformation pattern to the southwest of the Xianshuihe fault. The deformation signal from 2017–2019 may be contaminated by the atmospheric noise. We attempted to adjust the input interferograms in the multitemporal stacking and the result did not change significantly. The 2019–2021 period shows the least amount of near-fault deformation. In Fig. 9, the gradient of the LOS velocity delineates the fault creep signal clearly for all three time periods. During the 2015–2017 period, there is a large blob of positive signal far away from the Xianshuihe fault that can be attributed to atmospheric noise. In the following inversion step, we downsampled the InSAR LOS data using a non-uniform scheme to mitigate the effect from the atmospheric noise.

After combining the LOS velocity from both the ascending and descending orbits we obtained the time-dependent fault-parallel velocity (Fig. 9). A clear linear feature along the Xianshuihe fault can be identified for all three time periods. The 2015–2017 and 2017–2019 periods show significant near-fault deformation that can be related to fault creep. The 2019–2021 period shows smallest near-fault deformation comparatively.

Because the time-dependent fault creep is a subtle signal to study and the atmospheric noise is significantly large, we compared the results from multitemporal stacking with coherence-based SBAS (Fig. 10). To first order we found that the LOS velocities derived by multitemporal stacking looks similar to the results derived by SBAS. For example, the LOS velocities along the descending orbit from 2019 to 2021 (Figs 9f and 10f) looks almost the same. In some cases, that is, descending orbit from 2015 to 2017 (Figs 9d and 10d), the fault creep signal from the multitemporal stacking is cleaner than what is derived by SBAS. To summarize the time-dependent deformation from multitemporal stacking is better than the results derived from SBAS. Thus, we used the LOS velocities from multitemporal stacking in the fault creep inversion.

We estimated the time-dependent surface creep rate along the Xianshuihe fault by fitting linear functions to the cross-sections of the fault-parallel velocities. The result is illustrated in Fig. 5(c). The spatial-temporal variation of the surface creep is complicated and there is no clear pattern we can follow. From 2015 to 2021, the surface creep rate changes spatially and temporally within the range from 0 to 6 mm yr⁻¹. In general, the surface creep rate is largest for the 2017–2019 period. At distance 40 km along the Xianshuihe fault (near Kangding), there is a temporal decay of surface creep rate where it went from 4 mm yr⁻¹ at 2015–2017 to 0 mm yr⁻¹ at 2019–2021. There is an earthquake doublet in Kangding in 2014 so this temporal decay might be related to the afterslip triggered by the Kangding earthquakes, in agreement with previous study (Li & Bürgmann 2021; Huang et al. 2023). At other locations, we did not find a systematic temporal decay.

We developed the time-dependent fault creep model by inverting the near-fault deformation imaged by multitemporal stacking. Figs 11–13 show the data, model and residuals from the time-dependent fault creep modelling. The across-fault profiles of the data/model comparison are shown in Fig. 14. Fig. 15 shows the distribution of deep creep rate of the Xianshuihe fault for the three time periods: 2015–2017, 2017–2019 and 2019–2021. It can be seen that the fault creep is not steady but varies drastically with time. During 2015–2017 period, the fault is creeping beneath Daofu at a rate of 5 mm yr⁻¹. From Bamei to north of Kangding most of the fault patch is creeping at 10–20 mm yr⁻¹, slightly larger than the geodetic fault slip rate. To the south of Bamei, there is a small fault area that is creeping at 30–40 mm yr⁻¹. During 2017–2019 period, the fault is creeping for a long distance (from Zhuwo to Bamei) at a rate of 5 mm yr⁻¹. From Bamei to north of Kangding the fault creep rate increases to 40 mm yr⁻¹, exceeding the fault slip rate by almost a factor of four. During the 2019–2021 period, the fault creep diminishes to 5 mm yr⁻¹ or less from Bamei to Kangding. From Luhuo to Daofu, the fault creep rate is low and it seems that the fault become partially locked. There is also significant temporal variability in fault creep rate near the northern end (Zhuwo) of the Xianshuihe fault.

5 DISCUSSIONS

Past studies from short-baseline survey and InSAR have revealed the along-strike variation of surface creep rate. We compare the averaged surface creep rate estimated by this study with previous studies. These studies all used the Sentinel-1 InSAR data with similar time interval (i.e. from 2014 to 2021) and employed similar methods to estimate surface creep rate. Here, we focus on the secular component of the creeping section. Qiao & Zhou (2021) found that the surface creep rate along the entire Xianshuihe fault varies between 0 and 6 mm yr⁻¹. Li & Bürgmann (2021) estimated the surface creep rates are 2.5 mm yr⁻¹ near the Daofu rupture and 5 mm yr⁻¹ near the Kangding earthquake doublet. In our study, the creep rate at Daofu is ∼ 2 mm yr⁻¹, and the creep rate near Kangding earthquakes is ∼ 5 mm yr⁻¹. Our study is in good agreement with Qiao & Zhou (2021) and Li & Bürgmann (2021). Their along-strike variation of the surface creep rate also shows a wave-like pattern, similar to our estimation (Fig. 5). Huang et al. (2023) identify a creeping section along the Kangding segment with a creep rate distributed from 0 to 12 mm yr⁻¹, a factor of two larger than our estimation.

Li et al. (2021) used GPS and InSAR data to develop an interseismic fault coupling model along the Xianshuihe–Zemuhe–Anninghe–Xiaojiang fault system. We found that the Xianshuihe fault is creeping at shallow depth from Zhuwo to Shimian. We also found that the Xianshuihe fault has a fast-creeping segment located between Bamei and Kangding. This finding is supported by the decoupling results from Li et al. (2021). Li & Bürgmann (2021) inverted an interseismic shallow slip model with Sentinel-1 InSAR observations, suggesting that there is a large creeping patch beneath Bamei and the creep rate reaches 25 mm yr⁻¹, larger than the tectonic loading rate. This result is similar to the aseismic slip model from Huang et al. (2023). They interpret this creeping patch to be caused by postseismic afterslip from 2014 Kangding earthquake doublet. Our secular fault creep model is in general agreement with the model from Li & Bürgmann (2021). In our model, the creeping rate reaches 20 mm yr⁻¹ slightly larger than the long-term fault slip rate (Fig. 9). Thus, there is limited elastic energy accumulating on this creeping segment.

The 2022 M_w 6.6 Luding earthquake occurred at the southern end of the Xianshuihe fault, between the town Moxi and Shimian. By comparing the location and distribution of the coseismic slip of the Luding earthquake and the features of the secular fault creep model we found that the Luding earthquake occurred at the partially locked portion of the Xianshuihe fault. Our finding is consistent with published studies (Guo et al. 2023): the Luding earthquake only broke a part of the locked Moxi segment of the Xianshuihe fault. The other part of the Moxi segment is still accumulating seismic moment, which will be released in the form of big earthquakes in the future.

In our study, there is a wide spread signal of fault creep along the entire Xianshuihe fault. The fault creep not only occurred at surface but also exists at depth in the upper part of the crust. As pointed out by previous studies, fault creep will slowly release elastic energy on the fault plane. Thus, the finding of secular fault creep is an indicator of reduced seismic hazard. To quantify this, we calculated the moment accumulation rate and the moment release rate. If we assume that the shear modulus of the crust is 30 Gpa, we found that the moment release rate from fault creep is 4.9 × 10¹⁷ N·m yr⁻¹. If we assume that the total length of the Xianshuihe fault is 360 km long and the locking depth is 12 km and we found that the slip rate of the Xianshuihe fault is 11.5 mm yr⁻¹, the moment accumulation rate of the Xianshuihe fault will be 1.5 × 10¹⁸ N·m yr⁻¹. After accounting for fault creep, the moment accumulation rate is reduced to 1.0 × 10¹⁸ N·m yr⁻¹. This moment accumulation rate is equivalent to an M_w 7.3 earthquake for every 100 yr if there is a rupture through the entire Xianshuihe fault.

The fault creep of the Xianshuihe fault is ultimately driven by long-term tectonic loading. However, the behaviour of fault creep is hardly steady. Based on evidence from trilateration, GPS and InSAR we can draw a conclusion that the aseismic slip distribution on the Xianshuihe fault changes dynamically with time (Li et al. 2021; Qiao & Zhou 2021; Huang et al. 2023). The aseismic slip might accelerate immediately after major earthquakes and then slowly decay with time. It might also migrate along the fault due to mechanisms such as: earthquake swarms and fluid migration. We evaluate our time-dependent fault creep model (Fig. 15) by comparing our model with the study from Huang et al. (2023). Huang's model shows that the aseismic slip rate of the Kangding segment decays from 40 mm yr⁻¹ in 2015 to 0 mm yr⁻¹ in 2021. Our model, however, shows a rather complex temporal evolution of the aseismic creep along the segment from Bamei to Kangding. The maximum creep rate first increases from 30 mm yr⁻¹ in 2015–2017 to 40 mm yr⁻¹ in 2017–2019 and then decreases to 5 mm yr⁻¹ in 2019–2021. This feature is consistent with the surface deformation pattern derived by multitemporal stacking (Fig. 9) and the temporal variation of the surface creep rate (Fig. 5c). Both our model and Huang et al. (2023) model confirm a relocking of the Kangding segment from 2019 to 2021, that is, the creep rate reduced from 2014 to 2021, but we did not find a simple temporal decay of fault creep rate after 2014, when the Kangding earthquake doublet occurred.

The sudden acceleration of aseismic slip rate from 2017 to 2019 demands an explanation (Fig. 14). Because there are no large earthquakes occurred from 2017 to 2019 in this region (Fig. 1) we can rule out the possibility of static stress triggering or afterslip. There are two possible ways to explain our model. First, aseismic slip could be triggered dynamically by seismic waves excited at far distance. Tymofyeyeva et al. (2019) found a slow slip event along the southern San Andreas fault triggered dynamically by a distant earthquake. The magnitude of the transient creep reached several millimetres and the duration of the creep event is several months. Similarly, the Xianshuihe fault could experience a dynamic stress perturbation from passing seismic waves and start to slip unsteadily. Unfortunately, we cannot identify a particular event that triggered the slow slip due to low temporal resolution of the InSAR data and a lack of creepmeter data. The second explanation is that the acceleration of creep rate could be caused by fluid migration. Geophysical evidence suggest that the fluid can migrate into crust from upper mantle leading to a weakened and creeping crustal fault (Becken et al. 2011). Overpressure from fluid injection is able to assist aseismic slip in hydraulic fractured region, which may explain the acceleration of aseismic creep observed in this study (Eyre et al. 2019).

6 CONCLUSIONS

In this study, we processed Sentinel-1 InSAR data from 2014 to 2021 by stacking long-term interferograms to derive surface velocity maps along the Xianshuihe fault in China. Both the ascending and descending LOS velocity are combined to obtain the fault-parallel and vertical velocity field. We found that the Xianshuihe fault show a clear surface creep signal and the surface creep rate reaches as high as 6 mm yr⁻¹. The length of the fault showing surface creep reaches ∼200 km. A large fraction of the interseismic velocity map can be explained by deep dislocation with fault slip rate of 11.5 mm yr⁻¹. Residuals of the deep dislocation from the interseismic deformation can be fit with shallow fault creep in the upper 12 km of the crust. One significant creeping patch is located between Bamei and north of Kangding.

We applied the multitemporal stacking method and the coherence-based SBAS method to the Sentinel-1 InSAR observation to study the time dependence of surface deformation around the Xianshuihe fault. The six-year time span of the InSAR data is divided into three periods: 2015–2017, 2017–2019 and 2019–2021. The temporal variation of the surface creep rate shows a complex pattern. The time-dependent fault creep modelling indicates that the fault creep rate from Bamei to Kangding accelerated to a maximum of 40 mm yr⁻¹ in 2017–2019 and decayed to 5 mm yr⁻¹ in 2019–2021. The aseismic slip on the Xianshuihe fault can vary significantly in a short time period (a couple of years) and over a short distance (tens of kilometres). From 2015 to 2021, the fully creeping segment of the Xianshuihe fault (from Bamei to Kangding) seems to transform to a partially locked segment. This study is important to earthquake hazard assessment because it implies that time dependence should be carefully examined when using geodetic data to estimate interseismic coupling or moment deficit along active faults. To better understand fault mechanism and better forecast earthquakes, we need to investigate more details of the time-dependence of fault creep for active faults in the future.

ACKNOWLEDGMENTS

The figures in this study were generated using the public domain Generic Mapping Tools software (Wessel et al. 2013). This research is funded by the National Key Research and Development Program of China (2023YFC3007305 and 2023YFE0101800) and the Special Fund of the Institute of Geophysics, China Earthquake Administration (grant nos DQJB23B24 and DQJB22X12). We thank Dr Hua Wang and an anonymous reviewer for their critics and suggestions on this manuscript.

DATA AVAILABILITY

Sentinel-1 data were provided by ESA and accessed from ASF.

REFERENCES