Illuminating the Spatio-Temporal Evolution of the 2008–2009 Qaidam Earthquake Sequence with the Joint Use of Insar Time Series and Teleseismic Data

: Inferring the geometry and evolution of an earthquake sequence is crucial to understand how fault systems are segmented and interact. However, structural geological models are often poorly constrained in remote areas and fault inference is an ill-posed problem with a reliability that depends on many factors. Here, we investigate the geometry of the Mw 6.3 2008 and 2009 Qaidam earthquakes, in northeast Tibet, by combining InSAR time series and teleseismic data. We conduct a multi-array back-projection analysis from broadband teleseismic data and process three overlapping Envisat tracks covering the two earthquakes to extract the spatio-temporal evolution of seismic ruptures. We then integrate both geodetic and seismological data into a self-consistent kinematic model of the earthquake sequence. Our results constrain the depth and along-strike segmentation of the thrust-faulting sequence. The 2008 earthquake ruptured a ∼ 32 ◦ north-dipping fault that roots under the Olongbulak pop-up structure at ∼ 12 km depth and fault slip evolved post-seismically in a downdip direction. The 2009 earthquake ruptured three south-dipping high-angle thrusts and propagated from ∼ 9 km depth to the surface and bilaterally along the south-dipping segmented 55–75 ◦ high-angle faults of the Olonbulak pop-up structure that displace basin deformed sedimentary sequences above Paleozoic bedrock. Our analysis reveals that the inclusion of the post-seismic afterslip into modelling is beneﬁcial in the determination of fault geometry, while teleseismic back-projection appears to be a robust tool for identifying rupture segmentation for moderate-sized earthquakes. These ﬁndings support the hypothesis that the Qilian Shan is expanding southward along a low-angle décollement that partitions the oblique convergence along multiple ﬂower and pop-up structures.


Introduction
Inferring earthquake slip distributions is useful for understanding the spatial and temporal evolution of crustal strain, in order to investigate how faults are segmented and interact, and to study the relationship of an earthquake with the structural tectonics of an area. However, earthquake fault inversion is an ill-posed problem. Its solution is limited by data resolution, model simplifications and biases, the non-linearity between the data and the unknown parameters, the impact of the topography and lateral variations of lithology, and/or significant dependence between inverted parameters e.g., [1][2][3][4]. The reliability of any earthquake fault inversion depends on many factors, including the segmentation and size of the rupture, the quality of the data, the uncertainties in the subsurface model, the chosen optimisation method to infer earthquake fault parameters, and the a priori inferences of the earthquake geometry or mechanisms. For these reasons, different research groups often find different solutions for the same earthquake [5].
One way to reduce the non-uniqueness of an inversion and the impact of data errors is to combine several data sets with different characteristics (e.g., space-based imaging, seismic data, pixel offsets, tsunami records) e.g., [6][7][8][9]. In particular, several studies e.g., [4,[10][11][12][13] have shown the benefit of complementing the high spatial resolution of near-field geodetic observations with the high temporal resolution of far-field seismic waveforms. Seismic waveforms carry information about the rupture onset time, its spatio-temporal evolution and duration, as well as the fault mechanism and seismic moment. Near-field surface data provide constraints on the fault location and mechanism, the amount of slip, and the extension of the rupture. Another way to limit biases in the inversion is to account for data and model parameter uncertainties due to the data noise and the imperfect knowledge of the medium through data weighting e.g., [1,[14][15][16], and to rigorously propagate errors through optimisation algorithms, such as with un-regularised Bayesian optimisation frameworks e.g., [2,9,13,[17][18][19][20] Differential Interferometric Synthetic Aperture Radar (DInSAR) detects surface displacements in the line-of-sight (LOS) direction between satellite acquisitions before and after an earthquake, and it is commonly used to infer earthquake fault parameters due to its high spatial resolution e.g., [21][22][23][24][25][26]. Despite the large benefits that are brought by the measurement of the surface displacements across wide swaths, DInSAR is also polluted by changes of atmospheric conditions between the two acquisitions, associated with the turbulent and stratified components of the atmosphere and the electron content in the ionosphere [27][28][29]. In addition, due to the revisit time of the satellites of more than a few days, most inferred earthquake fault models based solely on DInSAR observations are biased by early post-seismic displacement signals. Apart from few studies that account for early post-seismic signals e.g., [30,31], most of the fault inferences neglect the contribution of early post-seismic deformation in the inverse problems with the consequences that it is attributed to the co-seismic rupture. However, since the first surface displacement image of an earthquake following the 1992 Landers earthquake [32], numerous advances in the DInSAR technique have been made to derive multi-temporal time series (TS) of ground displacements by processing multiple SAR acquisitions. This increases the number of locations where displacement can be measured and reduces associated errors e.g., [33][34][35]. Advances have also been made in reducing atmospheric contributions by the use of atmospheric models e.g., [28,36,37], empirical methods e.g., [35,[38][39][40][41][42] or by data stacking and TS analysis e.g., [30,[43][44][45]. The latter is beneficial to derive the cumulative surface displacements (contrary to single interferograms) and separates the various inter-, co-and post-seismic phases of the seismic cycle. In addition, TS analysis allows researchers to check the consistency of each interferogram within the network, increases the signal-to-noise ratio via data redundancy and attenuates the effects of turbulent atmospheric noise by smoothing of the TS. Recent developments in InSAR methodologies as well as the new generation of SAR satellites (Sentinel-1, ALOS-2) with shorter revisit times and systematic acquisitions in both ascending and descending geometries, therefore, raise queries about the prospects of multi-temporal InSAR TS for subsurface fault slip inferences.
The North Qaidam thrust (NQT) system is located at the southernmost part of the Qilian Shan-Nan Shan and is marked by recent active seismicity within the Olongbulak ranges ( Figure 1). The ∼16 km deep 10 November 2008 and the ∼5 km deep 28 August 2009 Haixi M w 6.3 earthquakes occurred in close proximity (∼15 km) to each other within the Olongbulak Shan and were followed by a period of increased rate of seismicity in the region [46,47]. The 2009 earthquake was followed by a M w 5.4 and a M w 5.6 aftershock one and two days, respectively, after the main-shock, as many as 34 M w > 4 aftershocks within a month (USGS), and time-dependent after-slip of similar pattern to the co-seismic slip [48][49][50][51]. The two earthquakes represent particular examples of two ∼M w 6 earthquakes that occurred nearby in space and time, highlighting the need to constrain both along-depth and along-strike segmentation of thrusting events to better quantify seismic hazards [46]. Their segmentations and geometries also need to be placed in the context of the regional tectonic as the deformation style of the northeastern part of the Tibetan plateau is at the heart of the debates about its dynamic evolution. In a first class of deformation models, Burchfiel et al., 1989 [52] first postulated a basement décollement fault in the middle crust of the Qaidam basin that connect to the Qilian Shan-Nan Shan thrust belts in the north, absorbing progressively the left-lateral slip on the Altyn-Tagh Fault ( Figure 1). Subsequently, Metivier et al., 1998 [53], Meyer et al., 1998 [54] and Tapponnier et al., 2001 [55] suggested that the Qaidam basin may have been trapped of the surrounding reliefs in the middle-late Miocene as a result of the northward stepwise jumping of the intra-continental deformation. This first class of model proposes that the oblique convergence between the India-Eurasia collision is absorbed by crustal thickening along the major Kunlun or Qilian Shan suture zones, where strike-slip and thrust faults root at depth to extrude laterally and partition vertically the deformation causing mountain uplift on localised ranges ( Figure S1a) [58] suggested a plausible and simultaneous early Eocene onset of the Kunlun and Qilian ranges. This second class of model suggests that the strain from the India-Asia collision was rapidly transmitted to the northern edge of the modern Tibetan Plateau via preexisting and complex zones of weaknesses created during pre-Cenozoic [57,59,60] or via homogeneous deformation driven by mantle or crustal flow [61][62][63] ( Figure S1b). Those model contradictions leave thus open the debate on the geometry of the NQT at depth and its possible relations with the Qilian ranges, to the north [55,59]. The study of these earthquakes is also challenging from a methodological point of view. First, because of the deeper source of the 2008 earthquake, the data can be explained by conjugate south-dipping or north-dipping planes, as illustrated by the various solutions found by different authors [46][47][48]. Secondly, the larger gradients and more complex surface displacement patterns of the 2009 earthquake, in comparison to the 2008 earthquake, suggest a segmented and shallow rupture. These surface displacement features have been explained by varying numbers of fault segments in the literature [46,48,65], and this earthquake sequence, therefore, also represents a good case study to evaluate the level of complexity that is required in a fault model to explain the data.
In this study, we investigate the spatio-temporal evolution of the 2008 and 2009 earthquakes using seismological back-projection analysis and three overlapping InSAR time series (TS) data of the Envisat satellite covering the two earthquakes. We then integrate both geodetic and seismological data into a self-consistent fault model. The goals of this study are multi-folds. We first investigate the fitness of back-projection to describe the rupture characteristics of moderate-size earthquakes. We then inspect the benefits of an InSAR TS analysis approach for earthquake fault modelling and discuss the improved signal to noise ratio of such data sets in comparison to traditional co-seismic interferograms, as well as the additional constraints that are offered by the extracted post-seismic surface displacements from InSAR time-series data. We finally integrate our fault characteristic estimations to the local tectonic context to propose a structural model of the Olongbulak ranges.

Seismic Back-Projection Source Imaging
Time-domain seismic back-projection images the spatio-temporal evolution of earthquakes by mapping the location of coherent seismic wave emissions recorded at a seismic network e.g., [66][67][68][69]. It relies on constructive and destructive stacking of seismic waveforms. In comparison to kinematic fault inversion, it has the advantage of requiring no assumption of the fault geometry, and is, therefore, a useful complement to inverse modelling. The only model assumptions required are the travel times derived from a one-dimensional (1-D) Earth velocity model.

Data and Station Clustering
We first download all available broadband records from the international web service (FDSNWS) Geofon and IRIS [70,71], for the 2008 and 2009 Qaidam earthquakes. We consider waveforms from broadband stations between 20 • and 93 • distance from the source and only use stations that recorded both earthquakes ( Figure S2). We cluster the available stations using the k-means algorithm into many small virtual arrays (minimum number of five stations). Virtual arrays are formed assuming a minimum number of stations and a maximum aperture up to 3.5 • . The use of many small virtual arrays has the advantage of minimising the effect of velocity differences between stations in an array as well as the effect of radiation patterns and source directivity [72].

Muti-Array Back-Projection Method
We carry out the back-projection for point locations on a horizontal grid. We perform the back-projection separately for high-frequencies (HF) and low-frequencies (LF), in order to explore areas with changes of frequency content over time, which can indicate changes in slip velocity/acceleration [73] or rupture initiation points [74]. We perform the back-projection on moving time windows of 2 s steps and 10 s length for HF and of 8 s steps and 26 s length for LF. We bandpass-filter for the HF back-projection between the earthquake corner frequency (2/T r , with T r being the reported gCMT rupture duration [64]) and 1.5 Hz. For the LF back-projection we bandpass-filter between 0.01 Hz and the earthquakes' estimated corner frequency f c = 0.25 Hz.
We use a phase-weighted stacking method [75] and stack in separate time windows for P-and SH-phases. Waveforms are back-projected for each time step and each point on a given spatial grid.
Un-modeled path effects cause systematic travel-time shifts, which lead to a regularly biased source location. In order to reduce the impact of this effect, we calibrate travel-time shifts for each station using empirical corrections [76][77][78][79] based on nearby 10 November 2008 03:47:23 Mw 5 aftershock to calibrate the empirical travel-time and correct the travel time shifts of the P and SH waves of the two earthquakes ( Figure S3). As a result of the back-projection, we obtain the phase semblance (a measure of phase coherence) of the waveform records at each array for the specific time window and grid point considered. Semblance is a strictly positive function of waveform coherence and can be interpreted as a likelihood of seismic energy released [72]. We also calculate the beampower, which is an absolute measure of the stacked amplitudes [80]. Beampower is calculated by integrating the energy of all arrays.
We make use of bootstrapping and travel-time perturbations derived from the 1-D velocity model to access the uncertainty of the obtained back-projection results. For 100 bootstraps, we calculate the back-projection with different array weights and perturb the waveform records randomly by ±2 s for the P-phase window and by ±4 s for the SH-phase, extracted from the standard deviation of the normal distribution. To avoid azimuthal biases, we subdivide the azimuth into 12 sectors of 30 • , based on the azimuth of the earthquake epicenter to the array center, and each array is assigned to its corresponding azimuth sector. The semblance from each azimuth sector is normalised, such that each azimuthal sector has the same influence on the combined semblance.

Kinematic Fault Inversion
We seek to explore the range of possible fault geometries for the two M w > 6, 2008 and 2009 earthquakes that explain the surface displacements and the seismological observations with minimum prior constraints on fault geometry and fully accounting for data uncertainties. We optimise both near-field and far-field data with a Bayesian bootstrapping algorithm implemented in the open source Grond earthquake inversion framework [4,81,82] and part of the Pyrocko software (https://pyrocko.org/grond/docs/current/method/index.html).

Teleseismic Data
We select stations at epicentral distances between 20 • and a maximum of 93 • . Those distances are large enough, such that near-field wave terms are attenuated and small enough such that phase triplications do not occur and the direct P and S waves are recorded first. We visually check and select the data quality after the automated pre-selection. As a result, we retrieve for the 2008 and 2009 earthquakes good-quality recordings at 72 and 45 stations, respectively. The waveforms are trimmed and filtered between the corner frequency ( f c = 0.25 Hz) and 1.5 Hz, and they are then cut and filtered with a frequency taper falling with a cosine flank.
For the 2009 earthquake, we observe several waveform records contaminated by a signal from a M w 6.9 earthquake that occurred in the Banda Sea ∼1 min prior. Thus, we select stations within an azimuth of −105 • and 105 • and between distances of 0 • to 40 • , where the 2009 earthquake signals arrived before the Banda Sea earthquake signals. Additionally, we only use P-phases, because the SH-phases are potentially contaminated.

Insar Time Series (Ts) Data
We construct InSAR TS of surface displacements from the complete Envisat archive along two ∼300 km long overlapping descending tracks (T319, T047) and one ascending track (T455) ( Figure S4). Extensive details about the filtering and unwrapping procedure can be found in Daout et al., 2019 [45]. In contrast to the processing of only few co-seismic interferograms, here we take advantage of the regional extent of the data and the large number of interferograms to perform empirical corrections after unwrapping, which we validate with a time-series analysis approach. First, we correct the 215 interferograms for orbital and atmospheric errors. To do so, within non-deforming regions, we jointly estimate the LOS-elevation relationship with the azimuth and the LOS-range relationship, as follows: LOS(x, y) = az + bzy where LOS(x, y) is the LOS displacements for each interferogram, x, y, z are, respectively, the range direction, the azimuth direction, and the elevation, and a, b, c, d, e, and f the inverted coefficients. Inconsistencies of parameters a, b, c, d, e, and f are afterwards detected within the interferometric network by a least-squares inversion of all independent coefficients for each acquisition date and then re-estimated for all interferograms before applying the correction. After empirical corrections, the phase TS is constructed with the method of Doin et al., 2015 [35] and inconsistencies in the interferometric network, mainly due to phase unwrapping errors, are further detected and corrected with iterative and automatic approaches [83]. We finally extract the co-seismic and post-seismic surface displacements of the 10 November 2008 and the 28 August 2009 Haixi earthquake by decomposing pixel-by-pixel the TS data into a linear trend, two steps, and two logarithmic post-seismic functions (Figure 2a), such as: ). ( where t is the SAR acquisition time, v, the long-term velocity between 2003 to 2011, A, B, and τ, the amplitude of the Heaviside step functions (H), the amplitude of logarithmic functions, and the characteristic time of the logarithmic relaxation, respectively ( Figure 2a). Given the moderate size of the shallow earthquakes, the short-time period observation following them, the short-wavelength of the post-seismic signal, and the similar signs and patterns of the post-seismic and co-seismic surface displacements, we discard any viscoelastic relaxation or poro-elastic rebound in the shallow crust e.g., [84][85][86] and associate the deformation to afterslip [87,88]. We explore several post-seismic relaxation times, τ 2008 and τ 2009 , but notice no significant differences in the misfit due to the temporal sampling of few months of the data and the moderate magnitude of the earthquakes and associated after-slip ( Figure S5). However, the choice of the relaxation time does not influence the interpretation of the results, and choose a characteristic time of half a month that minimises the misfit of the decomposition, as we only interpret and invert later the modeled cumulative after-slip displacements. We impose that displacements from the logarithmic functions must have the same sense of direction as their associated coseismic steps, and solve the constrained inverse problem with a sequential least-squares algorithm. The resulting co-and post-seismic ground motion maps ( Figure 2b) show a very smooth signal, which validates the robustness of the pixel-by-pixel approach.
In addition to this TS decomposition approach, we find, for each of the three tracks, the image pair with the shortest temporal baseline spanning the 10 November 2008 earthquake, and produce filtered, unwrapped and empirically corrected interferograms for each pair.

Fault Inference Method
Synthetic seismic waveforms and static surface displacements in a layered medium are computed with the Green's Functions (GFs) calculated using QSEIS and PSGRN/PSCMP, respectively [89,90]. For seismological data, the layering of the Earth is defined by the AK135 1D velocity model [91], while, for InSAR data, we use the elastic moduli from the local Global Crustal Database based on empirical 1D velocity models derived from seismic reflection and refraction profiles [92] ( Figure S6).
Fault optimisations are performed with a bootstrapping approach running in parallel 200 evaluations of the model misfit based on different realisations of random Bayesian weights (far-field) and random synthetic noise (InSAR) [4,81,82] (https://pyrocko.org/grond/docs/current/method/index.html). In a first exploration phase of the optimisation, the model space is randomly sampled. Afterwards, the exploration domain within the model space is defined by the model parameter covariances that are computed from a list of low misfit model realisations. Only one model in the best-fitting model list is allowed to be changed per iteration. Additionally, the center points for the sampling distributions are randomly offset from the mean value of the list to avoid becoming stuck in local minima. By doing so, multi-modal solutions are efficiently explored simultaneously.

Fault Inference of the 2008 Earthquake
We optimize one rectangular fault for the coseismic slip in the 2008 earthquake and a second fault for the postseismic deformation that followed. We explore the time offset of the earthquake in comparison to gCMT solution (Time), the position in comparison to the Global Centroid Moment Tensor (Global CMT) solution (northing, easting, depth), the fault dimension (length, width), the fault orientation (strike, dip, rake), the amount of slip (slip), the nucleation starting point (nucleation X, nucleation Y), and the rupture velocity (velocity) (Tables S2 and S3). We define the fault position (northing, easting, depth) as the center of the top edge of the rectangular fault, while the nucleation starting point is set relative to the plane center (0,0 is the center). In addition, to these earthquake fault model parameters, we also invert for a bilinear ramp in the InSAR data in east and north directions to tie the data to the model.
We perform a series of tests in order to explore the likelihoods of north-dipping and south-dipping planes. As a first test, we separately explore both north-dipping (U (250, 310) • and U (0, 80) • prior distributions for strike and dip angles, respectively) and south-dipping solutions (U (80, 140) • and U (0, 80) • prior distributions for strike and dip angles, respectively) in agreement with the co-seismic and post-seismic displacement maps and look for independent solutions that show coplanarity between the co-seismic and the post-seismic slip. In a second test, we invert for the stack of 2008 co-seismic displacement interferograms, that include both co-seismic and post-seismic deformation, leaving as free parameters the dip angle of the fault plane (U (80, 310) • and U (0, 180) • prior distributions for strike and dip angles, respectively).

Fault Inference of the 2009 Earthquake
We model the 2009 earthquake with three rectangular faults with only few and loose a priori constraints on the fault locations based on InSAR surface displacements because of the asymmetric displacement patterns that illuminate three distinct strike angles and based on the following back-projection analysis.
We model three distinct and independent fault segments at the same time, with a full set of fault parameters as described for the 2008 earthquake. We impose that no overlap between fault segments can occur (Table S4). Prior distributions of the individual fault positions are chosen, such that they roughly describe a central segment and one segment to the east and west.
We additionally explored both stacks of 2009 co-seismic and post-seismic interferograms together with the seismic waveforms within the same inversion framework as the co-seismic fault optimisation based on TS data in order to also assess the fault geometry and slip responsible responsible for both 2009 co-and post-seismic displacements. We here chose to not invert for post-seismic surface displacements extracted from the TS analysis as the model is too complex (three independent fault parameters) for the available data set (InSAR post-seismic data only).  In space, we observe, for the 2008 earthquake, a smooth southeast-northwest elongation of high semblance in both HF and LF coherent energy emissions (Figure 4b, Figure S7). The highest average HF semblance is achieved using a 16 km deep grid of points for the 2008 earthquake back-projection. In contrast, the 2009 rupture process appears to be more shallow (∼9 km) and more spatially focused (Figure 4b, Figure S8). It is reflected by three distinct peaks of HF energy, which suggest rupture of three separate steeply dipping ramps. The rupture of the 2009 earthquake starts in the center area and propagates bilaterally. In the eastern part, HF is emitted during an intermediate time step (3-6 s), in an area that roughly matches a gap in the geodetic displacement data (Figure 2b) and it may be associated with the jump of the rupture across fault-segments and the activation of a third fault segment. The central area shows some back-projected energy during the intermediate time step (3-6 s). The western area shows back-projection energy in the last stage of the rupture, associated to the activation with a third fault segment.

2008 Fault Characteristics Estimations
Among north-dipping solutions for the 2008 earthquake ( Figures 5 and 6a), the models that best fit jointly the InSAR TS and teleseismic data feature a 16 ± 2.5 km long and 3.6 ± 1.4 km wide segment with a top edge center located at 11.9 ± 0.5 km depth. The best-fitting solution converges toward a 264 ± 8 • -striking plane dipping at a low-angle of 32 ± 2 • towards the north. Waveform fits for the 2008 earthquake north-dipping solution for five randomly selected stations are shown in Figure S9.
The north-dipping fault model that best explains the first four months of post-seismic deformation following the 2008 event (imaged using TS data) is a narrow segment, 18.5 ± 1.4 km long and 0.6 ± 0.3 km wide, with a top edge center located at 12.5 ± 0.8 km depth. Slip is 1.2 ± 0.4 m for the co-seismic rupture and also of 1.2 ± 0.3 m for the four-month post-seismic period on the narrow fault patch. While the strike of the co-seismic solution is oriented parallel to the southernmost branch of the Olongbulak Shan, the post-seismic solution has a strike of 287 ± 2 • , perpendicular to the overall ∼N20 • E regional shortening [93,94]. The dip angle is also slightly lower than the co-seismic solution, dipping at 20 ± 1.3 • towards the north.
In comparison, south-dipping solutions ( Figures 5 and 6b) converge towards a 96 ± 4 • -striking co-seismic plane and a 102 ± 2 • -striking post-seismic plane with high-angle dips of 60 ± 2.5 • and 70 ± 1.5 • , respectively. Slip is of 1.7 ± 0.5 m for the co-seismic rupture and of 0.05 ± 0.003 m for the 4-month post-seismic period on the 29.7 ± 0.5 km wide and 16.5 ± 1.4 km long fault patch. As for the north-dipping solution, the south-dipping solution presents a very low misfit, as shown in the comparison between the geodetic data and the best-fitting model ( Figure S10). While the north-dipping co-and post-seismic solutions are almost co-planar, south-dipping solutions indicate two distinct planes that are located almost ∼5 km apart.
Additionally, from the joint probabilistic inversion of the near-field stack of co-seismic interferograms and the far-field waveforms data emitted by the 2008 earthquake (Figure 7), we converge towards a 14 ± 1.8 km long and 5.8 ± 2 km wide rupture plane with a top edge centre located at 9.8 ± 1 km depth. Although both north and south-dipping solutions, which are sampled during the whole optimisation (Figure 7b), show low misfit values due to the ambiguity if the displacement patterns of the buried fault plane, the best-fitting solution is a north-dipping plane striking at an angle of 286 ± 80 • and dipping at an angle of 32 ± 10 • . There is good agreement between both geodetic and seismological observations and predictions for the three independent tracks (Figure 7a).

2009 Fault Characteristics Estimations
For the joint probabilistic three-fault inversion of the 28 August 2009 M w 6.3 earthquake based on InSAR and far-field teleseismic data ( Figure 5), the optimisation converges towards a central segment of 10.0 ± 0.5 km long and of 1.6 ± 0.5 km wide with a top edge centre located at 2.7 ± 1.2 km depth. The eastern segment is constrained with a length of 4.7 ± 0.5 km, a width of 1.8 ± 0.5 km, and a depth of 3.1 ± 1.5 km, while the third western segment is of 9.0 ± 1 km length, 2.5 ± 1 km width, and at 4.5 ± 0.5 km depth. All three segments have narrow widths and are, respectively, dipping with 71 ± 2 • , 57 ± 3 • , and 73 ± 3 • high-angles to the south. Slip on the fault segments amounts to 2.0 ± 0.4 m, 2.2 ± 0.3 m, and 1.9 ± 0.5 m. A summary of all posterior PDFs is available in the Figure 8, while waveform fits for the 2009 earthquake north-dipping solution for five random stations are shown in Figure S9b.
Additionally, from the joint probabilistic inversion of the near-field stack of co-seismic interferograms and the far-field waveforms data emitted by the 2009 earthquake (Figure 8 and Figure S11), we converge towards a middle-segment of 9.4 ± 0.5 km long and 4.1 ± 1.2 km wide with a top centre edge located at 2.5 ± 0.3 km depth. The eastern eastern segment is constrained with a length of 6.1 ± 1.0 km, a width of 2.8 ± 0.9 km, and a depth of 2.5 ± 0.2 km, while the third western segment is of 6.8 ± 0.6 km-length, 2.8 ± 0.4 km-width, and 4.0 ± 0.5 km-depth. The three segments all dip steeply to the south, at 73 ± 3 • , 61 ± 4 • , and 59 ± 7 • , respectively. Therefore, fault geometries are very similar to the three fault segments obtained from the inversion of the co-seismic surface displacements only, and we can, therefore, reasonably conclude that post-seismic slip occurred on the same fault geometries as the coseismic rupture. Slip on the fault segments is 1.5 ± 0.6 m, 1.3 ± 0.2 m, and 1.9 ± 0.4 m, respectively.

Benefits of InSAR Time Series (TS) for Fault Inference
Traditionally, an earthquake is imaged with an interferogram, which contains the phase delay between two acquisitions. However, the travel time difference between the two acquisitions does not only contain the surface displacements associated with the earthquake, but also multiple undesired signals, such as ionospheric and atmospheric delays, orbital ramps, and topographic phase errors. The interferometric phase may also include some transient post-seismic surface displacements, depending on the time span between the acquisitions prior to and after an earthquake. Therefore, in comparison to single interferograms, modeled TS maps present lower variances and covariances (Table S5). The final data variances from the TS data range from 0.2 to 3×10 −6 m 2 , while the correlation lengths are about 0.5-2 km.
The stack of co-seismic interferograms also helps to improve the unwrapping coverage and the signal-to-noise ratio. Stacking is a simple way of removing turbulent atmospheric phase patterns that are random in time and space, and improves the spatial coverage, as interferograms used for the stack have varying unwrapping coverage, depending on their temporal and geometrical baselines and their respective noises. Conversely, TS data only have information for pixels where the TS analysis and the time decomposition is possible, ie. when there is enough unwrapped interferogram for a given pixel for the inversion to be reliable. However, and contrary to a TS approach, a stack may also constructively add tropospheric signal if it is of the same sign in all stacked interferograms and, therefore, may also contain more correlated noise than differential interferograms, as illustrated by the higher correlation length values of the stacked data in Table S5. In addition, by using long-temporal baseline interferograms that cover the post-seismic period within the stack, the stacking operation also increases the amount of post-seismic deformation within the data. This effect is illustrated in the 2008 fault inference constrained with stacked interferograms as posterior solutions are spatially located between co-seismic and post-seismic solutions inferred with the TS data (Figure 7a).
We perform another north-dipping 2008 co-seismic fault inference jointly constrained with the DInSAR interferograms shown in Figure 3 and with teleseismic data in order to further evaluate the benefits of the TS approach. We then compare the posterior PDFs with those derived from the joint inversion of co-seismic TS maps and teleseismic data (Figure 9). The analysis shows that, as for the solution constrained with the stack of co-seismic interferograms (Figure 7), the fault parameters inference constrained with DInSAR are biased towards a solution that averages the co-seismic and the post-seismic solutions derived from TS data ( Figure 6). The solution is pushed further north, where post-seismic slip occurred, at 10 ± 0.2 km of the gCMT solution, and towards a 18.2 ± 2.2 km long and narrow 1.7 ± 0.4 km wide rupture plane with a top located at 12.1 ± 0.5 km depth. Posterior PDFs of the depth indicate higher values for the DInSAR optimisations than the TS optimisations, which is also in agreement with the fact that post-seismic slip occurred at a greater depth than the co-seismic slip ( Figure 6). The best-fitting values for slip, strike, and rake are 2.0 ± 0.4 m, 284 ± 14 • , and 87 ± 15 • , and their PDFs are wider than those that were obtained from the inversion of TS and teleseismic data. Slip for DInSAR is also expected to be higher due to the presence of early afterslip in the data. The DInSAR best-fitting 22 ± 4 • northward dip angle is ∼10 • lower than the dip angle reported in the gCMT catalog and the one of the northward-dip constrained by TS solutions, but again closer to the dip solution obtained from the post-seismic fault ( Figure 6). Posterior PDFs of the easting, depth, slip, strike, dip, and rake parameters have higher standard deviations when constrained from DInSAR than the TS maps, but, inversely, posterior PDFs for the northing and width fault parameters are narrower with DInSAR. This comparison shows that the inclusion of afterslip in the co-seismic modelling may introduce some large biases in the results towards the post-seismic solution. Those biases are proportional to the magnitude of the aseismic slip and the period between the two acquisitions of the co-seismic interferogram. Additionally, higher posterior uncertainties in the fault characteristics estimation are expected for the DInSAR optimisation in comparison to TS due to the higher noise in the co-seismic interferograms. However, wider PDFs for DInSAR is not, in this case, clearly highlighted for all parameters, which may be due to efficient empirical corrections performed on the DInSAR interferograms.
A pixel-by-pixel decomposition of the cumulative InSAR displacement TS has the limitation of providing a solution that may not have a spatial physical meaning and with reliability that depends on the number of acquisition dates and/or on the residual noise in the TS data. A multi-pixel TS decomposition approach [95] may tackle part of these limitations by accounting for the spatially correlated noise in the surface displacement data. However, the spatial consistency that was obtained in the decomposition (Figure 2) validates the pixel-by-pixel method and the reliability of the results. In addition, the simplicity of the inverse problem (contrary to a multi-pixel approach that demands the inversion of a large matrix function of the number of InSAR points) allows us to perform optimisation with some additional constraints, such as imposing post-seismic surface displacements of the same sign than co-seismic surface displacements. An alternative approach would be a decomposition of the InSAR time-series data by joint and direct inversion of slip on the co-seismic and post-seismic fault planes. Such an approach would impose a physical meaning of the decomposition contrary to a pixel-by-pixel or multi-pixel decomposition in time. However, it also requires a fixed fault geometry e.g., [30,43,48,50] as the joint optimisation of both fault source parameters and temporal basis functions for each InSAR point is computationally expensive and such a modelling tool is currently not available to our knowledge. Figure 9. Comparison between posterior Probability Density Functions (PDFs) of the 2008 earthquake parameters that were obtained from the optimisation of InSAR co-seismic TS data + teleseismic data (blue) and DInSAR interferograms + teleseismic data (coral red). Dashed vertical lines are best-fitting models.

Back-Projection Imaging for Moderate-Size Earthquakes
The derived 2008 earthquake's beam-power time function (Figure 4a) is very similar to the source time function from the SCARDEC database [96] (http://scardec.projects.sismo.ipgp.fr/). Of particular note, the method retrieves information about the source of the 2009 earthquake even if some waveforms might be contaminated by the 2009 Banda Sea Earthquake. Therefore, the method shows its efficiency in comparison to other approaches to suppress unwanted signals arriving from other azimuths through destructive inferences [67]. It is robust against signals from outside the defined source grid due to the stacking procedure.
The back-projection tool has also located the coherent emissions of both earthquakes with a remarkable agreement with geodetic data. It also shows, with the 2009 earthquake, promise as a method to infer the level of segmentation of a rupture without formal kinematic modelling, even for moderately sized earthquakes.

The 2008-2009 Qaidam Earthquake Sequence
By reconciling co-and post-seismic geodetic observations, back-projection imaging, kinematic fault inferences, and previous structural interpretations, we develop a model for the fault geometry of the Qaidam thrust system shown in Figure 10, with overall fault and fold geometry typical of foreland fold-and-thrust belts Daout et al., 2019 [45] modeled, from Envisat measurements across the Zolonbulak Shan, the short-term after-slip following the 17 April 2003 Delingha earthquake with a shallow segmented ramp-flat-ramp north-dipping structure that steepens as it approaches the surface under the Delingha anticline, and inferred a low-angle north-dipping fault for the 2003 earthquake (orange line in Figure 10). The north-dipping gCMT solution of the 2003 earthquake is very similar to the 2008 gCMT solution with a ∼30 • dip angle rupture plane at about 16 km depth (Table S1, red focal mechanism in Figure 10). In addition, the two ∼18 km and ∼21 km deep April 2004 right-lateral ∼M w 5.3 earthquakes and the ∼12 km deep October 2004 thrust ∼M w 5.5 earthquake (according to gCMT, Table S1) may belong to the same flower structure that ruptured in 2003 ( Figure 10). Static stress transfer within this flower structure may have triggered the 2004 sequence of earthquakes one year after the 2003 earthquake.
This study show that kinematic inversion alone cannot rigorously resolve the ambiguity between the north-dipping and south-dipping geometries due to the low misfits of both solutions. Biases may also arise from assumptions on the 1D crustal structure or the absence of topography in the model, although bootstrapping technique can efficiently estimate model parameter uncertainties by propagating data errors and by assessing modelling errors. In addition, the stations have individual time shifts to compensate for crustal structures, such as each station can shift its waveforms up to 4 s with respect to the earthquakes' onset time to achieve maximum fit. This compensates for unmodeled three-dimensional (3D)-path effects.
Here, we infer a shallow north-dipping low-angle structure that is responsible for the 2008 earthquake (green line in Figure 10), mainly constrained from the coplanarity between the co-and post-seismic north-dipping solutions ( Figure 5). The depth of the rupture (∼12 km deep) is in agreement with the width of the Olongbukak range (∼10-15 km), as north-dipping and south-dipping faults bounding the Olongbulak range are expected to merge as those depths ( Figure 5). The presence in the morphology of a major north-dipping structure in the southern front of the Olongbulak range also challenges a hypothetic south-dipping and high-angle 2008 rupture plane, coplanar to the 2009 rupture ( Figures S1b and S12). Such coplanar fault geometries would not create a differential vertical movement between the deeper part of the south-dipping fault (where the 2008 earthquake occurred) and the shallow part of the fault, within the Olongbulak Shan (where the 2009 earthquake occurred) [19] ( Figure S12b). The absence of differential vertical movements would, therefore, require no triple junction and shallow high-angle north-dipping plane, while this structure is identified in the southern front of the Olongbulak Shan. As a triple junction is not stable over long time scale and, as the current morphology involves multiples deformation stages, it is difficult to conclude from this observation. However, the main north-dipping vergence of the NQT and the current activity of its southernmost branches, as within the Olongbulak pop-up, support our view of a low-angle north-dipping plane, responsible for the 2003 and 2008 earthquakes, branching at the roots of the shallow high-angle thrusts of the NQ ranges and of the Qaidam basin.
Our results also indicate that the 2008 earthquake did not break the southern and north-dipping high-angle branch of the Olongbulak pop-up structure. Post-seismic slip is constrained down-dip of the co-seismic rupture on a very narrow plane with slip of the same order of magnitude than the co-seismic slip. Surface displacements might not be very sensitive to the amount of down-dip afterslip, and the amount of slip might, therefore, be over-estimated for the north-dipping solution. In this regard, we performed several tests forcing a higher fault width or lower fault slip. We observed that with such prior values, posterior PDFs tend to be asymmetric towards the boundary values, indicating that the data require such high slip and narrow fault width. The modelling, therefore, shows that the observed short-term post-seismic displacements cannot be explained by fault movement on the north-dipping shallow and high-angle branch of the Olongbulak pop-up structure, thus leaving an open question on the seismic potential of this branch of the Olongbulak range.
The 2009 rupture is the result of the vertical partitioning of the shortening along the high-angle and shallow south-dipping back-thrust of the Olongbulak Shan (blue line in Figure 10). The rupture started at ∼9 km depth along the middle segment and then propagated bilaterally along two shallower and higher-angle segments (Figures 4 and 5) that displace basin-deformed sediments of the Olongbulak range above the Paleozoic bedrock of the Zonlongbulak range [47,57,97]. Afterslip likely occurred on a similar fault patch surrounding the coseismic slip, in agreement with similar post-seismic and co-seismic surface displacement patterns measured with InSAR (Figures 2 and 3), the optimisation of the stack of long-and short-temporal baselines interferograms ( Figure S11), and previous studies [48][49][50]. The observed rupture complexity in both geodetic displacements and back-projected seismic emissions likely reflect the along strike and along depth lithological segmentation of the Olongbulak fault system.
The 2008-2009 Qaidam earthquake sequence illuminates a short-timescale and spatially confined example of regional tectonic and geodynamic processes at play. However, the inferred flower structure, in which high-angle thrusts and back-thrusts root into a low-angle décollement, may be typical of structures and deformation styles in the South Qilian Shan at greater scale. The low-angle north-dipping structure may connect to the roots of the Olongbulak pop-up structure, the Zongbulak Shan, and the Xietie Shan thrusts and may transfer the shortening rate from the north, along the Qilian Shan, to south, in Qaidam Basin, along the high-angle folds and thrust belts ( Figure 10, inset, Figure S1a). It supports regional tectonic models that involve a regional décollement that connects to the Qilian Shan suture or pre-existing Cenozoic lithospheric weaknesses further north [52,54,57,59] and the southward expansion of the South Qilian Shan [60,97].

Conclusions
In this study, we explore the benefits of an InSAR TS approach for earthquake fault inference and demonstrate the improved signal-to-noise ratio of such data sets in comparison to DInSAR, and the additional constraints offered by the post-seismic surface displacements that were extracted from the InSAR TS data. The 2008-2009 Qaidam earthquake sequence exemplifies typical earthquakes in remote areas, where field work is difficult and structural seismic profile or geological models are rarely available or poorly constrained, yet satellite-based geodetic and teleseismic data are easily available.
We use the open source Grond toolbox [81], the Pyrocko framework [98] and the Kite toolbox [40] for pre-processing and analysing surface displacements to infer the earthquake fault parameters from both geodetic and teleseismic data. Pyrocko provides a uniform framework for both near-field and far-field data and Green's function databases to speed up forward calculations. Our analysis reveals that the optimisation scheme successfully explores bimodal north-dipping and south-dipping solutions for the deep 2008 earthquake, which was further validated based on the co-planarity between the inferred co-seismic and post-seismic fault planes. It also shows that the tool allows the exploration of the rupture segmentations with three independent faults, where no fault geometrical parameters are fixed a priori. Besides, the back-projection tool, also based on the Pyrocko framework, has located the coherent emissions of both earthquakes with a remarkable agreement with geodetic data. It also constrains the segmentation of 2009 ruptures in accordance with the InSAR-based fault models and, therefore, shows promise as a useful tool, even for moderate earthquakes.
Supplementary Materials: Figure S1: Simplified sketches summarising the two opposite geodynamic deformation models of northeastern Tibet with their implications for the faults geometry of the North Qaidam thrust system, Figure S2: Station array locations used for the backprojection of the 2008 and 2009 earthquakes. Blue lines are distances in degree and red lines are the inner and outer circles of the station selection at 22 and 94, Figure S3: Empirical time shifts at virtual array stations for the P-phase for a) low-frequency and b) high-frequency backprojections that maximise the semblance of the reference event and are used for the backprojection, Figure S4: Computed interferograms for the three tracks. Triangles are SAR acquisitions with sizes according the their spatial extent, Figure S5: Time series of surface displacements from 2008 to 2011 for the pixel 2 of track 319 of Figure 2 (blue circles) with best-fitting estimations of long-term velocities, 2008, 2009 co-seismic offsets and logarithmic afterslip functions for three relaxations times (1 day in red, 15 days in blue and 45 days in green), Figure S6: Crustal velocity model of the elastic stratified medium used to compute the near-field surface displacements, Figure S7: Semblance mapping for each time-steps for the 2008 Qaidam earthquake, Figure S8: Semblance mapping for each time-step for the 2009 Qaidam earthquake, Figure S9: Waveform fits of the P-phase for the 2008 earthquake north-dipping solution (a) and the 2009 earthquake (b) for five random stations, which are common to both datasets. Restituted and filtered traces without tapering are in light grey while traces with tapering and processing are in dark grey, Figure S10: Comparison between data and model from the optimisation of one rectangular southdipping fault in agreement with the 10th November 2008 earthquake data, Figure S11: Posterior models for the 28th August 2009 earthquake obtained from the optimisation of three rectangular faults in agreement with a stack of co-seismic interferograms, Figure S12: Conservation of the kinematic motion across the Olongbulak pop-up (OLT) for the south-dipping scenario, Table S1: Summary of the gCMT solutions for the 8 Mw >5.2 earthquakes that occurred in Qaidam between 2003 and 2009, Table S2: Summary of the prior probabilities for the 2008 co-seismic rectangular fault inference U defines normal distribution, Table S3: Summary of the prior probabilities for the 2008 post-seismic rectangular fault inference, Table S4: Summary of the prior probabilities for the 2009 co-seismic rectangular fault inference, Table S5: Comparison of the variance-covariance estimations of the InSAR co-seismic Time Series (TS) data maps, co-seismic interferograms (IFG), and stack of co-seismic Interferograms are available online at http://www.mdpi.com/2072-4292/12/17/2850/s1.