# The interplay between seismic and aseismic slip along the Chaman fault illuminated by InSAR

M. Dalaison<fnref n="1" />, R. Jolivet<fnref n="1" />,<fnref n="2" />, E. M. van Rijsingen<fnref n="1" />,<fnref n="3" />, S. Michel<fnref n="1" />

## Key Points:

* Aseismic slip along the Chaman fault reaches 12 mm/yr with three large distinct aseismic sections, despite marked elastic strain increase.

* Three earthquakes are imaged along the central Chaman fault and two exhibit large post-seismic afterslip.

* Along-strike distribution of slip is compared with historical seismicity and fault geometry.

<footnote n="1">¹École Normale Supérieure, PSL University, CNRS U.M.R. 8538, Laboratoire de Géologie, Paris, France</footnote>

<footnote n="2">²Institut Universitaire de France, 1 rue Descartes, 75005 Paris, France</footnote>

<footnote n="3">³Department of Earth Sciences, Utrecht University, Utrecht, the Netherlands</footnote>
Abstract

The 700-km-long Chaman fault marks the western edge of the plate boundary between India and Eurasia. Although global plate models predict 2.3-3.6 cm/yr left-lateral motion between both plates, the fault is known to have hosted few earthquakes in historical times. Recent geodetic measurements attested the presence of aseismic slip locally. To detail the interplay between fast and slow slip along the Chaman fault, we build three InSAR time series of ground deformation covering the whole fault length over 5 years (2014-2019). We find that most of the active fault trace slips aseismically and continuously. From south to north, we identify three creeping fault portions: the Nushki, Central and Qalat segments of lengths between 80 and 130 km. The loading rate is 1.2±0.3 cm/yr for the two southernmost portions, while it is about 0.7±0.2 cm/yr for the Qalat segment. The Central segment and the nearby locked segments have hosted the largest known historical earthquakes on the Chaman fault, and three moderate magnitude earthquakes in our observation period. We image these earthquakes for which modeled slip at depth (Mw5-5.6), time series of surface slip and deformation patterns argue towards large triggered aseismic slip. The June 2018 event displays postseismic moment 3-15 times greater than coseismic moment. Over the two decades covered by geodetic observations, continuous or triggered aseismic slip dominates along most of the fault and co-locates with earthquakes. We observe that fault geometrical complexities delimit active segments and may be responsible for the kilometer-scale intertwining between seismic and aseismic events.

Plain Language Summary

The Chaman fault separates Indian and Eurasian tectonic plates moving at about 3 cm/year with respect to each other. This fault is known to have hosted very few earthquakes (sudden and rapid slip) in historical times, and appears to slip slowly. We measure slip on the Chaman fault with radar satellite images covering 2014-2019, and assess the average amount of slip attributed to earthquakes over the past century. We identify three 80-130-km-long portions of the 700-km-long fault that slip silently, without radiating seismic waves, at rates reaching 1.2±0.3 cm/year along the southern half, and 0.7±0.2 cm/year in the north. The largest historical earthquakes and three earthquakes of moderate size in our observation period occurred on the central part of the Chaman fault. We image the surface deformation induced by these three earthquakes and find that they were followed by significant slow slip on the fault, especially the third event in June 2018. Since the 2000s space-based measurements indicate that the slow mode of slip dominates on the Chaman fault plane and co-locates with earthquakes. Changes of fault orientation may delimit portions of the fault which either slip slowly or break into fast earthquakes.

# 1 Introduction

The potential occurrence of earthquakes depends on the available elastic energy stored in the crust, which in turn depends on external driving forces and fault properties. While some faults are locked most of the time and release energy by rapid (m/s) slip (i.e. earthquakes), some slip slowly (cm/yr) [e.g. Steinbrugge et al., 1960, Jolivet, Simons, et al., 2015, Bürgmann, 2018, Jolivet & Frank, 2020]. Slow slip is known to influence the budget of slip locally [e.g. Çakir et al., 2012, Maurer & Johnson, 2014, Jolivet, Simons, et al., 2015, Dal Zilio et al., 2020], hence plays a role in tuning the magnitude [e.g. Michel et al., 2018, 2019], as well as the initiation, propagation and arrest of potential earthquakes [e.g. Kaneko et al., 2010, Avouac, 2015]. Understanding the interplay between slow slip, shear loading and earthquakes is therefore fundamental for seismic hazard assessment.
Because slow slip does not radiate seismic waves, unlike earthquakes, global occurrences of slow slip remain incompletely documented and their ubiquity is an open question. In this study, we focus on the Chaman fault, in Pakistan and Afghanistan, where slow slip has been described [Fattahi & Amelung, 2016, Barnhart, 2017] alongside multiple moderate magnitude earthquakes (moment magnitude, M<sub>w</sub>, of ~5) [Ambraseys & Bilham, 2003a, Fattahi et al., 2015, Bilham et al., 2019]. This fault belongs to one of the longest continental strike-slip systems extending over a length of 1000 km and accommodates the relative displacement between India and the Helmand Block of the Eurasian Plate [Yeats, 2012] (Figure 1). However, a detailed characterization of the spatio-temporal dynamics of slip along the full length of the Chaman fault (CF) is lacking. In the following paragraph, we briefly summarize the current understanding of the kinematics of the CF and precise the goal of this study.

Global plate models predict that the Chaman fault system, including the CF together with the adjacent fold and fault belts to the east, acts as a transpressional plate boundary and accommodates left-lateral motion ranging from 23 to 36 mm/yr and convergence between 0 and 18 mm/yr [Ruleman et al., 2007, DeMets et al., 2010, Kreemer et al., 2014, Altamimi et al., 2017]. Global Positioning System (GPS) and Interferometric Synthetic-Aperture Radar (InSAR) suggest left-lateral motion is currently accommodated by the CF at about 12 mm/yr with significant along-strike variations [Mohadjer et al., 2010, Szeliga et al., 2012, Crupa et al., 2017, Fattahi & Amelung, 2016, Barnhart, 2017]. Consequently, it is currently accepted that the CF accommodates about 30% of the total relative plate motion between latitudes 27°N and 35°N [Fattahi & Amelung, 2016, Barnhart, 2017]. Remaining relative displacement is probably accommodated by structures east of the CF, such as the Ghazaband fault and the Gardez fault zone [Szeliga et al., 2012, Fattahi & Amelung, 2016, Bilham et al., 2019] (Figure 1).

Interestingly, some observations along this major plate boundary challenge the classic elastic cycle of loading and unloading of a single fault through earthquakes. Indeed, a long-known characteristic of the Chaman fault system is the scarcity of seismicity, which has been reliably documented in the region since the 1880s [Ambraseys & Bilham, 2003a], and, specifically, the very few earthquakes with magnitude greater than 6, which include the M<sub>w</sub>6.6 Chaman (1892) and M<sub>w</sub>7.7 Quetta (1935) earthquakes [Dewey et al., 2006, Ambraseys & Bilham, 2003b,a, Ambraseys & Douglas, 2004, Szeliga et al., 2009, Bilham et al., 2019]. Along the same fault system further north, the city of Kabul was destroyed by an earthquake of magnitude ~7.2 in 1505 [Ambraseys & Bilham, 2003a, Yeats, 2012]. So few large magnitude events cannot explain the relative motion accommodated by the CF and even less the whole plate boundary relative motion [Bernard et al., 2000]. Common explanations of the apparent lack of large earthquakes include an important role of (i) aseismic slip and microseismicity below the detection threshold of the local bare seismic network or of (ii) non-rigid deformation within the wide fault system. Distributed deformation on many structures (in support of (ii)) is suggested by detailed fault traces inferred from satellite imagery and current seismicity [Ruleman et al., 2007], both outlining intricate fault patterns with many subsidiary faults and discontinuities over a >80-km-wide region.

Here, we evaluate the role of aseismic slip and microseismicity on strain release along the CF (hypothesis (i)), focusing on the interplay between seismic and aseismic slip. A 125-km-long creeping segment, called the Nushki creeping segment, has been identified along the CF between latitude 29.28°N and 30.58°N [Fattahi & Amelung, 2016, Barnhart, 2017]. Additionally, the inferred shallow locking depth (0-7 km) along the CF is thought to impede the occurrence of large earthquakes [Szeliga et al., 2012, Barnhart, 2017]. Moreover, the abnormally large and long (1 year) afterslip following at least two M<sub>w</sub>5 earthquakes [Furuya & Satyabala, 2008, Fattahi et al., 2015] and the recent creep event presented in Bilham et al. [2019] highlight the need to frame observations temporally, and consider the role of transients in long-term behavior. Detailed time series are needed to
untangle the dynamics of the fault system and depict a better image of the fault frictional properties. We therefore propose to examine the distribution of aseismic slip along the whole CF length, using high resolution time series of deformation derived from InSAR. We cover the 2014-2019 period, building on existing geodetic observations for older periods.

We consider more than 5 years of Synthetic-Aperture Radar (SAR) acquisitions (2014-2019) on ascending and descending orbits by the Sentinel 1 constellation to characterize spatio-temporal variations of aseismic slip and investigate the recent distribution of seismicity. In the following, we first summarize the structural and geological properties of the CF and detail how we select earthquakes in existing catalogs of seismicity. Then, we detail an innovative systematic data assimilation approach to process InSAR data and deal with five moderate magnitude earthquakes that occurred during our observation period. Next, we carefully extract slip and slip rate along the fault trace from InSAR-derived displacement and velocity maps. Also, through inversion of InSAR-derived earthquake surface deformation, we model slip at depth for three earthquakes affecting the CF and study the spatio-temporal signature of the resulting surface slip. Finally, we describe the intriguing interplay between earthquakes and aseismic slip and discuss a segmentation of the CF, relating slip behavior with fault geometry.

## 2 Tectonic context and seismicity

### 2.1 Activity and geometry of the fault

In the topography, the CF sharply delimits the western border of young mountain ranges (Eocene to now) from 28° to 32.5°N: the Kirthar range elongated along a NNE axis and the lobate Sulaiman range north of 30.5°N. With a more subtle topographic and geomorphologic signature, the fault extends north to the Afghan capital, Kabul (~34.5°N) [Crupa et al., 2017, Wheeler et al., 2005]. The CF cuts through structural trends and is a major geological boundary between the Miocene-Eocene Khojak flysch unit to the east and the Afghan block essentially made up of meta-igneous, volcanic and plutonic rocks [Jones et al., 1961, Lawrence & Yeats, 1979, Lawrence et al., 1981]. The CF is a long-lasting feature that has accommodated several hundred kilometers of sinistral offset and continues till present time as attested by young stream beds offsets and recent earthquakes [Beun et al., 1979, Tapponnier et al., 1981, Sborshchikov et al., 1981, Lawrence et al., 1992, Wheeler et al., 2005, Ul-Hadi et al., 2013]. Considering such cumulative offset, the CF reached significant level of structural maturity in the sense of Manighetti et al. [2007].

Precise fault trace mapping for the Chaman fault zone north of 29°N has been done by the U.S. Geological Survey using remote-sensing imagery [Ruleman et al., 2007]. The northeast-trending left-lateral fault zone is associated with smaller-scale thrust faulting indicating transpressional motion (Figure 1). From south to north, the main fault trace progressively rotates eastward from 15 to 37° azimuth and displays a 100-km-long restraining and releasing bend in its central portion, close to the city of Chaman (30.9°N). Compressive features mostly focus west of the CF within and north of this restraining bend. East of the CF, the Kirthar and Sulaiman fold and thrust belts host thrust earthquakes and are most likely accommodating a significant part of the 0-18 mm/yr of compression expected from the CF azimuth. The CF geometry at depth is poorly known but according to the parameter exploration by Barnhart [2017] the best hypothesis remains that the fault plane is vertical with negligible fault-normal displacement at least between 29° and 30.6°N.

As reference we draw a unique and continuous Chaman fault trace from Ruleman et al. [2007]'s mapping, following exact fault coordinates. The assumption of a continuous fault trace is consistent with fault mapping from 28° to 32.5°N, while mapped faults north of 32.5°N tend to be more discontinuous. However, the uniqueness of the fault trace is a
challenge as several fault strands are mapped often within 1 km of each other. To get
a fine scale quantification of azimuth variation we differentiate the fault trace using Char-
trand [2011]'s algorithm. We will discuss fault azimuth variations and their relationship
with fault slip at a scale larger than 10 km with specific caution north of 32.5°N in Sec-
tion 5.3.

## 2.2 Seismic record

The Chaman fault system seismicity is scarce and distributed over the Kirthar moun-
tain range (Figure 1). We explore the seismic record with a double interest: first, we pre-
cise their long-term (~100 years) contribution to fault slip along the CF, second, we eval-
uate their contribution to the observed deformation measured by InSAR over 2014-2019.
We compile available seismic events in the Harvard Global Centroid Moment Tensor (GCMT),
U.S. Geological Survey (USGS) and International Seismological Center (ISC) catalogs
[Dziewonski et al., 1981, Ekström et al., 2012, ISC On-line Bulletin, 2020]. Because seis-
mic stations are rare and far away, local earthquakes are poorly located with uncertain
depths and magnitudes. For instance, the 2005 M<sub>w</sub>5 earthquake imaged by Furuya &
Satyabala [2008] is located about 27 km (ISC) and 30 km (USGS) away from its actual
epicenter relocated by InSAR. We take the ISC as a reference for its re-location and com-
pilation of numerous magnitude estimates [Bondár & Storchak, 2011]. We use M<sub>w</sub> es-
timates, when available, or transform m<sub>b</sub> to M<sub>w</sub> using the global linear relation by Scordilis
[2006]. In order to compare the amount of slip released seismically to aseismic slip on
the CF, we collect events since 1900 occurring within 30 km of the Chaman fault trace.
As a result, our catalog going back to 1900 includes 139 events of M<sub>w</sub> between 3.8 and
6.5.

Among these events, some must affect surface displacements measured by InSAR, hence
must be accounted for when inferring surface deformation rates. Therefore, we identify
earthquakes in the spatio-temporal frame of the Sentinel 1 acquisitions (blue frames in
Figure 1) that are likely to influence inferred surface displacements. Identified events must
have a focal depth inferior to 25 km, an estimated M<sub>w</sub> superior to 5 and have to be dis-
tinguishable from each other given the InSAR temporal sampling of 6-12 days. Out of
the 16 and 170 events in the GCMT and ISC catalogs, respectively, we identify four events
or groups of events corresponding to our criteria with M<sub>w</sub> between 5.1 and 5.7 (Table
S1). In addition, we consider the earthquake on 27 June 2018 that clearly shows within
our InSAR data, even though it is not referenced in the GCMT catalog and the ISC database
indicates a m<sub>b</sub> between 3.8 and 4.2, depending on sources (Figure S1).

Therefore, we consider the deformation related to five earthquakes or group of earthquakes
during our observation period taking the properties (time, location and magnitude) of
the largest earthquake of the group. Each group contains earthquakes within 20 days
and 40 km of the largest one, considering uncertainties in location and timing arising from
seismic catalogs and from InSAR time series, respectively. Earthquakes sequences will
be detailed in Section 4.3.1. Among these five events, the two earliest events (3 August
2015 and 21 March 2016) are off the CF towards the southeast, whereas the three most
recent events on 13 May 2016, 10 July 2016 and 27 June 2018 occurred close to the CF
near the town of Chaman (Figure 1).
<figure><img src="image_1.png" /></figure>

**Figure 1.** Tectonic setting and seismicity of the Chaman Fault system. Focal mechanisms are from the GCMT catalogue [Dziewonski et al., 1981, Ekström et al., 2012]. Earthquakes later than 2014 are in red. Gray dots are historical and recent (2000 BC to 2004 AC) seismic events for which magnitude estimates are known (radius is relative magnitudes) [Dewey et al., 2006]. Orange lines are fault traces [Ruleman et al., 2007]. Background shading is elevation gradient from SRTM [Farr et al., 2007].
## 3 Method

In this section, we describe how we use SAR images to obtain a picture of tectonic deformation in the Chaman fault zone, and more specifically along the Chaman fault (CF) itself. First, we detail the way SAR images are combined to produce a network of interferograms (Section 3.1). Second, we explain how we derive maps of ground deformation, velocities and associated uncertainties with our innovative and efficient time series analysis method [Dalaison & Jolivet, 2020] (Section 3.2). Third, we explain how we combine these maps to obtain a 2-D deformation field and to measure surface slip along the CF (Sections 3.3 and 3.4). Fourth, we depict the inversion procedure used to get an image of slip at depth for three earthquakes which occurred on the CF during our observation time span (Sections 3.5).

### 3.1 Interferometric processing

We process all available Sentinel 1A-B wide swath SAR images on two ascending (tracks 42 and 144) and one descending (track 151) tracks until December 31st 2019. The viewing geometry of track 151 makes it nearly insensitive to fault parallel displacement (Figure 1) but this second direction of measurement is key to differentiate horizontal from vertical displacement. We use the ISCE package (JPL/Caltech, winsar.unavco.org/isce.html) to corregister SAR images, as well as to compute, multilook, filter and unwrap interferograms [Gurrola et al., 2010]. We build 456, 482 and 424 interferograms out of the 118, 120 and 108 acquisitions over tracks 42, 144 and 151, respectively. Interferograms are built systematically between each acquisition and the next 4 acquisitions. In addition, we build interferograms with longer temporal baselines (several months to years) for track 144 and 151 to tackle the loss of coherence over winter in the northern, mountainous part of the area. Interferometric pairs are shown in Figure S2. We apply multilooking for a final pixel size of 85 m along the satellite flight direction and of 62 m in the range direction. We apply a phase preserving Gaussian filter [Goldstein & Werner, 1998] on the interferograms and correct the interferometric phase from tropospheric delays using the ERA-5 global atmospheric reanalysis [PyAPS; Jolivet et al., 2011, Jolivet, Agram, et al., 2014]. We unwrap areas with a minimum coherence of 0.6 using a branch cut algorithm [Goldstein et al., 1988] and subsequently correct potential unwrapping errors [CorPhU; Benoit et al., 2020]. To correct for residual orbital errors, we remove a bilinear polynomial ramp from each interferogram. This ramp corresponds to the best fitting ramp of the phase change north of the CF, a region considered as tectonically stable (Figure S3). The end result of this processing is a network of unwrapped interferograms co-referenced to a small zone of 5 by 5 pixels considered as fixed. Each interferogram is the relative phase change between two dates in line of sight (LOS) converted to millimeters, thus, it contains the effect of ground deformation but also significant residual atmospheric delays.

### 3.2 Iterative time series analysis incorporating coseismic displacements

We apply a Kalman Filter time series analysis (KFTS) [Dalaison & Jolivet, 2020] on each track to reconstruct the evolution of the interferometric phase ($\phi_k$) at each time step ($t_k$) together with the parameters of a time dependent model and associated uncertainties. Unwrapped interferograms are being ingested sequentially to build the temporal evolution of the phase and refine the parametric model at each time step. The model is used to make a forecast at future time steps, which is then refined by including new interferograms. This data assimilation approach is comparable to the New Small Baseline Subset method as it builds a model of deformation used when data is missing [Berardino et al., 2002, Doin et al., 2011]. However, our method has the benefit of allowing fast and efficient update of pre-existing time-series as new interferograms become available, as well as providing detailed error propagation throughout the process.
The parametric model describes the evolution of deformation with a linear combination of predefined functions of time from which we optimize the L coefficients $a_i$ ($0 \le i < L$). A well designed model conditions the accuracy of the estimated mean ground velocity [Dalaison & Jolivet, 2020]. Notably, we must account for instantaneous phase changes caused by earthquakes.

Over the observation period, we identify five earthquakes (sometimes associated with smaller nearby earthquakes) from the GCMT and ISC catalogs likely to produce significant signal in the time series of phase change (Section 2.2). With this a priori information, we include a Heaviside step function of time ($H_j$) centered on the date of each earthquake in our parametrized model. Consequently, our parametrized model of deformation includes a constant term ($a_0$), a secular rate (velocity $a_1$), a seasonal oscillation modeled as the sum of a sine and cosine with a period, $T_{yr}$, of one year, and N Heaviside functions $H_j$ centered at the time of the N earthquakes affecting the area (here $N = 5$). Therefore we optimize $a_i$ for all i in 0 to N+3 so that

$$ \phi_k = a_0 + a_1 t_k + a_2 \sin \left( t_k \frac{2\pi}{T_{yr}} \right) + a_3 \cos \left( t_k \frac{2\pi}{T_{yr}} \right) + \sum_{j=4}^{3+N} a_j H_j(t_k) + \gamma_k \quad (1) $$

where $\gamma_k$ represents the mismodeling error. The associated a priori standard deviations are 25 mm for $a_0$, 9 mm/yr for $a_1$ and 8 mm for $a_2$ and $a_3$. Those values, associated with a null a priori, are chosen in such a way that they reflect the expected range of variation [Dalaison & Jolivet, 2020].

Because the spatial extent of earthquake related deformation is limited, including a Heaviside function to pixels far from the epicenter might lead to over-fitting. Therefore, to limit the spatial extent of the Heaviside functions, $H_j$, we force the earthquake amplitude $a_j$ to stay at zero far from the earthquake location by setting a null a priori on $a_j$ with zero uncertainty. In practice, we defined the a priori variance of $a_j$ for each pixel as a two-dimensional Gaussian function centered on the earthquake location so that the variance decreases smoothly as a function of the Euclidean distance to the earthquake (Text S1).

As the standard deviation of the null a priori on $a_j$ limits the explored range of values, it should be greater or equal to the order of magnitude expected for $a_j$ [Dalaison & Jolivet, 2020]. In agreement with the amount of displacement and size of the rupture predicted by empirical scaling laws for $M_w \le 5.6$ [Wells & Coppersmith, 1994, Petersen et al., 2011], we choose to set the maximum a priori standard deviation for $a_j$ to 30 mm and the characteristic width of the spatial Gaussian to 9 km (Text S1 and Figure S4). This implies that we look for earthquake related displacement within a radius of ~40 km, a region much larger than the typical area affected by shallow $M_w \le 5.6$ earthquakes [Savage & Burford, 1973, Petersen et al., 2011] so that the model does not put too much weight on imperfect earthquake location.

KFTS includes two additional tunable parameters. First, the standard deviation of the disclosure error, $\sigma_e$, measures interferometric disclosure, which originates from interferometric processing, mainly multilooking [De Zan et al., 2015] since unwrapping errors are corrected for. Second, the standard deviation of the mismodeling error, $\sigma_\gamma$, quantifies the difference between phase change evolution and the parametric time dependent model, and, thus, accounts for temporally decorrelated signal in InSAR thought to arise from turbulent tropospheric delays as well as for unmodeled processes. Effectively, it controls the distribution of $\gamma_k$ in Equation 1. Following guidelines in Dalaison & Jolivet [2020] we set $\sigma_e$ to 0.02 mm and $\sigma_\gamma$ to 10 mm. The fact that $\sigma_e \ll \sigma_\gamma$ guarantees that interferometric phase reconstruction from interferograms is prioritized over fitting the much more uncertain model.
### 3.3 Combining InSAR tracks

From the time series analysis on each of the three InSAR tracks we obtain a map of ground velocity in the line of sight (LOS) (i.e. parameter $a_1$ for each pixel of each track). Where both ascending tracks overlap, the difference in velocity estimates is minimized subtracting a bilinear ramp to the velocity field of track 144. We combine ascending and descending LOS velocities, to get vertical ($v_V$) and fault-parallel horizontal velocities ($v_H$), assuming zero fault-perpendicular velocity. The observed velocities ($a_1$)$_n$ relates to $v_H$ and $v_V$, for every $n^{th}$ InSAR track imaging the given pixel [e.g. Wright et al., 2004, Tymofyeyeva et al., 2019, Lindsey et al., 2014] as

$$ (a_1)_n = \cos(\theta + \beta_n) \sin(\alpha_n)v_H + \cos(\alpha_n)v_V, \quad (2) $$

where $\alpha_n$ and $\beta_n$ represent the varying incidence and azimuth angles of the LOS vector, respectively, and $\theta$ the local fault azimuth. We solve this equation minimizing the least-squares criterion. The azimuth of the CF varies between about 0° at 28°N latitude to 40° near Ghazni with significant local variations (Section 2.1 and Figure S5). To determine local fault azimuth, we resample the fault trace with regular spacing (100 m) and smooth its azimuth value with a Gaussian filter of 10 km length.

### 3.4 Measuring fault slip

We quantify the amount of slip along the CF, and evaluate the corresponding along-strike and temporal variations. In order to obtain this space and time vision, we work with both the time series of phase-change in LOS and the fault-parallel velocities. Effectively, we extract 500 m-wide fault-perpendicular profiles every 200 m along the fault [Kaneko et al., 2013, Jolivet et al., 2013].

In a first step, we adjust a screw dislocation model in an elastic half-space [Savage & Burford, 1973] to the 30-km-long profiles in fault-parallel velocities, $v_H$. We combine a shallow dislocation extending from the surface to a given depth, $D_C$, with a semi-infinite dislocation below a depth, $D_S$. The model is an approximation of the surface deformation produced by a vertical strike-slip fault creeping near the surface, above the depth $D_C$, at a constant rate $C$, and slipping at a constant rate $S$ below the locking depth, $D_S$. It writes as

$$ v_H(x) = -\frac{C}{\pi} \arctan\left(\frac{D_C}{x - X_f}\right) - \frac{S}{\pi} \arctan\left(\frac{x - X_f}{D_S}\right) + Bx + A \quad (3) $$

where $x$ is the distance along the profile, $A$ a constant, $B$ a ramp and $X_f$ the fault location [Segall, 2010]. We explore the values of $S$, $D_S$, $C$, $D_C$, $X_f$, $A$ and $B$ through Bayesian sampling of the parameter space using a Markov-Chain-Monte-Carlo algorithm [Salvatier et al., 2016]. The $A$ and $B$ terms account for the reference of our velocity map. The $a$ priori probability density function of $S$ and $C$ is uniform between 0 and 30 mm/yr (left-lateral slip), while the $a$ priori of $D_C$, $D_S$ and $X_f$ are bounded Gaussians. The fault location is allowed to vary by ±1.2 km from our mapped continuous fault trace. The creep extent, $D_C$ is within 0 and 8 km with an $a$ priori of 1±3 km and it has to be inferior to the locking depth, $D_S$, which is within 0.01 and 9 km with an $a$ priori of 2±5 km (see Text S2 and Table S2 for more details about the parametrization). We do not allow for greater locking depths as our profiles are short in order to avoid interactions with the nearby Rigestan desert to the west (Figure 2). Displacement features of wavelength larger than 10 km will be approximated by the additional linear term. As a result, slip below 8 km is not well constrained and we arbitrarily increase the associated uncertainty by 4 mm/yr.
We also aim to measure shallow fault slip without any underlying model. We consider the phase or velocity difference on each side of the fault as a direct measure of the surface expression of slip along the fault. We measure the phase step across the fault by subtracting the mean values within 500 m and 1.5 km of the fault trace on each side (e.g. profiles in Figure 3). We apply this procedure to the phase-change, $\phi_k$, at each time step, to the velocity field ($a_1$ or $v_H$), and to the associated velocity standard deviations (associated measures are referred to as $d\phi_k da_1$, and $dv_H$, respectively). We carefully computed surface slip rate uncertainties accounting for both the propagation of the previously estimated standard deviations on each velocity estimate and the scattering of velocity values spatially (Text S3). By considering points within 1.5 km of the fault trace, we capture slip at or close to the surface (down to 1.5 km depth) as predicted by Equation 3. Measures of surface slip based on points further away from the fault trace (e.g. within 1 and 5 km) show very similar patterns with larger slip amplitudes, but time series estimates are more affected by spatially correlated noise.

## 3.5 Slip inversion for 3 earthquakes

The three most recent earthquakes ($M_w$ 5.6, $M_w$ 5.1 and $m_b$ 4.1) in our time series occurred close to the CF, near the city of Chaman and were imaged by both ascending and descending InSAR data. Corresponding surface displacement fields (Section 3.2) are sharply contrasted on each side of the CF trace demonstrating that these events ruptured the CF itself close to or reaching the surface (Figure 2). We use these earthquake displacement maps to infer slip at depth along the CF using a constrained least squares inversion [CSI library; Elliott et al., 2016]. We build a three-dimensional fault plane with triangular elements following the CF trace between latitudes 30.3°N and 31.1°N, dipping vertically down to 10 km (Figure S6). We assume pure sinistral displacement on a vertical fault, in agreement with focal mechanisms and with Barnhart [2017]. We downsample the reconstructed InSAR map of earthquake displacement ($a_{6,7,8}$ in Equation 1) and associated uncertainties using a quadtree algorithm based on model resolution [Lohman & Simons, 2005], with element size between 600 m and 10 km (Figure S7). We compute the Green's functions relating slip on the fault to surface observations for point sources in a homogeneous elastic half-space with a Poisson ratio of 0.25 [Zhu & Rivera, 2002]. We use data uncertainties obtained from the output of KFTS to build the diagonal data covariance matrix, $C_d$. The *a priori* model Covariance, $C_m$, is based on a decreasing exponential of the distance between fault elements [Radiguet et al., 2011]. We adjust defining parameters of $C_m$, that is, the amplitude of correlation ($\sigma_m$) and the characteristic length scale ($\lambda$), while we fix the normalizing distance to 1 km (see L-curve in Figure S8). We find that the combination of $\sigma_m=3$ mm and $\lambda=2$ km results in a model that is both close to the data (i.e., small misfit) and physically sound and smooth (i.e., no large values of parameters for negligible drop in misfit). Additionally, we vary the value of $\sigma_m$ from this reference value to regularize the inversion behavior along fault segments not covered by data. We add 1 mm to $\sigma_m$ for fault elements below the spatial footprint of data down to the bottom of the fault and remove 1 mm outside (Figure S9).

# 4 Results: InSAR-derived deformation

## 4.1 Ground velocity

The InSAR-derived velocity fields display signals with local and regional length scales related to tectonics, hydrology and human activity (Figure 2). Earthquake related deformation is treated separately thanks to our parametric decomposition of the signal (Section 3.2). The standard deviation on ground velocity in LOS is close to 1 mm/yr for most pixels, except on the edge of the Rigestan Desert, in the Indus plain, above lakes and around the high peaks of the Hindu Kush (north of our study area), where many interferograms could not be unwrapped due to low coherence of the signal. Accounting for
earthquakes in our parametric model of deformation increases uncertainties in local ve-
locity estimates to ~2 mm/yr due to existing trade-offs between adjusting a ramp and
a step function in the time series of interferometric phase change (Figure S10).

Both velocity fields along ascending tracks show a ~100 km wide gradient of deformation perpendicular to the fault of about 20 mm/yr in the LOS direction. This gradient seems larger along two roughly north-northeast-striking regions : the Chaman fault (CF) and the central Kirthar range, west of the CF. We interpret this signal as the accumulation of interseismic strain across a transform plate boundary. In this study, we focus on strain accommodation by the CF only. In addition to the long-wavelength signal (~10 km), sharp discontinuities across the CF are visible. In the fault-parallel projection of velocity (Figure 3), such contrast is of the order of 1-2 cm/yr, with significant along-strike variations in surface slip rate. Dislocation model adjustment along 30 km-long selected profiles in Figure 3 confirms this near-field strain accumulation rate (S from 9 mm/yr in profile C to 20 mm/yr in profile A) and suggest very shallow locking depths along the fault plane (D<sub>S</sub> from 1 km in profile D to 5 km in profile A).

Strong local minima in velocity in both ascending and descending tracks are mostly subsidence from anthropogenic water pumping and, to a lesser extent, mining. They are clearly retrieved in the vertical decomposition of the signal with rates close to 15 mm/yr in Kandahar and near Qalat and up to 200 mm/yr in Quetta and Pishin basins, a region known as the 'fruit orchard' of Pakistan (Figures S13 and S14). Other studies retrieved comparable subsidence rates [Szeliga et al., 2012, Kakar et al., 2016, Ahmad et al., 2017].

## 4.2 Slip along the Chaman fault

Figure 4 displays the spatio-temporal distribution of surface slip in the LOS, $da_1$, along the CF with slip accumulating over time. The relative displacement caused by the three earthquakes between 30° and 31°N clearly stands out. Interestingly, these events occurred on fault portions that appear to be slipping beforehand at rates between 1 and 3 mm/yr. It is most obvious for the 2018 event, where ~10 mm of LOS surface slip accumulated in the 3.5 years before the earthquake. We refer to this fault portion which creeps and hosted three earthquakes between roughly 30.4°N and 31°N as the Central segment of the CF. Moreover, we see pronounced surface slip on each side of the Central segment. In the south, we recognize the Nushki creeping segment [Barnhart, 2017] between about 0 and 80 km north of the city of Nushki, Pakistan, with a surface slip rate of about 1.4±0.5 mm/yr in LOS. Another portion of the fault between 200 and 350 km north of Nushki slips with a maximum rate of 3.4±0.8 mm/yr in LOS (Figure 5b). We refer to this fault portion as the Qalat segment. On both portions, the temporal evolution of slip is linear and continuous at the temporal scale of SAR acquisitions (several days). We do not identify any obvious transient slip events apart from post-seismic transients (see Section 5.2).

Left-lateral slip rate close to the surface in Figure 5c is quantified either measuring the
offsets in velocity on each side of the fault, $dv_H$, or using the median of the modeled creep
rate $C$ (Equation 3). Slip rate at depth (or loading rate) in Figure 5d is modeled slip
rate, $S$, below the locking depth $D_S$. These measures are inferred from the horizontal,
fault-parallel velocity ($v_H$ in Equation 2). In the following, we detail along-strike vari-
ations in slip rates using the distance from the city of Nushki as reference. A strong neg-
ative slip rate (i.e. apparent right-lateral slip) estimate around 180 km comes from a 10-
km-wide region of low velocity anomaly west of the fault. This region is collocated with
the sharp releasing bend of the CF (Figure 5a) suggesting local extension and subsidence.
Additional subsidence may originate from irrigation and human settlements clustering
around the Dori river (Figure S18). Therefore, the validity of the assumption of fault-
parallel horizontal motion is questionable at this location (170-185 km) and slip estimates
cannot be trusted.
Highest surface left-lateral (i.e. positive) slip rates ($dv_H \sim 8$ mm/yr) are measured south of 180 km, on three fault segments separated by null surface slip rate within uncertainty: the 80-km-long Nushki creeping segment, a 30-km-long segment around 120 km and a 20-km-long segment around 165 km. Measures of surface slip rates north of 180 km are lower but still positive with a notable peak at $5.6 \pm 1.7$ mm/yr between 275 and 300 km. Variations of $C$ and $dv_H$ are close, however, $C$ tends to reach lower slip rates because of the larger uncertainty arising from the use of a model [e.g. Duputel et al., 2014] and because of the existing trade-off between $C$ and $S$ when $D_C$ and $D_S$ estimates are close to each other, such as along the Nushki segment (Figures 5d and S23). North of 380 km, measured surface slip rate is zero $\pm 1.7$ mm/yr, an uncertainty arising from high near-fault variability in velocity measures. South of Nushki, low and slightly positive $dv_H$ is recorded with an uncertainty of 0.5 mm/yr while the mean interquartile range of $C$ is 0.5-2 mm/yr; indicating no certain surface slip.

Previously identified creeping segments have locking depths, $D_S$, of $4 \pm 2$ km and a locked depth range, $D_S - D_C$ of $2 \pm 2$ km suggesting that the whole width of the fault plane slips with a change in rate at about $D_S$. The 20-km-long portion between the Nushki and Central segments displays a $D_S - D_C$ of $3.5 \pm 1.3$ km while in between 190 km and the Qalat segments $D_S - D_C$ is $4 \pm 2$ km. Moreover, positive locked depth range is also modeled around 140 km ($3.4 \pm 1.1$ km, $D_S \sim 4$ km), within the Central segment, and between 310 and 325 km ($5 \pm 2$ km) in what we named the Qalat segment (Figures 5d and S21). Modeled slip rates at depth, $S$, are consistently non-zero along the fault portion between -140 and 380 km with limited along-fault variations. Between -80 and 165 km $S$ are $12 \pm 3$ mm/yr, while it is $7 \pm 2$ mm/yr between 190 and 380 km. Therefore, the change in fault azimuth around 180 km appears as a discontinuity in the CF loading rate. North of 380 km, the median model is 1 mm/yr of slip (0.6-2 mm/yr interquartile range) below $5.6 \pm 2$ km depth with large uncertainties, suggesting that there is no significant loading of the fault at high latitude, in agreement with the absence of across fault gradient in velocity maps (Figures 2a,b and 3). Lateral variability in modeled slip north of 380 km probably arises from poor fault location and numerous basins along the hypothetically continuous fault trace. South of the city of Nushki, we obtain significant amount of slip at depth, especially between -130 and -90 km where it is in average $21 \pm 4$ mm/yr below $6 \pm 2$ km (e.g. Profile A in Figure 3). This corresponds to a strong across fault gradient in projected velocities (Figure 3). However, the near north-south orientation of the fault trace on this portion implies that uncertainties associated with fault-parallel velocities are larger than 5 mm/yr (Figure S15). Moreover, modeled slip is even more uncertain because of the elevated locking depth associated. Therefore, we conclude that the fault is locked at the surface south of Nushki down to at least 5 km, but the amount of strain accommodated is uncertain (probably 5-20 mm/yr).

As a summary, we delimit the $\sim 400$ km-long shallow slipping CF between the latitude of the city of Nushki, 29.5°N, and about 32.6°N ($\sim 380$ km). This length includes creeping segments that slip continuously at 4-9 mm/yr over 2015-2020 from the surface to about $4 \pm 2$ km, a depth below which higher slip rates are expected (roughly 7-12 mm/yr). Creeping segments are separated by 5 40-km-long segments locked down to at least 3-6 km depth. The releasing bend of the fault 180 km north of Nushki ($\sim 31.1°$N) marks the frontier between two fault portions of similar length, the southern portion slipping faster than the northern portion by about 3 mm/yr. North of 32.6°N, the CF seems inactive at our scale of observation, whereas the fault appears to slip below roughly 6 km south of the city of Nushki. In the following section, we further explore the slip behavior of the central CF between the creeping Nushki and Qalat segments where earthquakes have been identified.
<figure><img src="image_2.png" /><figcaption>**Figure 2.** Velocity field for ascending tracks (in (a)) and the descending track (in (b)) in the line of sight (LOS) direction with histograms of their associated standard deviations (std). The histogram's vertical axis is the number of pixels in a logarithmic scale (full std map is in Figure S10). Black line is the continuous trace of the Chaman fault used as reference. Filled black square markers locate cities named in Figure 1. Black frames delimit area in subplots (c)-(g) of the estimated coseismic displacement according to the ascending track 42 (see Figure S11 for descending track and Figure S12 for coseismic interferograms).</figcaption></figure>
<figure><img src="image_3.png" /><figcaption>**Figure 3.** Velocity field projected in fault parallel direction (Figure S5) with histogram of standard deviations (std) and 4 across fault profiles (see Figure S15 for the full std map and Figure S16 for a fixed direction of projection). Data with std >8 mm/yr is not displayed. Note that to the south (left of this plot) the fault azimuth is close to zero implying that projected InSAR velocities are more uncertain (± 6 mm/yr). Velocities where tectonic deformation is known to be secondary compared to anthropogenic activity or geomorphological processes are hatched. Profiles are taken at distances along the fault from Nushki of -120 (A), 40 (B), 100 (C) and 300 (D) km. Measured surface slip from points close to the fault (0.5-1.5 km in black) are outlined by double headed arrows (those values are compiled along the whole fault in Figure 5c). Adjusted dislocation model for a vertical left-lateral fault is shown by the red curve with the median estimate and interquartile range of parameters in Equation 3 written in red.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 4.** Temporal evolution of surface slip along the Chaman fault in line of sight (LOS) for ascending tracks. The color gradient represents cumulative surface slip of seismic and aseismic origin with respect to the first date (see Section 3.4 for method). Along most of the fault, positive slip indicates a left-lateral sense of slip (depending on the angle between fault azimuth and LOS). The timing and location of the 3 identified earthquakes are marked by black stars. Named portion of the fault with their approximate extent are indicated by labeled dashed arrows. The spatial reference (i.e., zero of the vertical axis on the left) is the latitude of the city of Nushki (29.556°N), Pakistan. For a variant see Figure S17.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 5.** Geometry, measured slip rates and seismicity along the Chaman fault. (a) Fault trace and its first derivative, quantifying the variation in fault azimuth along its length in radian. (b) Surface slip rate according to the across fault velocity offset, *da*₁, in the direction of the line of sight (LOS) for our two ascending tracks (T42 and T144) as well as from Fattahi & Amelung [2016] study (F&A 2016). Notice that the LOS varies from pixel to pixel for Sentinel, and that F&A 2016 data are acquired with Envisat in another (but close) geometry. (c) Surface slip rate projected in fault-parallel direction (Section 3.3) according to the velocity offset on each side of the fault, *dv*<sub>H</sub>, or the median of the modeled creep rate, *C* (Section 3.4). (d) Mean of the modeled slip rate *S* below the locking depth, *D*<sub>S</sub>, colored as a function of the locked depth range (difference of the median depths, *D*<sub>S</sub> − *D*<sub>C</sub>). While shading of *da*₁, *dv*<sub>H</sub> and *S* are the one-standard deviation range, shading around *C* is the interquartile range reflecting the asymmetry of estimated creep rate distributions. Slip rates in (b), (c) and (d) are low-pass filtered, to remove frequency higher than 1/10 km<sup>-1</sup> (Figures S19-S21). (e) Histogram of seismic moment releases from 1900 to 2019 in 10 km wide bins along the CF (Section 2.2). Bar colors indicate the number of events summed per bin. Equivalent M<sub>w</sub> are on the right axis (Figure S22 for other seismic catalogs). Our understanding of the fault segmentation based on aseismic slip rate lateral variations is sketched in (f).</figcaption></figure>
## 4.3 Seismic events

### 4.3.1 Surface displacement fields

Figure 2c-g displays earthquake related displacement fields imaged by ascending track 42 (see Figure S11 for descending track). Our time dependent model simplifies this displacement as a Heaviside function of time without specifically accounting for postseismic deformation but allows to extract surface deformation for quite small earthquakes [Liu et al., 2021]. Postseismic slip rate is typically the largest immediately after the earthquake and then decays rapidly with time [e.g. Perfettini & Avouac, 2004]. As a result, most of the postseismic deformation is included in the computed amplitude of the Heaviside function of time, implying that coseismic and postseismic slip are not distinguishable using our parametric description of interferometric phase change alone.

The earliest (August 2015) earthquake of our time series is a M<sub>w</sub>5.4 thrust event on the Ghazaband fault, with maximum uplift close to 8.5 cm in LOS and a rupture length of about 10 km. According to GCMT, it was preceded by a M<sub>w</sub>5.3 earthquake 19 days before and followed by a m<sub>b</sub>4-4.3 two days latter at distances within expected location uncertainties (see Figure S1 for a compilation of magnitude estimates). Surface motions associated with those smaller earthquakes may be included in our InSAR-derived deformation. These earthquakes are close to the northern termination of the 2013 Balochistan earthquake sequence and may relate to the induced stress change [e.g. Jolivet, Duputel, et al., 2014, Avouac et al., 2014]. All other studied earthquakes are essentially left-lateral strike-slip events. The M<sub>w</sub>5.7 earthquake that occurred on 21 March 2016 (Figure 2d) is located on the northern termination of the Ornach Nal fault and displays up to 4 cm relative displacement in LOS along a 20-km-long fault segment. It is in the spatio-temporal vicinity of nine other earthquakes of smaller magnitudes for a total moment equivalent to a single M<sub>w</sub>5.73. Because the August 2015 and March 2016 sequences are off the CF, we will not study them further.

Our images of surface displacement produced by the last three earthquakes clearly show that they occurred on the CF itself. The horizontal displacement fields from the combination of the two LOS (Equation 2) in Figure 6 are essentially made up of two lobes along the fault trace, a positive one west of the fault and a negative one to the east indicating left-lateral motion. The symmetry and shape of the deformation lobes are close to the expected shape of the quadrants of compression for an idealized strike-slip earthquake in a homogeneous medium. Deformation within the extensive quadrants is not retrieved due to the geometry of the ascending and descending LOS with respect to the fault orientation. The May 2016 (M<sub>w</sub>5.6) earthquake resulted in a maximum of 8 cm of relative horizontal displacement at the surface and is located right at the junction of a thrust structure to the west (i.e. the Spinatizha fault) and the CF main strand. Two months later, a smaller (M<sub>w</sub>5.1) earthquake occurred north of it with similar amplitude of displacement on a shorter fault segment. Note that ~5 mm of what looks like postseismic slip of the May 2016 event is included in the displacement field of the July 2016 event (Figure 6) due to their spatio-temporal proximity. In June 2018, a strike slip event produced about 2 cm of surface displacement along a 10 km long fault portion.

While the M<sub>w</sub>5.1 July 2016 earthquake is isolated, the other two individualized events seem to correspond to groups of earthquakes. The 13 May 2016, three earthquakes in three minutes occurred with m<sub>b</sub>4.9-5.3, m<sub>b</sub>4.7-4.8 and M<sub>w</sub>5.6 for a total seismic moment equivalent to about M<sub>w</sub>5.7 (Figure S1). The 27 June 2018, a second earthquake occurred two hours after the m<sub>b</sub> 3.8-4.2, but with m<sub>b</sub> 3.6 so that it appears negligible with respect to the main earthquake. Therefore, out of the three events on the CF, the May 2016 event only is treated as a sequence of earthquakes (with total M<sub>w</sub>5.7). We now examine the results of the source inversion for these three events.
<figure><img src="image_6.png" /><figcaption>**Figure 6.** Surface horizontal displacement induced by each of the three earthquakes close to the town of Chaman during their coseismic and postseismic periods combining ascending and descending InSAR observations (Equation 2). Fault azimuths are assumed constant and equal to 26.7°, 17° and 29° for each event respectively, from left to right. Interpreted direction of horizontal motion for the May 2016 earthquake are shown by white arrows and expected extensive quadrants are circled with white dashed lines. The geometry of the InSAR line of sights (LOS), imply a poor sensitivity to the motion in those extensive quadrants, likely to be as large as recorded motion in the compressive quadrants. This is also true for the July 2016 and June 2018 events. The epicenter of large historical earthquakes attributed to the Chaman fault are located by purple circles with the associated range of moment magnitude (M<sub>w</sub>) estimates and year of occurrence labeled in the left plot [Ambraseys & Bilham, 2003a]. Spina Tiza is a frontier post of Pakistan.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 7.** Modeled earthquakes slip on the fault plane and associated surface displacement fields for the earthquakes that occurred on 13 May 2016 (top), 10 July 2016 (middle) and 27 June 2018 (bottom). **Left:** Inferred moment magnitudes (M<sub>w</sub>) from modeled slip distributions are specified. Red arrows show the fault extent covered by the available deformation field for each earthquake and, thus, under which part of the fault, slip is constrained. The vertical exaggeration is 2. **Right:** Synthetics horizontal surface deformation fields comparable to observations in Figure 6 and residuals. The black line is the continuous Chaman fault trace used for inversion.</figcaption></figure>
4.3.2 Results of the slip inversion

Inverted slip at depth reaches the surface for all events with a maximum slip of ~10, 8 and 4 cm for the May 2016, July 2016 and June 2018 earthquakes, respectively (Figure 7). We observe a chronological decrease in area and amount of slip, consistent with the estimated magnitude from seismic station records (Table 1). The moment magnitude (M<sub>w</sub>) of the May 2016 and July 2016 events are 5.7 and 5.1, respectively. For the 2018 event, the National Earthquake Information Center (NEIC) estimates a m<sub>b</sub> of 4.1, equivalent to M<sub>w</sub> = 4.4±0.2 [Scordilis, 2006]. Our modeled slip on the fault, converted to M<sub>w</sub>, yields moment magnitudes of 5.6, 5.3 and 5.0 for the three earthquakes in chronological order, respectively.

The source of the May 2016 earthquake is a roughly circular patch centered at about 5 km depth. To first-order, slip of the July 2016 is also best represented by a roughly circular patch at shallow depth (0-5 km depth) where most of the surface deformation is imaged (around 30.65°N). 20 km to the south, a subtle slip signal is interpreted as post-seismic slip remaining from the May 2016 event, while no slip is observed in between above 5 km depth (Figures 6 and 7). This presumably postseismic signal reaches amplitudes of 1 cm on the fault plane and affects a large area. To improve our estimated M<sub>w</sub> of the May and July 2016 earthquakes, we isolate the sub-source moment induced by fault-slip between 30.6 and 30.8°N (dashed lines in Figure 7) and find a corresponding M<sub>w</sub> of 5.0. Re-attributing the postseismic moment to the May 2016 event would not change our M<sub>w</sub> estimate of 5.6, whereas the July 2016 M<sub>w</sub>, ignoring slip between 30.6 and 30.8°N, becomes 5.2.

The residuals between observed and modeled surface displacement are one order of mag-
nitude smaller than the modeled earthquake-related displacement for the July 2016 and
2018 events (Figures 7 and S24). Regarding the May 2016 event, residuals are large (+2
cm in horizontal) with a spatial distribution indicating (i) unmodeled deformation in the
south-west lobe and (ii) the difficulty to model the fact that the north-east lobe is off-
set to the east with respect to the fault trace. This complex spatial pattern may be the
consequence of a multi-fault rupture by the three earthquakes involved. Therefore, we
only model a fraction of the energy released during the coseismic and postseismic phase,
implying M<sub>w</sub> estimates from modeled slip of the May 2016 may be underestimated.
# 5 Discussion

## 5.1 A decadal perspective on Chaman fault slip

Fault slip along the 700-km-long Chaman fault (CF) is shown in Figure 5 together with a description of the continuous fault trace and the compilation of the seismicity recorded since the beginning of the 20th century. In our interpretation, we divide the CF into creeping, locked and uncertain (probably inactive) segments combining the vision given by fault parallel surface slip rates, modeled slip rates and locking depths only (Figure 5f). Because we do not *a priori* consider the fault trace geometry or the seismicity, our segmentation is purely based on the description of aseismic slip and does not pretend to provide units that may rupture in a single earthquake. In this section, we compare our slip rates and segmentation with the seismic record and past slip estimates. We identify two 80-km-long continuously slipping segments: the Nushki segment and the southern part of the Qalat segment (Q1), with mean surface slip rate of 6±1 and 4±2 mm/yr, respectively, measured within 3 km of the fault trace (Figure 5c). Three shorter creeping segments are also outlined, two in the Central segment with slip rates comparable to the Nushki segment and separated by a ~5-km-long locked portion, which significance is unclear, and one 50-km-long segment in the Qalat segment (Q2), with less clear evidence: a locking depth of 3.5±0.5 km and a slip rate at the surface of 2±0.8 mm/yr. These creeping segments are separated by ‘locked’ segments with null or uncertain near-fault slip rate, a locking depth exceeding 4 km and a non-zero modeled slip rate at depth. The fault portions north of the Qalat segment as well as between 165 and 190 km are labeled ‘uncertain’ as we cannot identify left-lateral strain on the fault with our data.

Our observed distribution of surface slip along the CF shows good agreement with previous studies and independent data sets. With the method described in Section 3.4, we extract slip rate estimates along the CF between 2004 and 2011 from Envisat velocities [Fattahi & Amelung, 2016]. The mean slip rates in 2004-2011 and the previously described slip rates over 2014-2019 exhibit very similar along-strike variations as shown in Figure 5b. Additionally, our conversion from LOS to left-lateral slip rate using both viewing direction of Sentinel 1 satellites matches values derived by Barnhart [2017] from previous SAR missions (Envisat and ALOS) with peak creep rate around 10 mm/yr (Figure 5c). Barnhart [2017]’s locking depth routinely shallower than 500 m along the Nushki segment cannot be directly compared to our value as the fault model is different. Therefore, we have robust observations representative of the aseismic surface slip rate along the CF which seems constant over at least the past 15 years, between ~28°N and ~31°N. In general, shallow aseismic slip (*C* or *dv*<sub>H</sub>) is significantly lower than the apparent loading rate (or slip rate, *S*) suggesting stress increases despite the occurrence of aseismic slip.

We compare this seemingly continuous aseismic slip rate with local slip from earthquakes. We consider the past-century seismic record (1900-2019) and assume that earthquakes within 30 km of the CF occurred on the fault itself (more details in Section 2.2). We observe a total seismic moment release of 1.8×10<sup>19</sup> N.m along the ~700 km of fault in 120 years (Figure 5e, Table S3). The frequency-magnitude relation [Gutenberg & Richter, 1944] indicates that the record of earthquakes with M<sub>w</sub><4.5 on the CF is incomplete, and leads to an estimated *b*-value of one (Figure S25). On this basis, the cumulative moment released associated with unrecorded M<sub>w</sub>0.1 to 4.5 over the same period of time is equal to 7.3×10<sup>18</sup> N.m. The combination of recorded seismicity and extrapolated microseismicity reaches then 2.5×10<sup>19</sup> N.m, a value equivalent to the moment released by a M<sub>w</sub>6.8 earthquake. Based on this estimate, assuming a shear modulus of 30 GPa and a locking depth of 3-6 km, the seismic moment released over 100 years corresponds to slip rates from 3.5 to 1.5 mm/yr, respectively. This suggests that seismicity along the CF contributes to less than 15% of the relative plate motion (about 30 mm/yr) over the last century. Over this time period, the recorded M<sub>w</sub>5.6 in May 2016 appears to be exceptional, as were the nearby M<sub>w</sub>~6.6 of 1975 and 1892 [Bilham et al., 2019]. Nonethe-
less, in the region, large historical earthquakes ($M_w \ge 7$) were recorded east of the CF within the fold and fault belt (e.g. 1931 $M_w$7.3 Mach and 1935 $M_w$7.7 Quetta earthquakes) [Ambraseys & Bilham, 2003a, Dewey et al., 2006]. Hence, assessing the partitioning of deformation between the CF and structures in the fold and fault belt is necessary to assess the seismic hazard of the region, and the CF specifically, but this is out of the scope of this study.

Although the total moment released by earthquakes is small, the along-strike distribution of seismicity and moment is key to assess the relationship between earthquakes and continuous aseismic slip. Historical records attest that earthquakes occurred on both the Nushki and Qalat creeping segments (Figure 5e). A seismic crisis between 1975 and 1978 struck the fault 30-100 km north of the city of Nushki with at least four $M_w > 5.5$ earthquakes [Lawrence & Yeats, 1979, Ambraseys & Bilham, 2003a]. Nonetheless, the largest event ($M_w$6.4-6.8 Spina Teza, also spelled Spin Tezha, earthquake) dating back to October 1975 falls on the locked segment at the northern limit of the Nushki segment with its epicenter at about 85 km. This segment hosted the epicenter of about 15 events since 1900 and, according to our model, accumulates elastic strain. The adjacent creeping segment to the north appears to be particularly seismically-active too and hosted the largest aftershock of the Spina Teza earthquake (with estimated $M_w$ between 5.8 and 6.5), the May 2016 $M_w$5.6 and the June 2018 event. Those two segments focusing most of the recorded seismicity are along the ~100-km-long restraining bend in the central CF (Figure 5a,e). Moreover, a nearby creepmeter (157 km north of Nushki) recorded a slow slip event starting on the 1st of March 2019, and accommodating 4.5 mm of left-lateral slip in 40 days [USGS; Bilham et al., 2019]. This suggest that what we imaged as continuous aseismic slip may include discrete transient slip accelerations. Further north, the only individual earthquake studied, a $M_w$5 in 2005 relocated with InSAR [Furuya & Satyabala, 2008, Fattahi & Amelung, 2016], occurred on the 20 km-long locked segment between Q1 and Q2. Considering that our compilation of seismic event generously includes events within 30 km of the fault trace, we have no direct evidence of any earthquake on the Qalat creeping segments itself. The peak of seismicity more than 420 km north of Nushki (~33°N) is on a portion of the CF on which we are not able to identify tectonic strain accumulation (labeled as ‘uncertain’). Looking at earthquakes in map view (Figure 1), the corresponding events are located around subsiding area southeast of the CF (south of the city of Ghazni and along the Gardez fault zone) with negative signal in LOS and deemed fault parallel motion (Figures 2a,b and 3).

To summarize, we measure surface slip and its lateral variations due to aseismic processes stable on a decadal scale, while recorded earthquakes have induced very limited displacement over the past century. We divide the fault into 5 creeping segments separated by locked segments according to the measured slip and locking depth. We find that the central restraining bend is the most seismically active section of the fault and hosts in close spatio-temporal relationship seismic and aseismic slip. In the following, we investigate the importance of postseismic signal for the three recent earthquakes imaged for the accommodation of tectonic stress before discussing mechanical implications.

## 5.2 Postseismic signal characteristics in comparison with coseismic deformation

Two studied $M_w$5 (21 October 2005) and 5.5 (19 October 2007) earthquakes within the Chaman fault zone display abnormally large and long postseismic signal with respect to what is usually measured for large ($M_w$>6) earthquakes [Furuya & Satyabala, 2008, Fattahi et al., 2015, Alwahedi & Hawthorne, 2019]. Postseismic slip lasted for more than a year with postseismic to coseismic moment ratio estimates of 0.7±0.1 [Fattahi et al., 2015], and 1.1 [Furuya & Satyabala, 2008], while it is typically less than 0.3 for large earthquakes (in California) [Alwahedi & Hawthorne, 2019]. We further analyze the spatio-temporal
<table><caption>**Table 1.** Properties of the three earthquakes imaged on the Chaman fault. InSAR M<sub>w</sub> are computed from modeled slip at depth. * This magnitude is equivalent to the sum of the seismic moment from the three earthquakes included in the event. ** This magnitude is a conversion using Scordilis [2006]'s relationship from a body-wave magnitude (m<sub>b</sub>). † Ascending line of sight (LOS) slip rate before the earthquake occurrence looking at a distance to the fault of 0.1-1 km. ‡ This is a lower bound estimate.</caption><thead><tr><th>Date</th><th>Latitude</th><th>LOS slip rate<sup>†</sup><br>(mm/yr)</th><th>seismic<br>M<sub>w</sub></th><th>InSAR<br>M<sub>w</sub></th><th>post/co-seismic<br>moment</th><th>post/co-seismic<br>slip<sup>‡</sup></th></tr></thead><tbody><tr><td>13 May 2016</td><td>30.63</td><td>1.2 ±0.5</td><td>~5.7*</td><td>≥5.6</td><td>>0.09</td><td>1.3</td></tr><tr><td>10 July 2016</td><td>30.78</td><td>1.1 ±0.5</td><td>5.1</td><td>5.2</td><td>0.4</td><td>0.4</td></tr><tr><td>27 June 2018</td><td>30.50</td><td>2.4 ±0.1</td><td>4.4 ±0.2**</td><td>5.0</td><td>3-15</td><td>0.8</td></tr></tbody></table>

footprint of our three earthquakes on the CF to explore potential mechanical characteristics of the fault zone.

We first look back at the moment magnitude (M<sub>w</sub>) estimates from the modeled slip on the fault plane (Section 4.3.2). Because the reconstructed surface displacement includes coseismic and postseismic slip, we can assess postseismic moment by subtracting the moment computed from seismic waves to our modeled moment. This postseismic moment would also include preseismic transient slip, if any. For the May 2016 event, the estimated M<sub>w</sub> of 5.6 is smaller than the total seismic M<sub>w</sub> equal to about 5.7. Nonetheless, residuals indicate that our model does not explain a large part of the observed deformation signal, suggesting that the M<sub>w</sub> from the model underestimates the InSAR M<sub>w</sub> (Section 4.3.2). Moreover, we identify a fraction of the May 2016 postseismic slip equivalent to a M<sub>w</sub>5 and, thus, the postseismic to coseismic moment ratio is at least 0.09. Regarding the July 2016 event, we estimate a M<sub>w</sub>5.2, while seismological catalogs indicate a seismic M<sub>w</sub>5.1, which converts into a postseismic to coseismic moment ratio of 0.4. The last recorded event of June 2018 is associated with a seismic M<sub>w</sub> of 4.4±0.2 much smaller than the M<sub>w</sub>5 we estimate, leading to a postseismic to coseismic moment ratio between 3 and 15 (Table 1). Therefore, the 2018 event appears mostly as an aseismic event, an idea further supported by the time series of surface slip at this location (Figure 8e), in which the step across the fault results from a month-long transient event rather than a clear cut in between two acquisitions.

We now consider the relationship between coseismic and postseismic deformation using time-series in the line of sight (LOS). From InSAR time-series, the closest measure of coseismic displacement is given by the relative phase change between 6 days framing the time of the earthquake and, thus, it also includes a few days of postseismic slip. We divide the subsequent fault slip into what we name an early postseismic phase lasting 2-3 months and a long-term postseismic in the 300-500 days after the early postseismic phase (pink and blue shadings in Figure 8a,c,e). Corresponding across-fault profiles for each period are in Figure 8b,d,f. Time series of fault displacement on each side of the fault, averaged between distances of 0.1-1 km of the fault trace, display characteristics logarithmic decay of slip in the months to years following the earthquakes which we interpret as afterslip (Figure 8a,c,e; method in Section 3.4) [Perfettini & Avouac, 2004, Thomas et al., 2017].

Profiles in Figure 8 allow to compare the amount of slip and its spatial wavelength for the coseismic, early and postseismic phase. The spatial wavelength of deformation qual-
<figure><img src="image_8.png" /></figure>

**Figure 8.** Temporal and spatial footprint of coseismic and postseismic deformation according to track 42 for earthquakes on the Chaman fault. **Left:** time series of differential displacement across the fault in LOS at the earthquake location on a 1-2 km long fault segment. The distance in the upper left corner refers to the distance along the fault to Nushki (Figure 4). Displayed velocities refer to the best-fitting slip rate before the earthquake. The differential parametrized model from KFTS on each side of the fault is the black line. The dark blue vertical dashed line indicates the timing of the earthquake and shaded regions correspond to the periods over which the displacement in the plots on the right are measured. **Right:** Profiles of displacement in LOS reflecting the spatial distribution of slip during the coseismic phase as well as during the immediate and longer postseismic phases. Coseismic deformation is taken as the displacement recorded between the two SAR acquisitions framing the earthquake occurrence (dark blue). The local extrema used to compute peak-to-peak displacement are marked by white symbols outlined in black (pentagones for minima, circles for maxima). Profiles cover contiguous but mutually exclusive periods.

itatively locates slip at depth, as predicted by Equation 3 (i.e. the wavelength of surface deformation roughly equals the depth of slip). First, the almost identical shape of the coseismic and early-postseismic profiles for the event in May 2016 (Figure 8b) argues towards an overlap of co- and postseismic slip with comparable amplitudes. Second, the steep across fault gradient of the July 2016 coseismic and early postseismic, indicates that fault slip reached the surface during the earthquake and during the early postseismic phase. However, postseismic signal is about three time smaller than co-seismic for this event. Third, coseismic displacement related to the 2018 event contrasts with the sharp postseismic signals lasting over a few months indicating the earthquake did not reach the surface and triggered postseismic afterslip near the surface. Taking the specific profiles in Figure 8b,d,f and their highlighted optima, we measure peak-to-peak postseismic to coseismic surface slip ratios considering together what we named ‘early’ and ‘long-term’ postseismic and subtracting the preseismic slip rate on our period of post-
seismic observations (Table S4 compiles all values). We find slip ratios of 1.3, 0.4 and 0.8 for the May 2016, July 2016 and 2018 events, respectively (Table 1). Because our InSAR coseismic images include the first few days of postseismic expected to have the highest slip rates according to logarithmic decay predictions, our postseismic to coseismic surface slip ratios are lower bounds.

Our data for the May 2016 sequence and June 2018 earthquake show that transient postseismic slip may accommodate as much slip as coseismic slip along the CF and, thus, is key in tectonic strain accommodation. In particular, the June 2018 event imaged with InSAR appears as predominantly aseismic with its postseismic moment being 3-15 times larger than the $M_w$ estimated from seismic waves. Our computation of a postseismic to ‘coseismic’ slip ratio of 0.8 only, indicates that our InSAR coseismic computed using an acquisition on the 7th of July, 10 days after the earthquake occurrence, includes most of the aseismic slip. The significance of the aseismic slip following the May 2016 event is supported by the postseismic to coseismic surface slip ratio of 1.3 for this earthquake. On the other hand, the July 2016 earthquake appears to scale like classical earthquakes, with post- to coseismic ratios in moment and slip equal to 0.4.

Finally, these three events are close to the rupture of the 1892 $M_{w}$6.5-6.7 Chaman earthquake [Bilham et al., 2019] and the 1975 Spina Teza earthquake sequence (with a main-shock of $M_{w}$6.4-6.8) [Ambraseys & Bilham, 2003a, *ISC On-line Bulletin*, 2020]. The extent of the 1975 Spina Teza surface rupture is not precisely known, but it likely terminated close to the southern termination of the 1892 rupture where the May 2016 earthquake locates (Figure 6). This is consistent with either a complex fault structure that would have stopped the rupture [Wesnousky, 2006] or with local stress concentration left in the wake of the largest events recognized on the CF. The spatio-temporal proximity of the July 2016 earthquake and its inverted slip pattern on the fault plane connected to the May 2016 event (Figure 7) indicate that this second event could have been triggered by the first one. Regarding the June 2018 event, the high postseismic to coseismic ratio could lead us to consider this transient deformation event either as a classic co- and postseismic sequence, or as a slow slip event triggered by a small earthquake (with $m_{b}$ estimates between 3.8 and 4.2). We therefore consider as a followup question whether other triggered or spontaneous events can be found along this fault segment, like the one recorded by a creepmeter in 2019 [Bilham et al., 2019].

## 5.3 Variations in slip behavior, implications for rheology and fault system orientation

We observe a spatial segmentation in slip behaviors along the fault with complex interactions between seismic and aseismic slip. What appears as continuous aseismic slip includes the collocated occurrence of earthquakes or slow slip events, and our identified long-lasting Nushki creeping segment has hosted earthquakes in the 1970s. Therefore, the traditional divide between seismically-active locked segments and creeping segments does not fully apply here.

Fault slip is the frictional response of a fault to stress loading. Classically, considering a rate-and-state formalism, aseismic slip is promoted by rate-strengthening materials, such as serpentine or various clay minerals [Dieterich, 1979, Marone, 1998]. Lawrence & Yeats [1979] mapped serpentinite bodies in the Chaman fault zone at least between 30.55°N and 30.9°N, which hints at velocity strengthening regime along the fault section that, paradoxically, hosted the May and July 2016 earthquakes as well as the $M_{w}$6.5-6.7 1892 earthquake. Hence, the fault plane along this section must also have rate-weakening patches to allow earthquakes to nucleate. However, evidence of spatially overlapping coseismic and postseismic slip challenges this vision. First, we infer shallow continuous slip in the first few kilometers below the surface before the three earthquakes, which are also seen to have reached the surface. Second, the May 2016 postseismic slip probably locates
on the coseismically ruptured circular patch at 5 km depth. Evidence comprise the post-seismic slip identified in the inversion of the July 2016 event in Figure 7 and the comparable coseismic and early-postseismic profiles in Figure 8b. Actually, in a rate-weakening regime, geometrical complexities alone promote a variety of slip rates in a continuous spectrum, from earthquakes to very slow events on the same fault segment [Romanet et al., 2018]. Furthermore, low effective normal stress from high pore fluid pressure could favor aseismic slip [Scholz, 1998]. We therefore argue that the coexistence of aseismic and seismic slip along the CF could be explained by a rate-weakening regime within a geometrically complex fault zone including complex fluid circulations. Such hypothesis now remains to be tested against numerical models.

Fault geometry and local topography are intrinsically related to local fault activity. The CF geometry is only known through its fault trace (Figure 5a, Section 2.1). The Nushki segment appears particularly straight compared to the rest of the fault, and can also be recognized by its low mountain front sinuosity [Crupa et al., 2017]. Together with other geomorphological indexes, this low roughness of the fault trace is interpreted as a sign of a tectonically more active fault segment by Ul-Hadi et al. [2013] and Crupa et al. [2017], an idea confirmed by the elevated surface slip rate along the Nushki segment. At smaller scale, the two peaks in surface slip (and seismic moment release) within the central bend segment are along rather straight fault portions (i.e. the derivative of the fault trace is flat) between 100-130 km and 147-170 km (Figure 5). Furthermore, sharp azimuth variations at about 95, 135 and 310 km north of Nushki coincide with three of the observed locked segments (zero surface slip and greater locking depth) separating the creeping segments. This is compatible with the idea that fault trace variations act as barriers to slip propagation [e.g. Wesnousky, 2006, Manighetti et al., 2015, Jolivet, Candela, et al., 2015, Van Rijsingen et al., 2019] or that they are the locations of significant off fault deformation [e.g. Okubo et al., 2019]. Actually, numerous subsidiary faults concentrate at bends in the fault system [Ruleman et al., 2007]. Those are mainly thrust faults, and seem to densify north of 31.3°N where the misorientation of the CF is accentuated by the eastward tilt of the fault.

In Section 4.2, we outlined the difference in slip regime south and north of the releasing bend at 180 km. Compared to the Nushki segment, the Qalat segment appears as geometrically more complex and more oblique to plate motion. This obliqueness of about 30° implies that the differential plate displacement projected in fault normal direction rises from 0±2 mm/yr south of 30°N to 15±2 mm/yr on the Qalat segment, according to DeMets et al. [2010] and Altamimi et al. [2017] (Figure S26). Thus, the zero fault normal displacement is likely to be a big approximation in this zone where both rigid plate rotation and geomorphological feature of non-rigid deformation evidence non-negligible fault normal motion with respect to the ~30 mm/yr fault-parallel motion.
## 6 Conclusion

We precisely describe slip patterns along the Chaman fault from InSAR time series covering 2015 to 2019. By integrating observations in space and time we outline five continuously creeping segments with maximum left-lateral slip rate in the shallow fault portion (< 2 km) reaching 5-10 mm/yr. This segmentation includes two prominent 80-km-long segments, notably the Nushki creeping segment which locates just south of the Central fault portion, also about 80-km-long, that hosted the most and largest earthquakes in the past-century. This description agrees with previous InSAR measurements from ALOS and Envisat. We estimate an upper bound for the seismic slip contribution to strain accommodation of about 3 mm/yr on average over 120 years. In our observation period, we observe the surface displacement induced by five slip events, and model the source at depth for the three events located on the Central Chaman fault. Significant aseismic slip is found in close proximity to those three earthquakes as pre-seismic slip and large afterslip, for a total induced aseismic strain release close to the one resulting from the earthquakes itself. For the M<sub>w</sub>5.6 May 2016 earthquake, part of the afterslip overlaps with the co-seismic rupture. Finally, we relate the variations in slip behavior along the Chaman fault at regional scale with the corresponding length-scales expressed in the fault geometry. The most striking feature is the correlation between the change in fault azimuth north of the Central portion and the change in left-lateral loading rate at depth from 12±3 mm/yr south to 7±2 mm/yr north. As a whole, the Chaman fault accommodates about one third of the differential plate motion. In order to assess long-term seismic hazard, the next step is to quantify the distribution of loading between the Chaman fault and nearby active faults of the plate boundary.

## 7 Acknowledgements

This work received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Geo-4D project, grant agreement 758210). Dr Angelique Benoit provided codes and methods at the foundation of this work. Dr Heresh Fattahi shared ENVISAT velocity maps and uncertainties used in Fattahi & Amelung [2016]. Raw data used in this study are freely available online. Synthetic Aperture Radar images are from the PEPS platform (scihub.copernicus.eu). The digital elevation model is from NASA EarthData. ERA-5 global reanalyses of atmospheric data are distributed by the ECMWF. Processed data represented in figures are archived on Zenodo (DOI:10.5281/zenodo.5221208). We thank both reviewers, the associate editor and the editor in chief for the extensive and complete reviews which greatly helped improving the initial manuscript. Many figures were created using Generic Mapping Tools [Wessel & Smith, 1998].

## References

Ahmad, W., Choi, M., Kim, S., & Kim, D. (2017). Detection of land subsidence due to excessive groundwater use varying with different land cover types in quetta valley, Pakistan using ESA-sentinel satellite data. *Nat Hazards Earth Syst Sci Discuss*, 1–21.

Altamimi, Z., Métivier, L., Rebischung, P., Rouby, H., & Collilieux, X. (2017). ITRF2014 plate motion model. *Geophysical Journal International*, 209(3), 1906–1912.

Alwahedi, M. A., & Hawthorne, J. C. (2019). Intermediate-magnitude postseismic slip follows intermediate-magnitude (Mw 4 to 5) earthquakes in California. *Geo-physical Research Letters*, 46(7), 3676–3687.

Ambraseys, N., & Bilham, R. (2003a). Earthquakes and associated deformation in northern Baluchistan 1892-2001. *Bulletin of the Seismological Society of America*, 93(4), 1573–1605.
Ambraseys, N., & Bilham, R. (2003b). Earthquakes in Afghanistan. *Seismological Research Letters*, *74*(2), 107–123.

Ambraseys, N., & Douglas, J. (2004). Magnitude calibration of north indian earthquakes. *Geophysical Journal International*, *159*(1), 165–206.

Avouac, J.-P. (2015). From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle. *Annual Review of Earth and Planetary Sciences*, *43*, 233–271.

Avouac, J.-P., Ayoub, F., Wei, S., Ampuero, J.-P., Meng, L., Leprince, S., ... Helmberger, D. (2014). The 2013, mw 7.7 balochistan earthquake, energetic strike-slip reactivation of a thrust fault. *Earth and Planetary Science Letters*, *391*, 128–134.

Barnhart, W. D. (2017). Fault creep rates of the Chaman fault (Afghanistan and Pakistan) inferred from InSAR. *Journal of Geophysical Research: Solid Earth*, *122*(1), 372–386.

Benoit, A., Pinel-Puysségur, B., Jolivet, R., & Lasserre, C. (2020, 03). CorPhU: an algorithm based on phase closure for the correction of unwrapping errors in SAR interferometry. *Geophysical Journal International*, *221*(3), 1959-1970. doi: 10.1093/gji/ggaa120

Berardino, P., Fornaro, G., Lanari, R., & Sansosti, E. (2002). A new algorithm for surface deformation monitoring based on small baseline differential SAR interferograms. *IEEE Transactions on Geoscience and Remote Sensing*, *40*(11), 2375-2383.

Bernard, M., Shen-Tu, B., Holt, W., & Davis, D. (2000). Kinematics of active deformation in the Sulaiman lobe and range, Pakistan. *Journal of Geophysical Research: Solid Earth*, *105*(B6), 13253–13279.

Beun, N., Bordet, P., & Carbonnel, J.-P. (1979). Premières données quantitatives relatives au coulissage du décrochement de Chaman (Afghanistan du sud-est). *Comptes rendus hebdomadaires des séances de l'Académie des sciences*, *288*, 931-934.

Bilham, R., Kakar, N. U., Kakar, D. M., Wang, K., Bürgmann, R., & Barnhart, W. D. (2019). The 1892 Chaman, Pakistan, earthquake. *Seismological Research Letters*, *90*(6), 2293–2303.

Bondár, I., & Storchak, D. (2011). Improved location procedures at the international seismological centre. *Geophysical Journal International*, *186*(3), 1220–1244.

Bürgmann, R. (2018). The geophysics, geology and mechanics of slow fault slip. *Earth and Planetary Science Letters*, *495*, 112–134. doi: 10.1016/j.epsl.2018.04.062

Çakir, Z., Ergintav, S., Özener, H., Dogan, U., Akoglu, A. M., Meghraoui, M., & Reilinger, R. (2012). Onset of aseismic creep on major strike-slip faults. *Geology*, *40*(12), 1115–1118.

Chartrand, R. (2011). Numerical differentiation of noisy, nonsmooth data. *International Scholarly Research Notices*, *2011*.

Crupa, W. E., Khan, S. D., Huang, J., Khan, A. S., & Kasi, A. (2017). Active tectonic deformation of the western indian plate boundary: A case study from the Chaman fault system. *Journal of Asian Earth Sciences*, *147*, 452–468.

Dalaison, M., & Jolivet, R. (2020). A Kalman filter time series analysis method for InSAR. *Journal of Geophysical Research: Solid Earth*, *125*(7), e2019JB019150.

Dal Zilio, L., Lapusta, N., & Avouac, J.-P. (2020). Unraveling scaling properties of slow-slip events. *Geophysical Research Letters*, *47*(10), e2020GL087477.

DeMets, C., Gordon, R. G., & Argus, D. F. (2010). Geologically current plate motions. *Geophysical Journal International*, *181*(1), 1–80.

Dewey, J. W., Bergman, E. A., Hopper, M. G., & Sipkin, S. A. (2006). *Seismicity of Afghanistan and vicinity* (Open-File Report No. 1185). U.S. Geological Survey.

De Zan, F., Zonno, M., & López-Dekker, P. (2015). Phase inconsistencies and multiple scattering in sar interferometry. *IEEE Transactions on Geoscience and Remote
Sensing, 53 (12), 6608–6616.

Dieterich, J. H. (1979). Modeling of rock friction: 1. experimental results and con- stitutive equations. *Journal of Geophysical Research: Solid Earth, 84*(B5), 2161– 2168. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10 .1029/JB084iB05p02161 doi: 10.1029/JB084iB05p02161

Doin, M.-P., Guillaaso, S., Jolivet, R., Lasserre, C., Lodge, F., Ducret, G., & Grandin, R. (2011). Presentation of the small baseline NSBAS processing chain on a case example: The Etna deformation monitoring from 2003 to 2010 using Envisat data. In *Proceedings of the fringe symposium* (pp. 3434–3437).

Duputel, Z., Agram, P. S., Simons, M., Minson, S. E., & Beck, J. L. (2014). Ac- counting for prediction uncertainty when inferring subsurface fault slip. *Geophysica l Journal International, 197*(1), 464–482.

Dziewonski, A., Chou, T.-A., & Woodhouse, J. (1981). Determination of earth- quake source parameters from waveform data for studies of global and regional seismicity. *Journal of Geophysical Research: Solid Earth, 86*(B4), 2825–2852.

Ekström, G., Nettles, M., & Dziewoński, A. (2012). The global CMT project 2004– 2010: Centroid-moment tensors for 13,017 earthquakes. *Physics of the Earth and Planetary Interiors, 200*, 1–9.

Elliott, J., Jolivet, R., González, P. J., Avouac, J.-P., Hollingsworth, J., Searle, M., & Stevens, V. (2016). Himalayan megathrust geometry and relation to topography revealed by the gorkha earthquake. *Nature Geoscience, 9*(2), 174–180.

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., ... others (2007). The shuttle radar topography mission. *Reviews of geophysics, 45*(2).

Fattahi, H., & Amelung, A. (2016). InSAR observations of strain accumulation and fault creep along the Chaman fault system, Pakistan and Afghanistan. *Geophys. Res. Lett.*, 43, 8399–8406. doi: 10.1002/2016GL070121

Fattahi, H., Amelung, F., Chaussard, E., & Wdowinski, S. (2015). Coseismic and postseismic deformation due to the 2007 M5.5 Ghazaband fault earthquake, Balochistan, Pakistan. *Geophysical Research Letters, 42*(9), 3305–3312.

Furuya, M., & Satyabala, S. (2008). Slow earthquake in Afghanistan detected by In- SAR. *Geophysical Research Letters, 35*(6).

Goldstein, R. M., & Werner, C. L. (1998). Radar interferogram filtering for geophys- ical applications. *Geophysical Research Letters, 25*(21), 4035–4038.

Goldstein, R. M., Zebker, H. A., & Werner, C. L. (1988). Satellite radar interferome- try: Two-dimensional phase unwrapping. *Radio science, 23*(4), 713–720.

Gurrola, E., Rosen, P., Sacco, G., Seliga, W., Zebker, H., Simons, M., & Sandwell, D. (2010). InSAR scientific computing environment. In *2010 american geophysical union meeting*.

Gutenberg, B., & Richter, C. F. (1944). Frequency of earthquakes in California. *Bul- letin of the Seismological Society of America, 34*(4), 185–188.

*Isc on-line bulletin*. (2020). International Seismological Centre. doi: 10.31905/ D808B830

Jolivet, R., Agram, P. S., Lin, N. Y., Simons, M., Doin, M.-P., Peltzer, G., & Li, Z. (2014). Improving InSAR geodesy using global atmospheric models. *Journal of Geophysical Research: Solid Earth, 119*(3), 2324–2341.

Jolivet, R., Candela, T., Lasserre, C., Renard, F., Klinger, Y., & Doin, M.-P. (2015). The burst-like behavior of aseismic slip on a rough fault: The creeping section of the Haiyuan fault, ChinaShort note. *Bulletin of the Seismological Society of America, 105*(1), 480–488.

Jolivet, R., Duputel, Z., Riel, B., Simons, M., Rivera, L., Minson, S., ... others (2014). The 2013 Mw 7.7 Balochistan earthquake: seismic potential of an accre- tionary wedge. *Bulletin of the Seismological Society of America, 104*(2), 1020– 1030.
Jolivet, R., & Frank, W. (2020). The transient and intermittent nature of slow slip. *AGU Advances*, 1(1), e2019AV000126.

Jolivet, R., Grandin, R., Lasserre, C., Doin, M.-P., & Peltzer, G. (2011). Systematic InSAR tropospheric phase delay corrections from global meteorological reanalysis data. *Geophysical Research Letters*, 38(17).

Jolivet, R., Lasserre, C., Doin, M.-P., Peltzer, G., Avouac, J.-P., Sun, J., & Dailu, R. (2013). Spatio-temporal evolution of aseismic slip along the Haiyuan fault, China: Implications for fault frictional properties. *Earth and Planetary Science Letters*, 377, 23–33.

Jolivet, R., Simons, M., Agram, P., Duputel, Z., & Shen, Z.-K. (2015). Aseismic slip and seismogenic coupling along the central San Andreas fault. *Geophysical Research Letters*, 42(2), 297–306.

Jones, A., et al. (1961). Reconnaissance geology of part of west Pakistan. *A Colombo Plan Cooperative Project. Toronto, Canada: Government of Canada*.

Kakar, N., Kakar, D. M., Khan, A. S., & Khan, S. D. (2016). Land subsidence caused by groundwater exploitation in Quetta valley, Pakistan. *International Journal of Economic and Environmental Geology*, 7(2), 10–19.

Kaneko, Y., Avouac, J.-P., & Lapusta, N. (2010). Towards inferring earthquake patterns from geodetic observations of interseismic coupling. *Nature Geoscience*, 3(5), 363–369.

Kaneko, Y., Fialko, Y., Sandwell, D., Tong, X., & Furuya, M. (2013). Interseismic deformation and creep along the central section of the North Anatolian fault (Turkey): InSAR observations and implications for rate-and-state friction properties. *Journal of Geophysical Research: Solid Earth*, 118(1), 316–331.

Kreemer, C., Blewitt, G., & Klein, E. C. (2014). A geodetic plate motion and global strain rate model. *Geochemistry, Geophysics, Geosystems*, 15(10), 3849–3889.

Lawrence, R. D., Khan, S. H., & Nakata, T. (1992). Chaman fault, Pakistan-Afghanistan. *Ann. Tectonicae*, 6, 196–223.

Lawrence, R. D., & Yeats, R. S. (1979). Geological reconnaissance of the Chaman fault in Pakistan. *Geodynamics of Pakistan*, 351–357.

Lawrence, R. D., Yeats, R. S., Khan, S. H., Subhani, A. M., & Bonelli, D. (1981). Crystalline rocks of the Spinatizha area, Pakistan. *Journal of Structural Geology*, 3(4), 449–457. (Tectonics) doi: https://doi.org/10.1016/0191-8141(81)90044-4

Lindsey, E. O., Fialko, Y., Bock, Y., Sandwell, D. T., & Bilham, R. (2014). Localized and distributed creep along the southern San Andreas fault. *Journal of Geophysical Research: Solid Earth*, 119(10), 7909–7922.

Liu, F., Elliott, J., Craig, T., Hooper, A., & Wright, T. (2021). Improving the resolving power of InSAR for earthquakes using time series: A case study in Iran. *Geophysical Research Letters* (48), e2021GL093043.

Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. *Geochemistry, Geophysics, Geosystems*, 6(1).

Manighetti, I., Campillo, M., Bouley, S., & Cotton, F. (2007). Earthquake scaling, fault segmentation, and structural maturity. *Earth and Planetary Science Letters*, 253(3), 429–438. doi: https://doi.org/10.1016/j.epsl.2006.11.004

Manighetti, I., Caulet, C., De Barros, L., Perrin, C., Cappa, F., & Gaudemer, Y. (2015). Generic along-strike segmentation of a far normal faults, e ast a frica: Implications on fault growth and stress heterogeneity on seismogenic fault planes. *Geochemistry, Geophysics, Geosystems*, 16(2), 443–467.

Marone, C. (1998). Laboratory-derived friction laws and their application to seismic faulting. *Annual Review of Earth and Planetary Sciences*, 26(1), 643–696.

Maurer, J., & Johnson, K. (2014). Fault coupling and potential for earthquakes on the creeping section of the central San Andreas fault. *Journal of Geophysical Research: Solid Earth*, 119(5), 4414–4428.
Michel, S., Avouac, J.-P., Jolivet, R., & Wang, L. (2018). Seismic and aseismic moment budget and implication for the seismic potential of the parkfield segment of the San Andreas fault. *Bulletin of the Seismological Society of America*, **108**(1), 19–38.

Michel, S., Gualandi, A., & Avouac, J.-P. (2019). Similar scaling laws for earthquakes and cascadia slow-slip events. *Nature*, **574**(7779), 522–526.

Mohadjer, S., Bendick, R., Ischuk, A., Kuzikov, S., Kostuk, A., Saydullaev, U., ... others (2010). Partitioning of India-Eurasia convergence in the Pamir-Hindu Kush from GPS measurements. *Geophysical Research Letters*, **37**(4).

Okubo, K., Bhat, H. S., Rougier, E., Marty, S., Schubnel, A., Lei, Z., ... Klinger, Y. (2019). Dynamics, radiation, and overall energy budget of earthquake rupture with coseismic off-fault damage. *Journal of Geophysical Research: Solid Earth*, **124**(11), 11771–11801.

Perfettini, H., & Avouac, J.-P. (2004). Postseismic relaxation driven by brittle creep: A possible mechanism to reconcile geodetic measurements and the decay rate of aftershocks, application to the Chi-Chi earthquake, Taiwan. *Journal of Geophysical Research: Solid Earth*, **109**(B2).

Petersen, M. D., Dawson, T. E., Chen, R., Cao, T., Wills, C. J., Schwartz, D. P., & Frankel, A. D. (2011). Fault displacement hazard for strike-slip faults. *Bulletin of the Seismological Society of America*, **101**(2), 805–825.

Radiguet, M., Cotton, F., Vergnolle, M., Campillo, M., Valette, B., Kostoglodov, V., & Cotte, N. (2011). Spatial and temporal evolution of a long term slow slip event: the 2006 Guerrero slow slip event. *Geophysical Journal International*, **184**(2), 816–828.

Romanet, P., Bhat, H. S., Jolivet, R., & Madariaga, R. (2018). Fast and slow slip events emerge due to fault geometrical complexity. *Geophysical Research Letters*.

Ruleman, C., Crone, A., Machette, M., Haller, K., & Rukstales, K. (2007). *Map and database of probable and possible Quaternary faults in Afghanistan* (Open-File Report No. 1103). U.S. Geological Survey.

Salvatier, J., Wiecki, T. V., & Fonnesbeck, C. (2016). Probabilistic programming in python using PyMC3. *PeerJ Computer Science*, **2**, e55.

Savage, J., & Burford, R. (1973). Geodetic determination of relative plate motion in central California. *Journal of Geophysical Research*, **78**(5), 832–845.

Sborshchikov, I., Savostin, L., & Zonenshain, L. (1981). Present plate tectonics between Turkey and Tibet. *Tectonophysics*, **79**(1-2), 45–73.

Scholz, C. H. (1998). Earthquakes and friction laws. *Nature*, **391**(6662), 37–42.

Scordilis, E. (2006). Empirical global relations converting MS and Mb to moment magnitude. *Journal of seismology*, **10**(2), 225–236.

Segall, P. (2010). Earthquake and volcano deformation. In (chap. 2 Dislocation Models of Strike-Slip Faults). Princeton University Press.

Steinbrugge, K. V., Zacher, E. G., Tocher, D., Whitten, C., & Claire, C. (1960). Creep on the San Andreas fault. *Bulletin of the Seismological Society of America*, **50**(3), 389–415.

Szeliga, W., Bilham, R., Kakar, D. M., & Lodi, S. H. (2012). Interseismic strain accumulation along the western boundary of the Indian subcontinent. *Journal of Geophysical Research: Solid Earth*, **117**(B08404). doi: 10.1029/2011JB008822

Szeliga, W., Bilham, R., Schelling, D., Kakar, D. M., & Lodi, S. (2009). Fold and thrust partitioning in a contracting fold belt: Insights from the 1931 Mach earthquake in Baluchistan. *Tectonics*, **28**(5).

Tapponnier, P., Mattauer, M., Proust, F., & Cassaigneau, C. (1981). Mesozoic ophiolites, sutures, and large-scale tectonic movements in Afghanistan. *Earth and Planetary Science Letters*, **52**(2), 355–371.

Thomas, M. Y., Avouac, J.-P., & Lapusta, N. (2017). Rate-and-state friction properties of the Longitudinal Valley fault from kinematic and dynamic modeling of
seismic and aseismic slip. *Journal of Geophysical Research: Solid Earth*, **122**(4), 3115–3137.

Tymofyeyeva, E., Fialko, Y., Jiang, J., Xu, X., Sandwell, D., Bilham, R., ... others (2019). Slow slip event on the southern San Andreas fault triggered by the 2017 Mw 8.2 Chiapas (Mexico) earthquake. *Journal of Geophysical Research: Solid Earth*, **124**(9), 9956–9975.

Ul-Hadi, S., Khan, S. D., Owen, L. A., & Khan, A. S. (2013). Geomorphic response to an active transpressive regime: A case study along the Chaman strike-slip fault, western Pakistan. *Earth Surface Processes and Landforms*, **38**(3), 250–264.

Van Rijtsingen, E., Funiciello, F., Corbi, F., & Lallemand, S. (2019). Rough subducting seafloor reduces interseismic coupling and mega-earthquake occurrence: Insights from analogue models. *Geophysical Research Letters*, **46**(6), 3124–3132.

Wells, D. L., & Coppersmith, K. J. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bulletin of the seismological Society of America*, **84**(4), 974–1002.

Wesnousky, S. G. (2006). Predicting the endpoints of earthquake ruptures. *Nature*, **444**(7117), 358–360.

Wessel, P., & Smith, W. H. (1998). New, improved version of Generic Mapping Tools released. *Eos, Transactions American Geophysical Union*, **79**(47), 579–579.

Wheeler, R. L., Bufe, C. G., Johnson, M. L., Dart, R. L., & Norton, G. (2005). *Seismotectonic map of Afghanistan, with annotated bibliography* (Open-File Report No. 1264). U.S. Geological Survey.

Wright, T. J., Parsons, B. E., & Lu, Z. (2004). Toward mapping surface deformation in three dimensions using InSAR. *Geophysical Research Letters*, **31**(1).

Yeats, R. (2012). Active faults of the world. In (chap. 8 India, the Himalaya, Mainland China, and Central Asia). Cambridge university press, New York.

Zhu, L., & Rivera, L. A. (2002). A note on the dynamic and static displacements from a point source in multilayered media. *Geophysical Journal International*, **148**(3), 619–627.