Neotectonics of the Western Suleiman Fold Belt, Pakistan: Evidence for Bookshelf Faulting

: The Suleiman Fold-Thrust Belt represents an active deformational front at the western margin of the Indian plate and has been a locus of major earthquakes. This study focuses on the western part of the Suleiman Fold-Thrust Belt that comprises two parallel NW–SE oriented faults: Harnai Fault and Karahi Fault. These faults have known thrust components; however, there remains uncertainty about the lateral component of motion. This work presents the new observation of surface deformation using the Small Baseline Subset (SBAS), Interferometric Synthetic Aperture Radar (InSAR) technique on Sentinel-1A datasets to decompose displacement into the vertical and horizontal components employing ascending and descending track geometries. The subsurface structural geometry of this area was assessed using 2D seismic and well data. In addition, geomorphic indices were calculated to assess the relative tectonic activity of the area. InSAR results show that the Karahi Fault has a ~15 mm right-lateral movement for descending and ~10 mm/for ascending path geometries. The Harnai Fault does not show any lateral movement. Seismic data are in agreement with the InSAR results suggesting that the Harnai Fault is a blind thrust. This work indicates that the block between these two faults displays a clockwise rotation that creates the “bookshelf model”. The well data show that the ﬁrst 414 m represent the Eocene Ghazij Shale. The Paleocene Dunghan Limestone follows with a thickness of 254 m. Other stratigraphic intervals drilled include Moro, Mughal Kot, Parh, Upper Goru Formations that are Late Cretaceous in age, and the Lower Goru and Sembar Formation of the early Cretaceous. The deepest stratigraphic horizon drilled in Ziarat-01 well was the Middle Jurassic Chiltan Formation at a depth of 833 m.


Introduction
Suleiman Fold-Thrust belt (SFTB) is a broad, lobate, and asymmetric active foreland deformational feature, more than 400 km wide at the western margin of the Indian plate ( Figure 1). It is formed as a result of oblique convergence between the Indian and Eurasian plates and provides a unique opportunity to study oblique margin processes [1,2]. The regional structural trend in the Central and Northwestern Himalayas is east-west and switches to the north-south in the Western Himalayas. The Indian plate is moving at a rate of~35 mm/year and~38 mm/year in the northwestern and northeastern margins, respectively [3]. Most of the present-day compressional tectonics are associated with the frontal thrusts and are well studied in the Central Himalayas in Nepal [4,5], in Northwestern Himalayas in India [6], and across the Salt Range Thrust in Pakistan [7,8]. However, the transpressional tectonics at the western margin of the Indian plate is poorly understood, although important in understanding active plate margin processes and calculating interseismic strain accumulation for an earthquake-prone area.
The continental-scale left-lateral strike-slip Chaman Fault marks the western boundary of the SFTB and the Indian plate. Suleiman Ranges (SR) merges with the Waziristan ophiolite belt with a pronounced structural reentrant known as Tank Reentrant in the north. In the south, it is bounded by the Sibi reentrant [9]. Both of these reentrants are formed  shows the location of the Sentinel 1 ascending and descending track geometry used in InSAR processing; the blue frame shows the location of the seismic lines. The red and cyan rectangular areas represent Vf index locations. The yellow circles illustrate the earthquake of magnitude 5 or above. Solid arrows indicate the GNSS velocity vector obtained from Szeliga et al. [13], showing south-southeastern movement across Harnai and Karahi Fault. All stations movement are relative to the stable Indian plate as defined in Altamimi et al. [28] and are plotted using a Web-Mercator projection. The GNSS station SANJ across the Karahi Fault has a movement rate of~20 mm/year, while station HRNI across the Harnai Fault has a movement rate of~10 mm/year [14]. (B) Three-dimensional perspective view of the 30 m high-resolution digital elevation model (SRTM) showing the location of the seismic lines with topography. . Earthquakes with a magnitude greater than 5 and their focal mechanism solution derived from USGS are plotted to show their distribution across the SFTB [26]. A high density of earthquakes is plotted in the western Suleiman near to Quetta Syntaxis/Sibi reentrant. The foreland propagation of the deformation front is also prominent. SR, Suleiman Ranges.

Geotectonic Setting of the Study Area
The SFTB is sub-divided into the eastern, central, and western parts. A prominent left-lateral strike-slip Kingri fault divides the eastern and central parts of the SFTB with a sharp contrast in their structural trend ( Figure 3). The eastern part (Suleiman Ranges) trends north-south, parallel to the direction of Indian plate convergence, while the central part (Suleiman Lobe) trends east-west, perpendicular to the Indian plate convergence [29]. The structural geometry is thin-skinned with basal detachment interpreted either in Permo-Triassic strata [29] or in Precambrian-Cambrian evaporates [30]. Folds describe the dominant structural geometry in Suleiman Ranges and the Suleiman Lobe and is devoid of thrusting. However, western SFTB extends from the Muslim Bagh ophiolites to the Sibi Trough and shows dominant fold and thrust geometry [30]. Continental-scale left-lateral strike-slip Chaman Fault marks the western boundary of the Indian plate and the SFTB, while Punjab Foredeep marks the eastern limit of Suleiman Ranges (SR). Present-day seismicity is well-distributed along the SR's deformation front and in the western part [12].
The structural elements of the western SFTB include two parallel, NW-SE striking, transpressive faults, i.e., Harnai and Karahi faults ( Figure 3). The surface trace of the Karahi fault is approximately 92 km and it offsets folds with a right-lateral sense of slip. In comparison, the surface trace of the Harnai Fault is approximately 76 km. At the surface, both of these faults do not merge with the Ghazaband Fault, which is the splay of the Chaman Fault; however, their subsurface connectivity can be interpreted with better Remote Sens. 2021, 13, 3593 6 of 27 seismic data. The surface stratigraphy comprises molasses deposits of the Siwalik Group, a mechanically weak unit of Ghazij shales, and a component limestone unit of the Kirthar and Dunghan Formation ( Figure 4) [26,31]. The drainage pattern shows prominent deflection across the nose of the fold, with some streams flowing parallel to the faults.

InSAR
InSAR is a powerful technique used to monitor surface deformation with centimeterto millimeter-level accuracy [20]. The displacement is observed relative to the radar line of sight (LOS) by calculating the phase offset between the incident and return waves. The phase difference between the two images over the same area with different time interval produces interferograms which represent deformation during that time frame [32]. The satellite with a shorter wavelength (C Band, Sentinel 1) is more sensitive to displacement than the satellite with a larger wavelength (L Band, ALOS-PALSAR). In this study, 45 ascending (path 144, frame 193) and 45 descending (path 78, frame 491) Sentinel-1A Terrain Observation with Progressive Scans SAR (TOPSAR) images (C-band) that cover the area between Karahi and Harnai faults with a twelve-day revisit period (Figure 3) were used (Table 1). Temporal coverage for both ascending and descending track geometry ranges from 20 September 2018 to 1 March 2020, and from 24 September 2018 to 17 March 2020, respectively. Both ascending and descending datasets were processed with the help of the SARscape module of the ENVI software from Harris Geospatial Solutions. The processing steps are summarized in Figure 5. InSAR SBAS technique forms interferograms with the short spatial and temporal baseline to minimize decorrelation (slow varying phase pixels) and has been successful in monitoring deformation in non-urban environments [33]. A maximum temporal baseline of 150 days was selected since there is no vegetation in the study area to cause surface temporal decorrelation (Table 1). This temporal baseline together with the perpendicular baseline threshold of 45% improved the coherence of the stacked image pairs (interferograms) (Figure 6). High-resolution (12 m) Tandem-X DEM were obtained from the German Aerospace Agency to account for topographic corrections. To remove the ramp and phase constant during reflattening and refinement, reference ground control points (GCP) were selected away from actively deforming areas on highly coherent pixels (coherence > 0.7) with no phase jumps and were uniformly distributed throughout the study area. In addition, the filtered interferograms with topographic fringes were removed to provide a powerful basis to select GCP's away from displacement fringes. Following this, the two inversion models, i.e., First inversion and Second inversion, were used to estimate the rate of displacement. In the first step, a linear displacement model was used to obtain the first velocity estimate and retrieve residual topography, while the second step removed atmospheric fringes by filtering. Finally, displacement time series were generated and geocoded from slant range to ground range (UTM zone 42N) using the reference DEM ( Figure 5). . SAR interferograms showing spatial-temporal baseline distribution, where each acquisition is represented by a diamond associated with an ID number; the yellow diamond represents reference acquisition and green diamond represents the secondary time position plot of all interferograms generated by S1 data showing a perpendicular baseline between different image's acquisition according to their relative position. Table 1. Sentinel-1A data specification. The product type is L1 single look complex (SLC) with ascending and descending geometry and interferometric wide swath beam mode (IW). The following acronyms were used in the table: polarization (Pol), incident (Inc), azimuth (Az), range (Rg), wavelength (WL), maximum perpendicular baseline (Bp) and maximum temporal baseline (Bt).
The acquired seismic data were preprocessed by DGPC up to time-migrated images. The preprocessing includes demultiplexing, geometry, trace editing, amplitude recovery, despike, f-k filter, surface consistency deconvolution, refraction statics, first velocity analysis, first residual statics predictive deconvolution, second velocity analysis, second residual statics, third velocity analysis, final vertical moveout analysis, final update RMS velocity, post-stack migration using stacking velocity, display filter, and gain. All the seismic and well data were imported into the IHS Markit Kingdom software (Academic Version 2018) for structural and stratigraphic interpretation. Seismic attributes analysis (Pseudo-relief, Average Energy, and Trace Envelope) was carried out to assist in the interpretation of the seismic lines.

Tectonic Geomorphology
Two different geomorphological indices, i.e., Hypsometric Integral (HI) and Valley Floor Width to Height Ratio (Vf) were used to assess the relative tectonic activity in the study area.
A hypsometric integral is usually determined by calculating the cumulative height and the cumulative area below that height for individual watersheds, followed by selecting the area under that curve to measure the hypsometric integral [34]. HI is used to assess the stages of landscape evolution from youth through mature to old, and gives an idea about the tectonic evolution. Strahler [34] proposed that different drainage basins represented different shapes of hypsometric curves and classified these drainage basins in terms of their geomorphic evolution ( Figure 7). In general, lower values describe older eroded landscapes in the late stages of geomorphic evolution. However, higher values of the HI usually indicate that less of the uplands have been eroded and may indicate a younger landscape, perhaps generated by active tectonic processes ( Figure 7). This geomorphic index defines the distribution of elevations across the study area. Relative height against relative area graph (normalized to 1) can show three geometries. Convex hypsometric curves represent relatively "young" weakly eroded regions, S-shaped curves describe moderately eroded regions, and concave curves illustrate relatively "old" highly eroded regions [35][36][37] (Figure 7).
In this study, individual watersheds were extracted using 12 m high-resolution DEM. The area of the drainage basins was calculated using the spatial analyst tool in ArcMap 10.7. The hypsometric integral was computed using the approach of elevation-relief ratio introduced by Pike and Wilson [38]. The ratio of Valley Floor Width to Height (Vf) was introduced by Bull [39] to determine the differences between V-shaped and U-shaped valleys ( Figure 8). Vf is described as: where Vfw is the valley floor width, Esc is the elevation on the valley floor, Erd and Eld is the elevation on the right and left side of the valley, respectively. According to Bull [39], the high Vf values indicate a low uplift rate that means dominant erosional processes and relatively tectonic stagnation, while low Vf values indicate deep erosion force with tectonic uplift. In other words, this index differentiates between broad-floor canyons, with relatively high values of Vf, U-shaped, and V-shaped valleys with relatively low values of Vf (Table 2) [39][40][41].

Present-Day Deformation
Within the study area, there are two main faults named the Karahi and Harnai Faults. Movement along these faults compartmentalized deformation into three areas/spots ( Figure 9A). Deformation spot 1 gives information about the movement of the northern portion of the Karahi Fault. In contrast, deformation spot 3 provides information on the southern portion of the Harnai Fault movement. Deformation spot 2 incorporates the movement of both faults and displays how the movement affected this interior block.
Results from the descending satellite path, illustrated in Figure 9A, display positive values (red color), signifying the surface movement towards the satellite, which can be interpreted as uplift, eastward movement, or both. On the other hand, negative values (blue color) represent the surface moving away from the satellite and can be interpreted as subsidence, westward movement, or both ( Figure 9A). Additionally, the radar can only measure displacement in the radar's line-of-sight (i.e., east-west and uplift/subsidence) however, it is insensitive to north-south displacement because of its geometry.
The descending track time series shows the ground LOS displacement at five different locations within the study area from September 2018 to March 2020. Active structures with an offset by Harnai and Karahi Faults were selected for their displacement time series. Five Regions of Interest (ROI) were chosen for each location from the velocity map, and the displacement data for the ROI's were plotted on the y-axis against the time on the x-axis in the velocity-time graph ( Figure 9A). It can be seen that the areas covered by ROI-1 and ROI-2 on the velocity map are red, showing an average of 10 to 15 mm of movement towards the satellite. The LOS movement in these areas is more pronounced from August 2019 to September 2019, reaching a value of 20 mm. Areas covered by ROI-3, ROI-4, and ROI-5 show movement away from the satellite ( Figure 9A).    Figure 9B shows the ascending track geometry in which the satellite moves from south to north. In this image, the northernmost section of the Karahi Fault displays negative velocity values (blue color), indicating that the area is moving away from the satellite LOS. Movement of the Harnai Fault was not detected in the ascending velocity map. In contrast to the descending results, positive values (red color) shown in the ascending results can be interpreted as uplift, westward movement or both while negative values (blue color) can be interpreted as subsidence, eastward movement or both ( Figure 9B).
The ground LOS displacement was analyzed at three different locations in the study area using the ascending velocity map ( Figure 9B). Three ROIs were selected for each deformation spot (1, 2, and 3) from the velocity map and time series were generated. It can be seen that the area of ROI-1 and the surroundings have been experiencing a slow LOS movement initially for some time but then subsided to the previous position. However, the areas of ROI-2 and ROI-3 have been experiencing subsidence or movement away from the satellite, which is noticeably very high for ROI-3 during recent times, reaching up to 15 mm ( Figure 9B).
To get the final horizontal and vertical velocity maps, we performed the displacement decomposition in ENVI SARscape with the help of ascending and descending track geometry. Since the satellite is sensitive to both vertical and horizontal components of motion, the tool decomposes the displacement and calculates the angles in the satellite LOS. The two angles calculated are described as Azimuth Line of Sight (ALOS) and Inclination Line of Sight (ILOS). The horizontal component of motion is described as east-west movement ( Figure 9C). The eastward movement is depicted by the positive values (red color), while the westward movement is depicted by the negative values (blue color). According to this information, deformation at spot 1 suggests dominance towards the east direction, as shown by the positive values or red color, and this is also supported by Figure 9C. The plot of ROI-1 illustrates around 20 mm eastward movement, while the plot of ROI-2 does not show any displacement ( Figure 9C). The Harnai Fault area displays a green color, which indicates a stable velocity component ( Figure 9C) and clear evidence for no movement along the Harnai Fault. It is evident that the area of ROI-5 has not shown any displacement; however, the area of ROI-6 shows slight subsidence of approximately 5 mm ( Figure 9C). Figure 9D illustrates the vertical component of movement as either subsidence or uplift. The uplift is characterized by positive values (red color), while subsidence is characterized by negative values (blue color). According to Figure 9D, for deformation spot 1, there is a slight uplift on the northwest portion shown as red. This uplift is shown on the plot of ROI-1 with 10 mm movement rate ( Figure 9D). The plot of ROI-2 shows stable movement throughout the observation time ( Figure 9D). These two ROI-1 and ROI-2 show the northern and southern parts of the Karahi Fault, respectively. This could provide information about the thrust component of the Karahi Fault. Additionally, the eastern part of deformation spot 1 shows a subsidence deformation zone that could explain the boundary of Fault A. The areas of ROI-5 and ROI-6 do not show any vertical displacement ( Figure 9D). Figure 2B shows the location of 2D seismic lines and one well in the study area. At the surface, unconsolidated alluvium and sediments of the Ghazij Formation are exposed. In order to understand subsurface stratigraphy and details of faulting, 10 seismic lines were interpreted, but only these three fault lines were illustrated in this paper (Figures 10-12). Three of these lines show the Harnai fault and 3D and 2D surfaces generated using these 10 seismic lines and well data. Figure 10 shows the 2D seismic interpretation of line ZRT04-02. For stratigraphic interpretation, Ziarat-01 well was tied to seismic line ZRT04-02 using the well to seismic tie function in Kingdom. The interpreted horizons were then transferred to other seismic lines for understanding structural complexity. Additionally, time to depth conversion was made by using well-log velocity information. Ziarat-01 well was drilled to the depth of 1050 m 1050 meters ( Figure 10A). The well data show that the first 414 meters represent the Eocene Ghazij Shale. The Paleocene Dunghan Limestone follows with a thickness of 254 meters. Other stratigraphic intervals drilled include Moro, Mughal Kot, Parh, Upper Goru Formations that are Late Cretaceous in age, and the Lower Goru and Sembar Formation of the early Cretaceous. The deepest stratigraphic horizon drilled in Ziarat-01 well was the Middle Jurassic Chiltan Formation at a depth of 833 meters.

Geophysical and Subsurface Investigation
The Moro Formation consists of 27 m-thick sandstone, while Mughal Kot Formation comprises interbedded shale and sandstone with some limestone intervals. The Parh Formation consists of White Limestone. The Goru Formation contains an assortment of lithologies, including marl, shale, limestone, and quartzose sandstone. The succeeding formation, the Sembar Formation, consists mainly of shale. The total thickness of these formations is 165 meters (from Moro to Sembar Formation). reflectors on the seismic section ( Figure 10A-D) that provide reliable interpretation f seismic lines and also make a straightforward correlation between each seismic line. Ho ever, other formation tops cannot be identified except for these three formations due the low vertical resolution of these lines. Therefore, Mughal Kot, Parh, Lower and Upp Goru, and Sembar formations tops have not been illustrated in the seismic section (Figu 10).  The 3D and 2D perspectives of the tops of the Dungan, Moro, Chiltan formations, newly identified two blind thrusts, and Harnai Fault are displayed in Figures 12 and 13. An anticlinal structure can be seen at the location of the Ziarat-01 well. This elliptical anticline structure is east-west oriented, and its depth increases towards the east direction. This prospect has been drilled with the Ziarat Gas Field discovery in 2005 [27].  At the top of Dungan Limestone, the Moro and the Chiltan formations show strong reflectors on the seismic section ( Figure 10A-D) that provide reliable interpretation for seismic lines and also make a straightforward correlation between each seismic line. However, other formation tops cannot be identified except for these three formations due to the low vertical resolution of these lines. Therefore, Mughal Kot, Parh, Lower and Upper Goru, and Sembar formations tops have not been illustrated in the seismic section ( Figure 10).
The 3D and 2D perspectives of the tops of the Dungan, Moro, Chiltan formations, newly identified two blind thrusts, and Harnai Fault are displayed in Figures 12 and 13. An anticlinal structure can be seen at the location of the Ziarat-01 well. This elliptical anticline structure is east-west oriented, and its depth increases towards the east direction. This prospect has been drilled with the Ziarat Gas Field discovery in 2005 [27].
InSAR results presented in Figure 9 demonstrate

Tectonic Geomorphology
To determine the role of recent tectonic activities on the drainage basins, hypsometric integral and valley floor width to height indices were calculated for all the watersheds across the Harnai and Karahi faults. Fourteen watersheds for the Harnai fault and 10 water-sheds for the Karahi fault were selected to illustrate hypsometric curves (Figures 14 and 15). While the average hypsometric integral ratio of watersheds for the Karahi Fault is 0.469 (Table 3), the average hypsometric integral for the Harnai Fault is 0.392 (Table 4). Analysis of the hypsometric integral is based upon digital elevation models and utilization of all watershed results shown in Tables 3 and 4. Overall, higher values of the hypsometric integral are illustrated by a convex shape curve indicating a young stage of drainage evolution. These values are generally > 0.5. Intermediate values tend to be both concave and convex, displayed by an S-shaped curve with values ranging between 0.4 and 0.5, which is considered a mature basin. Lastly, lower values < 0.4 tend to have concave-shaped curves ( Figure 7B).
In general, lower values describe older eroded landscapes in the late stages of geomorphic evolution. However, higher values of the index usually indicate that less of the uplands has been eroded and may indicate a younger landscape, perhaps generated by active tectonic processes. Therefore, watersheds related to the Karahi fault, which demonstrated a majority of "S"-shaped curves contain both young and mature basins. However, watersheds relating to Harnai Fault demonstrate a majority of concave curve shapes, indicating a mature basin ( Figure 16).
A total of 26 watersheds were extracted and analyzed for their Valley Floor Width-to-Height Ratio. Ten watersheds were evaluated for their Vf value for the Harnai Fault and show an average value of 0.83, indicating a relatively moderate level of tectonic activity. For the Karahi Fault, 14 watersheds were analyzed and they show an average Vf value of 0.32, which is quite low in comparison to the Harnai Fault and represents high tectonic activity in this area.

Discussion
The western part of the Suleiman Fold Belt is a structurally complex area represented by a series of thrust and strike-slip faults with varying geometry and kinematics. Previous researchers have made a few efforts in an attempt to describe the sense of movement of the parallel, NW-SE trending Karahi and Harnai faults, and interpret the style of defor-

Discussion
The western part of the Suleiman Fold Belt is a structurally complex area represented by a series of thrust and strike-slip faults with varying geometry and kinematics. Previous Our InSAR results, seismic interpretation, geomorphological indices, and published GNSS velocities suggest that shearing has played a significant role in the deformation of the western Suleiman Fold Belt (Figure 9). The shear zone has undergone simple shear deformation (Figure 17). This may be related to the oblique convergence of the Indian plate moving at a direction of N05E relative to the Eurasian plate [42,43]. Oblique convergent boundaries have resulted in a 3D distribution of deformation (simultaneous interaction between contractional and strike-slip motion), observed in the shear zone of the western Suleiman fold belt [44,45]. Our results agree with Teyssier's strike-slip partitioned transpression model shown in Figure 17 [44]. The model is effective in understanding the relationship between plate motion, instantaneous strain axes, and degree of strike-slip partitioning. The shear zone is strongly affected by the movements along the Karahi and Harnai faults. While the Karahi Fault displays a sense of strike-slip and thrust motions on In-SAR results (Figure 9), the Harnai Fault only shows thrust motion on seismic profiles (Figures 10-12).
The compressional forces perpendicular to the Harnai and Karahi faults seem to dominate when compared to lateral motions along these faults. As a result of the shearing of this zone, possible left-lateral strike-slip movement (Fault A and Fault B) compartmentalized this zone into distinct blocks such as Block A and Block B. If the Karahi Fault is interpreted as a right-lateral strike-slip fault, then Harnai Fault would be considered a solid block, and the shear zone of Fault A and Fault B can be interpreted as left-lateral strike-slip faults. These interpretations support the "bookshelf model" with the interpretation of shearing in this region having a clockwise rotation ( Figure 17).
The bookshelf model, proposed initially by Mandl [46], is the behavior of parallel faults that relate to a horizontal stack of books being displaced and consequently named 'bookshelf faulting' [46][47][48]. Mandl [46] described two basic modes of mechanism: the "dilational" mode and the "domino" mode. While he proposed the dilational mode as two parallel faults that tend to open up and possibly become paths for fluid migration, the domino mode is described as reclining faults that remain closed. In both scenarios, the formation and rotation of the faults may occur in one tectonic event [46]. The bookshelf model is generally applied to describe the antithetic accommodation of deformation affected by direct shear.
In conforming to the model, the southern block of the Harnai Fault is assumed to be fix and rigid. InSAR and seismic results support this assumption. The average rightlateral movement of the Karahi Fault is at a rate of around 15 mm/yr ( Figure 9B), while the Harnai Fault does not show strike-slip movement on InSAR datasets ( Figure 9A-C). Since the Karahi Fault has higher displacement than the Harnai Fault, the net movement should result in a fixed southern block. The interpretation of this study is similar to the rotational block model suggested by Szeliga [13] and Huang et al. [18]. According to this model, the shear zone of the western Suleiman Fold Belt is compartmentalize into blocks that have rotated in a clockwise rotation due to right-lateral shear [13]. One distinct characteristic of this model is the presence of parallel strike-slip faults (Fault 1 and Fault 2) which resemble a horizontal stack of books being displaced and hence the name "bookshelf faulting" (Figure 18). Anticlinal structures are visible on the surface, indicating an uplift on the InSAR dataset between Karahi and Harnai faults ( Figure 18D). Another embedded anticlinal structure was observed on the seismic section in the northern portion of the study area. These actively deforming anticlines mostly have a northeast-southwest direction, suggesting the "bookshelf model" with clockwise rotation.
Seismic data allowed us to understand the geometry of Harnai Fault in the Ziarat Block. Jadoon et al. [27] recognize the structures in the Ziarat Block as an in-sequence passive roof duplex and out of sequence deformation with the segmentation of the duplex structures. While Jadoon et al. [27] show Fault 1 and Fault 2 as blind thrust faults, they did not depict the Harnai Fault in their interpretation. They also describe Fault 2 as the main detachment fault. However, in this work, we described all these faults (Harnai, Fault 1, and Fault 2) as blind thrust faults with no surface trace or disruption. The area between the Harnai Fault and Fault 1 illustrates Block A and the area between Fault 1 and Fault 2 is interpreted as Block B. These block structures and existing faults may also indicate the presence of the "bookshelf model" (Figure 18).
Hypsometric integral results in this study, following Perez's [37] and El Hamdouni's [40] classifications, indicate that the Harnai Fault is comparatively less active or is in a more mature stage than the Karahi Fault, which is tectonically more active (Tables 3 and 4). In addition, the northern part of the Karahi Fault can be interpreted as having relatively high tectonic activity due to the lower Vf values. In contrast, the southern part of the Harnai Fault is interpreted as having lower tectonic activity because of higher Vf values (Tables 3 and 4). Overall, both the Hypsometric integral and Valley Floor Width to Height Ratio results support the InSAR findings, suggesting higher displacement along the Karahi Fault than the Harnai Fault.

Conclusions
This study integrates seismic reflection, InSAR SBAS analysis and geomorphological analysis to picture subsurface geometry and map surface deformation in the western part of the Suleiman Fold-Thrust Belt along the Harnai and Karahi Faults. Our findings are summarized below;

•
The Karahi Fault is a right-lateral strike-slip fault, confirmed by the ascending and descending path geometry of Sentinel-1 datasets. This conclusion is also supported by the GNSS dataset. This right-lateral movement of the Karahi Fault could potentially be the reason for the clockwise movement. • An average right-lateral movement of~15 mm/yr is calculated for Karahi Fault for descending and~10 mm/yr for ascending track geometries, respectively. Also, the SANJ GNSS station close to the Karahi Fault shows the 20 mm/yr rate of displacement in the southeast direction.

•
No displacement is observed on the velocity map for the Harnai Fault and is interpreted as an aseismic.

•
The 2D seismic interpretation shows the thrust component of the Harnai Fault. However, the lateral component of motion was not identified on the 2D seismic for both faults.

•
Fault A and B observed on the seismic sections might be evidence of the compartmentalization between the Karahi and Harnai faults, and we interpreted this to be best explained by a "bookshelf faulting" model.