The Himalayas are the resultant of the colliding interaction of the Indian and Eurasian plates during the late Cretaceous to early Tertiary periods (Wells, 1984; Yeats et al., 1984; Yeats and Hussain, 1987; McDougall and Khan, 1990; McDougall and Hussain, 1991; Pogue et al., 1992; Smith et al., 1994; Beck et al., 1995; Jaswal et al., 1997; Yin and Harrison, 2000; Khan and Glenn, 2006). Researchers have long been intrigued by this collision product to investigate the traits that change during the earthquake cycles. Several devastating earthquakes have resulted from this collision (Abbas et al., 2022). The suture zone formed by the Indian plate subducting beneath the Eurasian plate is dominated by numerous thrusting features in the southern region (Pudsey, 1986; Pogue et al., 1992; Bender and Raza, 1995). This region has also experienced the introduction of significant strike–slip structures due to the shoulder-rubbing nature of the collision between two plates at some points (Yeats et al., 1984; Jaswal et al., 1997; Yin and Harrison, 2000; Khan and Glenn, 2006; Abbas et al., 2022). Although thrusting and folding account for the majority of the deformation in the northwest Himalayas (McDougall and Khan, 1990; McDougall and Hussain, 1991; Jaswal et al., 1997; Khan and Glenn, 2006), the existence of strike–slip motion in the area makes the region’s tectonic system challenging to comprehend. The Punjab foreland basin along the northwestern Himalayan frontal thrust system is underlain by the Salt, Surghar, and Trans Indus ranges as a result of the progressive deformation brought on by the continent–continent type of collision (Baker et al., 1988; Blisniuk et al., 1998; Chen and Khan, 2009). Due to progressive tectonic processes of thrusting and folding, known structural reentrants in the region have displayed relatively substantial lateral variation in the deformation throughout the ranges that are under the process of ongoing growth (Houseman and England, 1986; Pudsey, 1986; Yeats and Hussain, 1987; Baker et al., 1988; Pogue et al., 1992; Smith et al., 1994; Blisniuk et al., 1998). However, a deep analysis of the components that counterbalance the frontal thrust system is lacking (Chen and Khan, 2010).

A major strike–slip structure in this region is the Kalabagh Fault (Yeats et al., 1984; Baker et al., 1988; McDougall and Khan, 1990; Blisniuk et al., 1998; Chen and Khan, 2010), shown in Figure 1. This structure truncates the Salt Range on the western end (Baker et al., 1988; McDougall and Khan, 1990; Abbas et al., 2022), one of the largest salt-hosting formations in the world. This fault has been studied from different perspectives recently. Yeats et al. (1984) claimed that Ghundi underwent Quaternary deformation. Long-term displacement along the Kalabagh Fault by piercing point data was used by McDougall and Khan (1990) to compute the average slip rate. Chen and Khan calculated the slip rate in the Kalabagh Fault Zone using two descending SAR images of ERS-1 and ERS-2 satellite data in 2010. Because of the available data source constraints and the use of a single viewing geometry, Chen and Khan (2010) could only compute the slip rate in the northern section of the Kalabagh Fault. Sentinel-1A, the most recent mission of the Interferometric Synthetic Aperture Radar (InSAR), is run by the European Space Agency and provides much more precise data to evaluate the creep rate across the fault (Savage and Lisowski, 1993; Lee et al., 2003; Titus et al., 2006; Cavalié et al., 2008). Sentinel-1A provides improved uncertainty of about 1 mm/yr over 100 km and better temporal and spatial coverage (Bürgmann et al., 2000; Hsieh et al., 2011; Jolivet et al., 2012; Zheng et al., 2017; Ng et al., 2018; Rosi et al., 2018) for the estimation of the slip rates across the whole Kalabagh Fault accurately.

Numerous studies have used InSAR data successfully for comprehending transient strain accumulation (Peltzer et al., 2001; Walters et al., 2013), creeping in the San Andreas Fault (Rosen et al., 1998; Johanson and Bürgmann, 2005), and slip rates on key faults of western Tibet (Chen et al., 2018). To estimate the slip rate accurately for active faults in a given area, which typically exhibit very slow slip rates (few mm/year), InSAR measurement can be utilized due to its exceptional signal-to-noise ratio (Peltzer et al., 2001; Wright et al., 2004; Xu et al., 2016; Yu et al., 2017).

More precise measurements of deformation and slip rate along the Kalabagh Fault are necessary in addition to the existing studies to comprehend the fault behavior. Previous studies on the Kalabagh Fault also did not account for atmospheric noise in their analyses. We employed the small baseline subset (NSBAS) technique and atmospheric correction through Generic Atmospheric Correction Online Service for InSAR (GACOS) to determine the rate map. Data of both ascending and descending tracks for 7 years, from 2015 to 2022, were analyzed to ensure comprehensive coverage and minimize the effect of atmospheric noise in the estimation of slip rates. In the end, we simulated the Coulomb stress changes due to the creeping of the KBF on nearby local strike–slip and thrust faults.

An established method for measuring ground deformation brought on by the buildup of interseismic strain is InSAR time-series analysis (Fialko, 2006; Jolivet et al., 2013; Biggs and Wright, 2020; Weiss et al., 2020). In several instances, InSAR has proven its capability to estimate as low slip rates as a few millimeters in a year (Fialko, 2006; Prati et al., 2010; Tong et al., 2013; Mousavi et al., 2015; Xu et al., 2016; Yu et al., 2017; Weiss et al., 2020). Several techniques and methodologies have been developed to compute the rate map of active fault systems (Berardino et al., 2003; Werner et al., 2003). In the current work, the deformation along the Kalabagh Fault was estimated using the SBAS–InSAR approach. This technique has significantly improved measurements of crustal deformation, for example, by establishing the spatio-temporal baseline of the chosen interferograms, which produces coherent deformation time series (Castellazzi et al., 2016). The SBAS InSAR time-series technique lessens the phase decorrelation effects by adding a filter to signals with strong spatial correlation but lower temporal correlation (Berardino et al., 2003; Hooper, 2008; Ferretti et al., 2011).

With full images of Pakistan being captured on average every 12 days, the Sentinel-1, C-band SAR satellites from the European Space Agency (ESA) offer extraordinary temporal data coverage for interferometry. Interseismic slip rates on various faults worldwide have been calculated using InSAR time-series methods (Fialko et al., 2002; Chen and Khan, 2009; Satyabala et al., 2012; Fattahi and Amelung, 2016). Though InSAR can estimate the slip rate and measure the locking depth of the significant KBF more precisely, it has not been utilized thus far to quantify the interseismic motion. Between 2015 and 2022, we used Sentinel-1A ascending orbit data from 175 imageries and descending data from 156 imageries, comprising 585 ascending and 554 descending interferograms, which were acquired from https://comet.nerc.ac.uk/COMET-LiCS-portal (last opened on 31-08-2022) to calculate the interseismic deformation across the KBF. As seen in Supplementary Figure S1, selected data pairs are characterized by a spatial baseline of 150 m. The temporal distribution of interferograms for the whole study period against a spatial baseline is shown in Supplementary Figure S1. The upper section shows the interferogram network for the ascending track T: 173 (LiCSAR Frame ID: 173A_05749_131313), whereas the lower section is for the descending track T: 005 (LiCSAR Frame ID: 005D_05797_131313).

The LicSBAS code (Morishita et al., 2020; Morishita, 2021) was employed for temporal analysis using coherence and unwrapped interferograms. It employs the SNAPHU algorithm to automatically unwrap the data (Chen and Zebker, 2002). We masked all pixels having a coherence of 0.1 before employing the time-series analysis to minimize the impact added by the unwrapping errors. We kept the unwrapped coverage to be less than 0.5 and the cut-off threshold for coherence to 0.06. After unwrapping the interferograms, resampling and geocoding were carried out through a better-resolution digital elevation model (DEM) provided by the SRTM. LicSBAS uses NSBAS (modified small baseline technique) and has the advantage of removing weak coherent interferograms through its loop closure algorithm (Ghorbani et al., 2022). The entire interferometric network is then inverted using the least-squares method for incremental displacements between the acquisition dates to estimate the displacement value of every pixel along the fault (KBF). Before estimating the rate map, interferograms were corrected for atmospheric noise by applying the GACOS service (Morishita et al., 2020).

Additionally, spatial filtering of 2 km (low-pass filtering) and temporal baseline variation of 47 days (high-pass filtering) employing Gaussian filter Kernel (Hooper and Zebker, 2007) were used to differentiate the noise from the temporal displacement variations. To eliminate the linear tendency in output time series, bilinear de-ramping was applied. In the end, the effect of topography error was removed from the final line-of-sight (LOS) velocity field. However, the deformation was estimated only in the LOS direction. Afterward, the LOS deformation vector was decomposed to fault-parallel and vertical components. Deformation along the LOS can be estimated by decomposing the LOS velocities into parallel, perpendicular, and vertical components to the fault velocities: D L O S = sin θ cos α   sin θ sin α − cos θ V E V N V U . (1)

DLOS marks the movement along the LOS toward the satellite, θ denotes the incidence angle, and the azimuth angle for the LOS vector is denoted by α. V_N denotes the vector decomposed in the north direction, V_E marks the vector in the East direction, and V_U denotes the vertical direction vector.

For resolving the horizontal and vertical components of the approximately north–south KBF, fault plane and fault geometry are prerequisites to compute the strike–slip and dip–slip values. In order to do so, the following relation (Chen and Khan, 2010; Lindsey and Fialko, 2013; Watson et al., 2022) of slipping motion and LOS is adopted: V a s c V d s c = sin γ sin θ a s c sin α a s c + cos γ sin θ a s c cos α a s c sin γ sin θ d s c sin α d s c + cos γ sin θ d s c cos α d s c − cos ⁡ ⁡ θ a s c − cos ⁡ ⁡ θ d s c V h o r i z o n t a l V U . (2)

The strike or azimuth of the horizontal displacement along the KBF is denoted by γ.

In order to assess the risk of local earthquake hazards and to better comprehend the Indian Plate’s northward movement, it is crucial to determine the slip rate presently posed by the KBF precisely. To calculate the interseismic deformation along the active fault, Savage and Prescott (1978) solved the following analytical formulation: V x = S π arctan x − x 0 d + a , (3) where the velocity of deformation caused by the fault is represented by V(x) and x denotes the distance from the fault. The fault centerline location is offset by x_o, the locking depth is marked by d, the current fault slip rate is denoted by S, and the offset between the profile and model is represented by a. Several assumptions are adopted to analyze the interseismic deformation of strike–slip faults; it is usually assumed that the fault has an indefinite length along the strike and is surrounded by a homogeneous elastic material. According to field research conducted by Yeats et al. (1984), it has been concluded that the Kalabagh Fault is a nearly vertical structure.

Previous studies have mentioned the creeping nature of the KBF. In order to estimate creeping in our studies, we used the following formulation (Segall, 2010; Hussain et al., 2016), in which we calculate the surface velocity V(x) at a particular offset x: V x = − S π arctan x + x c d 1 + C 1 π arctan x + x c d 2 − H x + x c + c , (4) where S is the slip, d₁ is the depth below which S occurs, the shallow creep rate is mentioned by C between the surface and depth d₂, H is the Heaviside function, and the static offset retained is c at another offset from the fault x_c.

The Coulomb stress model is deduced from Coulomb’s law of friction, which posits that the frictional force on a fault is directly proportional to the normal stress on the fault. By taking into account the deformation of an earthquake or external processes, the Coulomb stress model computes the change in the normal and shear stresses on a fault and determines whether such change is adequate to trigger slip on that particular fault (Lin and Stein, 2004; Borghi et al., 2016; Xiong et al., 2017; Hough and Bilham, 2018; Pope and Mooney, 2020). The Coulomb stress model is determined using the following equation: ∆ C F S = ∆ τ + μ ∆ σ .

In this equation, ∆τ is the shear stress and ∆σ is the normal stress, computed from known fault geometry, whereas μ represents the coefficient of friction.

All of the processing steps applied in InSAR data and the inversion scheme are summarized in the flowchart in Figure 2.

Figures 3, 4 show a map of the mean LOS velocities in the study area. The estimated velocities of both the ascending and descending tracks range from −20 mm/yr to + 20 mm/yr, with a maximum rate of approximately 20 mm/yr on the descending track in the Khisor Range and southeast of the KBF. Such large displacements in the LOS direction can be caused by several factors, including 1) the choice of reference point, which is relative to the measurement of InSAR observations; 2) topography-related error due to the mountainous region; 3) unwrapping errors; and 4) cumulative loop closure, which shortens the interferogram interval, which can result in signal uplifting due to the presence of phase biases. We estimate each error’s contribution and subtract it from the mean LOS velocities field to separate the signal from the errors. Collectively, noise levels on the KBF are relatively low in the range of −2.5 mm/yr to 2.5 mm/yr (see Figures 3, 4). Figure 5 depicts the corrected final mean velocity maps in the LOS direction for both ascending and descending tracks. These maps are georeferenced to the global coordinate system. The STD in the velocity field (Figures 5C, D) is estimated from the cumulative displacement employing the percentile bootstrap method (Efron and Tibshirani, 1986). The lower STD shows the best estimate of the displacement rate. Figure 6 shows the STD of the unwrapping phase in the spatial domain before and after the GACOS correction for each interferogram. Figure 6 represents the overall effect of the GACOS correction on the unwrapped interferograms. However, 384 out of 570 and 241 out of 360 interferograms were corrected after implementing the GACOS correction in the ascending and descending tracks, respectively. The GACOS correction improved the STD of interferograms on average by decreasing the STD from 3.18 cm to 2.68 and 2.91 cm to 2.49 for the ascending and descending tracks, respectively. Similarly, the STD decreased from 2.87 to 2.34 cm and 2.61 to 2.03 cm on median after the GACOS correction. The worsened STD interferograms after the GACOS correction may be related to the low accuracy of the GACOS turbulence data (Wang et al., 2019). The worst interferograms from the ascending track (Figure 6A) and descending track (Figure 6B) were removed before the estimation of the rate map. After removing atmospheric noise, the residual noise in the velocity measurements is significantly reduced. Finally, the velocity field obtained from both ascending and descending tracks was decomposed to components by selecting only those pixels that were available in both rate maps. Assuming a fixed strike of 165°, the LOS velocity is decomposed to horizontal and vertical components using equation 2. Figure 7 illustrates these components. The steep gradient of the LOS velocity across the Kalabagh Fault on both tracks is the primary feature observed in Figure 5. This gradient has a similar sign on both tracks, indicating right-lateral motion. The gradient reaches approximately 10 mm/year along the LOS, as illustrated in Figures 2–5.

Figure 8 shows cross-sections across the Kalabagh Fault, depicting the heterogeneous creeping along the fault’s strike. To avoid biases related to the topographic error, profiles are chosen at locations where the topography is almost identical on both sides of the fault. We plot a set of two 25 km-long AA/(on the northern segment of the KBF) and BB/(on the central segment of the KBF) and one 80 km-long CC/(on the southern segment of the KBF), cross-cutting the fault, as presented in Figure 8. Similarly, we plot the 25 km-long cross-section on the central section of Figure 8. Supplementary Figure S2 shows the profiles at the same fault segments without applying the topography correction. Each profile was kept consistently 2.5 km broad. The northern section of the KBF displays a significantly more pronounced horizontal velocity field discontinuity than the central section, which exhibits a higher fault-up velocity discontinuity. These observations are consistent with earlier InSAR and geological studies conducted by Chen and Khan (2010), McDougall and Hussain (1991), and McDougall and Khan (1990).

In order to analyze the creeping pattern along different segments of the fault, the LOS velocity components were inverted. This inversion enables us to quantify the creeping rate in the fault’s shallow region. We optimize the creep model parameters by considering a wide range of values for slip rates (ranging from 0 to 20 mm/yr), creep rates (ranging from 0 to 20 mm/yr), locking depths (ranging from 0 to 10 km), creep depths (ranging from 0 to 6 km), offset (ranging from −10–10 mm/yr), and offset from the fault (−20–20 km). These ranges were adopted based on previous studies conducted by Chen and Khan (2010), McDougall and Hussain (1991), and McDougall and Khan (1990). The observed signal may not always be perfectly aligned with the fault trace. In these cases, the offset to the location of the fault trace could be useful in order to get an accurate estimation of the creep rate and other parameters. The model parameters and their uncertainties are estimated through a Markov chain Monte Carlo approach. We use 1e6 iterations to define the posterior probability of the model using higher ranges of the model parameters. The maximum likelihood function by minimizing the weighted misfits is calculated between observations and models to rank the solutions. The first 20% of solutions are excluded as burn-in, whereas the remaining solutions are used to define the maximum a posteriori solution. The same approach was applied in previous studies on the estimation of creep or slip rates (Toda et al., 1998; Fattahi and Amelung, 2016). The velocity data on both sides of the fault were fitted with the best-fit 50 inversion models, as illustrated in Figure 9. We were unable to resolve the discontinuity over the KBF for the horizontal component of profile BB^/ because of negligible components of these directional creep components. However, we adopted the (Barnhart, 2017; Okada, 1992) forward model to fit the vertical displacement. This study is intended to assess the general mechanics of the KBF, including the variation of creep rate along the fault, locking depth, and creep depth. The horizontal velocity step in profile AA′ is the most noticeable, and it has the largest creep rate (4.8 ± 1.6 mm/year) within error bars. In contrast, profile CC′ consistently displays a lower creep rate than profile AA′. The models for profile BB′ did not show any substantial horizontal creep, but they predicted a vertical creep rate of approximately 10 mm/year (profile DD′). All investigated profiles on the fault have shown varying creep patterns on all three sections of the KBF. However, the results have shown good consistency with earlier studies carried out by Chen and Khan (2010) and suggest that the creeping occurs at depths ranging from 0.5 to 6 km (as shown in Figure 9). Our model did not indicate any significant slip rate along the Kalabagh Fault, which is in agreement with the region’s overall tectonics as the viscous decollement is present over Kohat–Potwar fold and thrust fault (Smith et al., 1994; Jaswal et al., 1997; Satyabala et al., 2012; Javed et al., 2022). Another noteworthy aspect of the vertical velocity field during the study period of 2015–2022 is the subsidence pattern, as shown in Figure 7. However, the maximum displacements in the horizontal component along Khisor Thrust and Salt Range Thrust (SRT) are not associated with the tectonics. This arises due to resolving the horizontal components to the geometry of the KBF. Additionally, we observed an asymmetrical pattern of horizontal and vertical velocities across the KBF. This could be due to the variation of fault dip with depth and heterogeneities of the elastic crust (Fialko et al., 2002; Fialko, 2006; Javed, 2017). Asymmetric velocity patterns have also been observed in other major active strike–slip faults around the world, including the San Andreas Fault (SAF) in the United States (Fialko et al., 2002; Fialko, 2006) and the Idrija Fault System in Italy (Javed, 2017). The modeling of asymmetric patterns, on the other hand, is beyond the scope of this work.

Finally, Coulomb 3.4 (Lin and Stein, 2004) was applied for computation of Coulomb failure stress changes (ΔCFS) to estimate the effect of creeping on local strike–slip and gentle dip (∼300) thrust faults. The strike of these other faults is taken similar to the KBF (i.e., strike ∼165), whereas the surface strike of the thrust geometry is adopted using rupture of local thrust faults (i.e., strike ∼ 250). Considering the creeping of the last 07 years, ΔCFS is computed more than 0.025 bar at 2.5 km and 5 km depths (see Figure 10). ΔCFS estimates positive on the northern and southeast regions of the KBF (see Figures 10C, D). However, the region of positive ΔCFS is limited to two small lobes on the right side of the KBF (Figures 10A, B).

A previous InSAR study on the interseismic deformation in the Kohat–Potwar region relies on a single pair of descending SAR images of ERS-1/ERS-2 across the Kalabagh Fault (Chen and Khan, 2010). In the context of InSAR processing, strain variations along the strike are often disregarded. However, this study demonstrates the potential for detecting such variations along faults using the extensive radar data archive provided by the Sentinel-1A satellite and applied an optimized processing approach that accounts for the reduction of atmospheric noise, similar to the Main Recent Fault in SW Iran (Watson et al., 2022). Despite this progress, limitations in our modeling persist, particularly regarding the simplified geometry and homogenous Earth crust.

Several studies have been conducted to explore the development of the KBF system and the northwestern Himalayan frontal thrust system’s progression over time (Yeats et al., 1984; Yeats and Hussain, 1987; McDougall and Khan, 1990; McDougall and Hussain, 1991; Smith et al., 1994; Bender and Raza, 1995; Khan and Glenn, 2006). The existence of salt has complicated the geometry and structural pattern of the frontal thrust system (Butler, 1987). The structural arrangement of the basement rock is the main factor influencing the Kalabagh Fault Zone’s development (McDougall and Khan, 1990). As suggested by numerous studies, the thrusting propagates forward on brittle detachment layers over the evaporitic decollement and brittle layers (Cotton and Koyi, 2000; Costa and Vendeville, 2002; Schreurs et al., 2002). However, a prior analysis of InSAR data demonstrates that the Chisal Algad stream experiences a strike–slip displacement rate of 6.2 mm/yr (Chen and Khan, 2010). In the last two Ma, McDougall and Khan (1990) calculated that a dislocation of approximately 7–10 mm/a has occurred near the Chisel Algad River, and a slow creeping of about ∼3.7 mm/a also persists along the Nammal Ridge (Chen and Khan, 2010). The findings of Chen and Khan (2009) indicate that the thrusting process has not always propagated consistently on a viscous décollement in the relation of the SRT and the Potwar Plateau. The SRT wedge, on the other hand, is creeping on a thick salt décollement. Based on their findings, the researchers concluded that current internal deformation outweighs the frontal deformation. Additionally, the outcomes of geomorphic and tectonic investigations support the various deformation patterns (Chen and Khan, 2009; Chen and Khan, 2010; Abbas et al., 2022). For instance, the western portion of the Salt Range is rising approximately 10 mm/yr (Abir et al., 2015). Additionally, Abbas et al. (2022) demonstrated the compression at the KBF’s central section. Our creeping measurement findings (i.e., 4.8 ± 1.6 mm/yr on horizontal velocity) are consistent with the work of Chen and Khan (2009) and support the idea that the northern region of the KBF is currently undergoing internal deformation. The absence of an up-dip component of the creep rate in that section supports the claim of lesser southward thrusting propagation by the Kohat–Surghar thrust wedge due to the more resistant nature of brittle detachment. Based on our InSAR findings, the center portion of the KBF has shown 10 mm/yr vertical creeping, confirming previous research’s observations of uplift (McDougall and Khan, 1990; McDougall and Hussain, 1991; Abbas et al., 2022). In addition to this, the KBF also showed the velocity-weakening region below 4 km depth on both the northern and southern sections of the fault. The occurrence of earthquakes, as reported by Abbas et al. (2022), could be the result of partial creeping at depth. However, estimation of the slip rate and looking depth is poorly constrained. In comparison to the center and northern sections of the KBF, the southern section of the KBF demonstrates a relatively lower creep rate. The predicted horizontally creeping rate of ∼4.2 ± 1.3 mm/yr is in good accordance with the earlier observation for this section (Chen and Khan, 2010). Several earlier studies concluded that the Potwar–Kohat plateau is accommodating 5–10 mm/yr deformation, but these estimations are based on geological studies, sparse GPS estimation, and limited ERS-InSAR data coverage (Wells, 1984; Yeats et al., 1984; Yeats and Hussain, 1987; McDougall and Khan, 1990; McDougall and Hussain, 1991). A recent study by Abass et al. (2022) showed the shallow duplex structure of Kohat accommodating 5 mm/yr. Similarly, Jouanne et al. (2020) used modeling of interseismic coupling of MHT beneath Northern Pakistan and illustrated that the Potwar Plateau is creeping at a rate of approximately 6 mm/yr, which is approximately 20% of the total deformation rate (i.e., 3 cm/yr) between the Indian and Eurasian plates. Our estimates also show that the creeping rate along the Kalabagh Fault is approximately 5 mm/yr, so creeping is the dominant phenomenon beneath the Kohat–Potwar plateau. Last, the asymmetrical distribution of horizontal and vertical velocities over the KBF is connected to the lateral variation in rigidity caused by the presence of a significant salt layer to the east of the KBF and complex fault geometry at depth. The impact of various parameters is highlighted by several authors in their studies (Efron and Tibshirani, 1986; Le Pichon et al., 2005; Lindsey and Fialko, 2013).

Various models have been utilized in the past to comprehend the interaction between the Indian and Eurasian plates. The crustal thickening and shortening are accounted for by the thin viscous model. On the other hand, crustal shortening taking place along large-scale strike–slip boundaries is addressed by extrusion tectonics models (Molnar and Tapponnier, 1977; England and Houseman, 1989; Harrison et al., 1992; Houseman and England, 1993). However, the precise kinematics and geometry of extrusions are not yet fully understood, as many structures, such as transform faults that facilitate the extrusions, are not identified or accurately mapped. Structural geometry disparities are frequent in the Himalayan orogeny, as noted by Chen and Khan (2010), and these variations significantly impact the deformation pattern, outcrop generations and movements, and sedimentary history of the orogeny, as discussed by Yin and Harrison (2000).

Finally, the interaction between the KBF and nearby local strike–slip and thrust faults are evaluated using the Coulomb stress changes. In the past, numerous studies demonstrated that ΔCFS >0.01 MPa is enough to trigger seismicity (Toda et al., 1998). To the east of the KBF, the decollement thrust–fold assemblages are the prominent features of the Potwar Plateau. However, the Indian Plate slips aseismically beneath the Kohat–Potwar thrust and fold belt (Satyabala et al., 2012), so large earthquakes are unlikely in this region due to viscous decollement (Satyabala et al., 2012). However, this region is characterized by a partial stick–slip phenomenon occurring on one or more locked portions of fault on the ramp for longer temporal gaps that are connected to the decollement. The increase in stress loading on adjacent fault segments (i.e., northern and southeast of the KBF) can have serious seismic implication.

Sentinel-1A ascending and descending data for the period of 2015–2022 are used to compute the rate map of the KBF, one of the distinctive features of the Indian and Eurasian plate collision. Our inversion showed a variation of creep rates at different segments of the fault. Horizontally and vertically decomposed components of velocity at three different cross profiles of the KBF zone are presented in this study. The fault creeping rates corresponding to the horizontal component at the northern segment and southern end of the fault showed creep rates of ∼ 4.8 ± 1.6 and ∼4.2 ± 1.3 mm/year, respectively, whereas the central section of the fault showed negligible horizontal component of velocity but a vertical component of ∼10 mm/year. The up-dip creep rate illustrates the rapid vertical uplift at the central fault segment of the KBF and is in good accordance with the thrusting and folding. Moreover, our model did not rule out an interseismic slip rate at depth, which could be linked to earthquakes in the past. Our findings of the heterogenous creep rate on the KBF will assist the understanding of the mechanics of Indian–Eurasian collision. This study will enable us to better address the seismic hazard assessment in the region.