The frequency of earthquakes that have occurred along a fault ultimately determines its seismic potential and is therefore important for hazard assessment. A direct way to reveal the number of surface rupture events is palaeoseismic trenching (Deng and Liao, 1996; Marco et al., 2005; Rockwell and Ben-Zion, 2007), but it is often limited by difficulties in logistics in the field. Since earthquakes are recorded on the earth’s surface in the form of offset geomorphic features, such offsets can be used as an alternative method to characterise repeating earthquakes and fault slip behaviour. Under the assumption that similar-aged offset geomorphic features experienced the same number of earthquakes, offset measurements appear as groups, which can be used as an indicator of repeating earthquakes (Sieh, 1978; McGill and Sieh, 1991). This method, now known as “cumulative offset probability distribution, COPD”, requires a large amount of offset measurements for statistical analysis (e.g., Zielke et al., 2015; Bi et al., 2018). During the last decade, high-resolution remote-sensing data promote the application of the “COPD” method in various places around the world (e.g., Zielke et al., 2010; Klinger et al., 2011; Li et al., 2012; Ren et al., 2016; Bi et al., 2018; Choi et al., 2018). For instance, Zielke et al. (2010) measured 149 offsets along the Carrizo segment of the San Andreas fault using light detection and ranging (LiDAR) data, and identified five peaks (5.3, 9.8, 14.9, 20.1, and 24.5 m, respectively) from the offset measurements, with the first peak corresponding to the slip of the 1857 M _w 7.9 Fort Tejon earthquake. Klinger et al. (2011) measured five slip peaks that are multiples of 6 m along the Fuyun fault in the south of the Altay Mountains, Northeastern Xinjiang, China, based on Quickbird images, and concluded that ruptures on the Fuyun fault exhibits characteristic slip behaviour. Nonetheless, a recent study by Lin et al. (2020) pointed out that the COPD method is only applicable to a few large mature strike-slip faults such as the San Andreas and Fuyun faults.

Surface offsets are important for estimating the magnitude of slip in past earthquakes. Several early studies, e.g., Arrowsmith et al. (1998); Avouac (1993); Colman and Watson (1983); Culling (1960); Jayangondaperumal et al. (2013); Jayangondaperumal et al. (2017), have noted that the morphology of a scarp is useful for constraining the evolution time since its formation. However, previous studies based on the COPD method depend solely on the offset itself, and ignore the shape of the geomorphic features. In this study, we propose a new method to constrain earthquake history by modelling scarp degradation. The method is applied to the Muztagh Ata and Tahman faults in the Chinese Pamir, providing new information about regional earthquake records.

The Pamir, located in the western end of the Indo-Asian collision zone, is a tectonically active region in central Asia (Arnaud et al., 1993; Brunel et al., 1994; Burtman and Molnar, 1993; Cao et al., 2013; Schurr et al., 2014; Sobel et al., 2011) (Figure 1) The deformation of the eastern Pamir is dominated by E-W extension along the Kongur Shan (Robinson et al., 2004; Yuan et al., 2013), a 250-km-long NW-SE-striking fault system, which can be divided into four segments from north to south: the right-lateral strike-slip Muji fault, the Kongur Shan and Muztagh Ata normal faults, the Tahman normal fault, and the Tashkurgan normal fault (Robinson et al., 2004; Thiede et al., 2013).

In this study, we focus on the Muztagh Ata and Tahman faults in the southern part of the Kongur Shan extensional system (Figure 1). The Muztagh Ata fault, a 70-km-long NW-SE-striking normal fault, forms the western boundary of the Muztagh Ata Mountain. Chevalier et al. (2015) combined geomorphic analysis and ¹⁰ Be dating of alluvial terrace risers to determine the late Quaternary vertical slip rate of the southern Muztagh Ata fault. They obtained a rate of ∼1.7 mm/yr, corresponding to a NE-SW extension rate of 1.3–2.0 mm/yr for a fault dip of 35°–45° (Chevalier et al., 2015) (Figure 1). The Tahman fault, extending 20 km north-south, is separated into two segments by a 0.5-km step-over zone (Li et al., 2011). The vertical slip rate was estimated to be 2.5±0.1 mm/yr from ¹⁰ Be dating of alluvial terraces (Li, 2014) (Figure 1).

Both the Muztagh Ata and Tahman faults have been seismically active since the Quaternary (Yuan et al., 2013; Li, 2014; Chevalier et al., 2015). The USGS earthquake catalogue shows that hundreds of M _w ≥ 2 earthquakes occurred in the region since 1900, including 46 M _w ≥ 5 earthquakes. The largest event that occurred in the Tashkurgan Valley is the 1895/07/05 M _s 7.0 earthquake, which reached a maximum intensity of IX (Li et al., 2011). Up to now, the fault scarps created by the 1895 M _s 7.0 Tashkurgan earthquake are well-preserved along the southern Muztagh Ata and Tahman faults (Li et al., 2011; Chevalier et al., 2015). Nonetheless, due to difficulties in palaeoseismic trenching and sampling in the area, the regional earthquake history remains largely unknown.

The scarps along the Muztagh Ata and Tahman faults are potentially useful for constraining the timing of surface rupture events. We use the well-established scarp degradation model (Culling, 1960, 1963; Kokkalas and Koukouvelas, 2005; Pelletier and Cline, 2007) to analyse the relationship between scarp morphology and its age. The derivation of the model can be found in many studies, (e.g., Andrews and Bucknam, 1987; Roering et al., 1999; Pelletier and Cline, 2007; Xu et al., 2021). Here, we briefly introduce the set-up of the model.

Sediment flux S(x) down surface topography h is related to local slope d h d x and a diffusion coefficient f(h): S x = f h ⋅ d h d x (1) where we take a nonlinear diffusion coefficient f(h) that follows the general form (Andrews and Bucknam, 1987; Pelletier and Cline, 2007): f h = k 1 − min c ; d h d x / S c n (2) where k is a constant, S _c is the critical slope, n is an empirical exponent, and c is a threshold to limit the maximum ratio between the hillslope and the critical slope S _c .

The continuity equation is used to relate sediment flux and the evolution of topography with time ( ∂ h ∂ t ) : − ρ s ∂ h ∂ t = ρ s ▽ ⋅ S x + ρ r C 0 (3) where ρ _s and ρ _r are the bulk densities of sediment and rock, and C ₀ is the subsidence/uplift rate.

Combining Equations 1–3, and ignoring the subsidence/uplift term ρ _r C ₀, we derive the nonlinear model for scarp degradation: − ∂ h ∂ t = ▽ ⋅ k 1 − min c , d h d x / S c n d h d x (4) The method is best suitable for normal fault scarps, the shape of which is less affected by gravity. Reliable estimates of the geomorphological parameters in the diffusion model (k, S _c , n, and c) need to be obtained first by local geological survey. Xu et al. (2021) have analysed the fluvial terrace riser degradation along the Kongur Shan fault ( ∼130 km north of our study area, yellow star in Figure 1) in great detail. They tested a variety of parameters in the diffusion model, and obtained the best-fit by comparing the ages derived from morphological modelling to those from geological dating. We use their best-fitting parameters (k = 1, n = 2, c = 0.99, S _c = tan 33°) in our study because the sediment transport processes should be similar given the close distance. The geometry of a scarp at its formation is determined by the initial slope angle and vertical offset. Also, a scarp degrades faster with a steeper slope. The typical vertical offset is associated with earthquake slip. We use an initial slope angle of 35° [average value from Li. (2014)] and a typical vertical offset of 2.5 m from a single event [from Li et al. (2011)].

Using Equation 4, we model the evolution of fault scarp with time. We assume repeating earthquakes every 2 kyrs, i.e., earthquakes occur at 0, 2, 4, 6 kyr with a vertical offset of 2.5 m (Figures 2A,D,G,J). This value is chosen based on the geomorphological dating in Chevalier et al. (2015), which suggests that the oldest geomorphic surface along the Muztagh Ata fault (∼20 kyrs) may have recorded ∼10 large earthquakes. We calculate the slope (first derivative of topography with respect to horizontal distance, d h d x ) (Figures 2B,E,H,K), as well as the second derivative d 2 h d x 2 in order to remove the regional topographic gradient of the natural earth’s surface (Figures 2C,F,I,L). As shown in Figure 2, the most significant feature is that the width of the scarp widens. We find that the second derivative of the scarp can be described by a modified Gaussian pulse function: g x = e 1 2 ⋅ x σ ⋅ exp − x σ 4 / 2 (5) where x is the distance, and σ is an exponential factor to be determined.

We use the Gaussian pulse function to define the width of fault scarp: the distance between two points with g (x) equals to 10⁻⁴ and −10⁻⁴ ( g (x) = 0 gives an infinite x, grey line in Figure 2C). The second derivative is advantageous over topography and slope (the first derivative) because it provides a mathematical description of the scarp morphology and therefore avoids artificial identification of scarp end points. The scarp width (also referred to as the bandwidth of the best-fitting Gaussian pulse function) increases from the initial value of 9.3 m after the first event, to 19.0 m after the second event, to 26.9 m after the third event, and to 34.7 m following the fourth event. We plot the bandwidth-time in Figure 3 and find that the changes in scarp widths can be fitted well with a 1st (black dash line) and 2nd order (red solid line) polynomial. We did not use the 1st order polynomial given that the degradation process is unlikely to be at a constant rate, but should slow down with increasing widths.

To investigate how the initial parameters in the model including earthquake recurrence intervals and offsets affect the time-width relationship, we simulated four different earthquake modes: 1) offset = 2 m, recurrence interval = 2 kyrs; 2) offset = 2.5 m, recurrence interval = 2 kyrs; 3) offset = 2.5 m, recurrence interval = 3 kyrs; 4) offsets = 2.5, 3.5, 2.5, and 1.5 m, with recurrence intervals of 3, 2, and 1 kyr, respectively. The time-width relationships of the four earthquake modes are shown in Figure 4.

Comparing modes 2 (red line) and 3 (green line), we noted that scarp width increases more rapidly with shorter recurrence intervals. For a scarp width of 35.2 m (the maximum scarp width observed in this study), the inferred ages from modes 2 and 3 are 6.0 kyrs before the present (BP) and 7.1 kyrs BP, respectively. Given that most scarps in our study area are narrower than 30 m, the age uncertainties caused by the initial recurrence interval in the model set-up should be within ∼1 kyr. The difference between modes 1 (purple line) and 2 (red line) in Figure 4 shows that larger offsets create wider scarps. Based on the time-width relationships, the age estimates of the widest scarp (35.2 m) are 7.6 kyrs BP (mode 1) and 6.0 kyrs BP (mode 2).

The actual earthquake rupture history is often more complicated, with different slips and recurrence intervals (Zielke et al., 2015). To test this scenario, we vary the offsets from 1.5 to 3.5 m, and recurrence intervals from 1 kyr to 3 kyrs. Interestingly, the time-width relationship of mode 4 is similar to that from mode 2, the long-term fault slip rates of which are equal (1.25 mm/yr). This means that the change in width of a fault scarp is ultimately determined by the fault slip rate. The simulations also suggest that although the time-width relationship varies with the assumed earthquake offsets and recurrence intervals, the age uncertainty for scarps of narrower than 30 m seems to be less than 1.4 kyr. By analysing a large amount of fault scarps, we may further bring down the uncertainty.

In this study, we proposed a new approach of using fault scarp degradation model to constrain earthquake history, with applications to the Muztagh Ata and Tahman faults in the Chinese Pamir. We mapped fault scarps based on WorldView-2 stereo imagery and determined the widths of the scarps using a modified Gaussian pulse function. The measurements show that the widths of the scarps are clustered in four groups, each group representing a surface rupture event. Given the empirical time-width relationship, we calculated the timing of the earthquakes, with the youngest cluster, 0.1±0.4 kyrs BP, corresponding to the 1895 Tashkurgan earthquake. Based on the inferred earthquake recurrence intervals and offsets, we estimated the vertical slip rate (2.2±0.3 mm/yr) of the Muztagh and Tahman faults, in agreement with the geological estimates of 1.7–2.5 mm/yr from Chevalier et al. (2015); Li. (2014). The consistent results suggest that the proposed method provides an innovative means for investigating earthquake history and can be potentially applied to other active faults where surface scarps are well-preserved.