Lineament mapping for a part of the Central Sulaiman Fold–Thrust Belt (SFTB), Pakistan

This study aims to analyze the lineaments using the field data and a remote sensing approach, to describe their relationship with the folds, faults, and regional tectonic stress of the central Sulaiman Fold-Thrust Belt. Joint data, from nine anticlines, has been collected using the scanline method and classified into three sets (JS1, JS2, and JS3) based on their geometrical relationships. Lineaments extracted from the 30 m digital elevation model have been classified subsequently into three lineament sets (LS1, LS2, and LS3) based on the azimuths from the corresponding joint sets. A very high correlation coefficient (rs = 0.97), between the azimuths of the field joints and the remotely sensed lineaments, has been observed which validates that the lineaments are the regional representation of the local field joints. The geometrical relationship of the lineament sets with the fold hinges indicates that the older LS1 and LS2 are strongly related to the regional folding episode, while the younger LS3 is a result of local shears. The chronological interpretation of the deformational events responsible for the lineament sets is constrained by the presence of the Kingri Fault, which induces a strike-slip component within the study area. Furthermore, the controls on the joint and lineament sets, established using multivariate statistics to decipher the effects of lithological and structural contrasts on the lineament density, reveal that an increase in the lineament density can be attributed to the competence and thickness of the rock units as well as the variable local stresses across the different folds. Based on the orientations of these lineament sets, the cumulative direction of the compressive event in the NW-SE direction (310–320) coincides with the regional stress direction of the SFTB.


Introduction
Lineaments are the significant lines on the earth's surface that reveal the hidden architecture of the basement (Hobbs 1904). The definition has been revised by Sabins (1996) where he defined lineaments as parts of a linear or curvilinear feature with a distinct relationship to a surface that may be an expression of a fault or other linear features symbolizing discontinuity. Different lineaments (e.g., faults, fractures, rock contacts, shear zones, dykes, etc.) can be classified based on their geometry, nature, and genesis and have been studied in the various disciplines of geosciences. Lineaments have been used in the petroleum industry (Rajabi et al. 2010), hydrogeological studies (Krishnamurthy et al. 1996;Tam et al. 2004;Galanos and Rokos 2006;Chandrasekhar et al. 2011), mineral exploration (Rameshchandra Phani 2014), geological structural mapping (Drury and Holt 1980;Nagal 2014), and seismology and related hazards (Stefouli et al. 1996;Ramasamy 2006).
In contrast to the concept, the methods and respective approaches to acquire, process, and present lineaments have evolved (Florinsky 2016), and include both the manual and digital techniques. Due to the ease in availability of the satellite data, digital techniques have become more frequent in the lineament extraction analysis. Some of these techniques are the Hough Transform (Wang and Howarth 1990), Segment Tracing Algorithm (Koike et al. 1995), Haar Transform (Majumdar and Bhattacharya 1988), and the Lineament Extraction algorithm of PCI Geomatica software (Hubbard and Thompson 2012;Thannoun 2013). The use of the automated lineament extraction coupled with a geographical information system (GIS) environment can provide an effective and reliable way to quantitatively analyze lineaments and generate the corresponding geospatial maps. However, the automated extracted lineaments must be correlated with the field data to enable the lineaments to be regarded as representative for the structural manifestation of an area (Mogaji et al. 2011).
The present study aims to analyze the lineaments within the Sulaiman Fold-Thrust Belt (SFTB) located in the Baluchistan province, Pakistan (Fig. 1a), between north latitudes 30°19′ 00″-30°43′00″ and east longitudes 69°03′00″-69°40′00″. Baluchistan is one of the least developed provinces and thus least explored in terms of natural resources. In terms of the availability of the geological data, detailed structural studies have been carried out in the SFTB Jadoon et al. 1994Jadoon et al. , 2019Siddiqui and Jadoon 2012). However, studies relating the field joints with the remotely sensed lineaments have not been reported yet. Furthermore, given the current hydrocarbon exploration activities in the study area, the study of lineaments can also serve as a useful input toward refining the already existing surface geological maps and can aid towards understanding the subsurface structural setup. The main aims of this study are (1) to evaluate the lineament architecture of the study area using joints data acquired from the field observations and the lineament data acquired from the remote sensing and (2) to statistically analyze the extracted lineaments with the field joints to evaluate their significance toward the regional geological makeup of the SFTB.

Geological setting
In the northwest and west Pakistan, the oblique collision of the Indian Plate with the Eurasian Plate in a transpressional zone allowed the development of the different fold-thrust belts including the Sulaiman Fold-Thrust Belt (SFTB) (Farah and DeJong 1979;Lawrence et al. 1981;Farah et al. 1984;Jadoon et al. 1992Jadoon et al. , 1994Jadoon and Khurshid 1996). The western boundary of the SFTB is marked by the Chaman Fault (Jadoon and Khurshid 1996). The Chaman Fault is a left lateral shear with a small convergence component (Quittmeyer and Jacob 1979;Prevot et al. 1980), and accommodates the NW relative plate motion of the Indian Plate as well as induces the thrust faulting in the SFTB to the southeast (Bernard et al. 2000;Szeliga et al. 2012). According to Rowlands (1978), Humayon et al. (1991), and Bernard et al. (2000), localized strike-slip faults such as Kingri Fault interact with the contractional structures of the SFTB. Based on the structural trends within its lobate geometry, the SFTB can be divided into the eastern, central, and western parts (Jadoon and Zaib 2018). The eastern boundary of the SFTB is marked by the NS oriented Sulaiman Fig. 1 a Regional tectonic map of Pakistan showing the location of the Sulaiman Fold and Thrust Belt (SFTB) with respect to the major tectonic elements found within the Indian Plate and the Eurasian Plate (after Jones et al. 1961;Kazmi and Rana 1982;Reynolds et al. 2015). The study area is bounded by the Main Frontal Thrust (MFT) to the east and by the transpressional Zhob Valley Thrust (ZVT) and the Ghazaband Fault (GF) to the west and northwest. b A detailed map of the SFTB with the Range, which has a monoclinal expression with eastward dipping molasse strata (Humayon et al. 1991), while the western and northern boundaries of the SFTB are marked by the Zhob Valley Suture Zone (Ahmad and Abbas 1979;Allemann 1979;Gansser 1979;Farah et al. 1984;Jadoon and Khurshid 1996), and Mesozoic sedimentary rocks intermittent with numerous faults, unconformities, and lateral facies change (Jones et al. 1961;Wells 1984;Jadoon et al. 1994;Jadoon and Khurshid 1996).
The evolution of the SFTB started during the Paleogene when the Indian Plate initiated its subduction underneath the Eurasian Plate. The basin architecture and evolution of the SFTB have been further influenced by the interaction with the Afghan Block toward the west during the Cenozoic (Treloar and Izatt 1993). The northern and eastern parts of the modern-day SFTB are dominated by north-south (N-S) trending structures comprising thrust faults, strike-slip faults, and asymmetrical folds (Khan and Scarselli 2021). Stratigraphically, the SFTB is composed of a passive margin sequence of the Mesozoic platform carbonates, clastics, and volcanics. The oldest units in the SFTB, predominantly the limestones ranging from Permian-Jurassic, have been deposited on the passive margin of the Indian Plate before its separation from the Afro-Arabian margin. The Cretaceous units also represent pre-collisional deposits forming the westward sloping shelf of the Indian Plate (Sultan 1997). The dominant Cretaceous units are the deep marine shales, mudstones, and shelf carbonate deposits (termed together as Parh Group), near-shore clastic deposits (Mughalkot Formation), and deltaic deposits (Pab Sandstone) (Khan and Clyde 2013). Approximately 10-km-thick sediments have been deposited in the SFTB during the pre-rift, rift, and drift phase of the Indian Plate (Jones et al. 1961;Hemphill and Kidwai 1973). The emplacement of the Paleocene to early Eocene Muslimbagh Ophiolites over the Cretaceous shelf sediments marks the initiation of a collisional sequence (Khan and Clyde 2013). Therefore, from Eocene onwards, the newly forming Himalayas contributed toward the younger siliciclastic sediments deposition in a shallow-water deltaic environment similar to the modern Indus delta-fan system (Eames 1951;Humayon et al. 1991;Treloar and Izatt 1993;Qayyum et al. 1996Qayyum et al. , 2001Kassi et al. 2009). The SFTB can be classified into the four distinct tectonostratigraphic zones (Kamzi and Abbasi 2008; Jadoon and Zaib 2018): (1) Permo-Triassic to Eocene deep marine to shallow-marine shelf strata, (2) Cretaceous ophiolites obducted during the late Paleoceneearly Oligocene over the platform marine sequence, (3) late Eocene to early Oligocene Khojak Flysch as evidence of early Himalayan uplift and marine clastic sedimentation in the converging Tethys ocean, and (4) late Oligocene to Holocene molasse deposits along the mountain front.
The thickness distribution and mechanical properties of these sediments highly influence the deformation pattern of the SFTB (Reiter et al. 2011). The SFTB has a thin-skinned deformation style with passive-roof duplexes, that have anticlinal stacks, fault-propagation, and fault-bend folds as dominant structures (Banks and Warburton 1986;Humayon et al. 1991;Jadoon et al. 1992Jadoon et al. , 1994. The northward subduction of the Indian Plate under the Eurasian Plate and the oblique interaction of the Indian Plate with the Afghan Block in the northwest SFTB cause the repeated activation of folds and faults within SFTB. The southern culmination of SFTB is termed as Sulaiman Lobe, which is inferred to have evolved as a thin-skinned feature over a weak decoupling zone (Davis and Engelder 1985;Jadoon 1991;Jadoon et al. 1992Jadoon et al. , 1993Jadoon et al. , 1994). The regional passive-roof sequence has incipient back thrusts and pop-ups present in the foreland. Therefore, the SFTB can be defined as an uplifted and deformed (thinskinned) passive margin sequence of the northwestern edge of the Indian Plate that is accreted to the eastern edge of the Afghan Block and being under thrusted by the Indian Plate (Reynolds et al. 2015).

Study area
The study area is located within the central SFTB, which is differentiated from the eastern SFTB based on the E-W trending structures perpendicular to the tectonic transport direction. A total of nine anticlines have been reported in the study area (Jones et al. 1961) and have been labeled from A1-A9 as shown in Fig. 2. A sinistral strike-slip Kingri Fault (Fig.  1b), formed as a result of the collision of the Indian and Eurasian plates, marks the boundary between the eastern and central SFTB. Kingri Fault also takes up the component of the N-S convergence along with left lateral strike-slip motion (Rowlands 1978;Humayon et al. 1991;Bernard et al. 2000). The Kingri fault is approximately 90 km, NNW-SSE trending fault, which terminates northwards in the imbricated Cretaceous strata within the eastern SFTB and is classified as a lateral ramp (Jadoon et al. 2019). There is a marked structural trend change (from N-S to E-W) toward the western and southwestern side of the Kingri Fault because of the southward movement of Sulaiman Lobe in a thin-skinned manner that produces drag at the edge of the fold belt. The effect of this drag is very evident in the southern and southwestern parts of the study area where the structural trends transition to a E-W orientation from the N-S orientation in the eastern and northeastern parts. The study area, like the SFTB, is also believed to have undergone thin-skinned tectonic deformation, related to a regional structural detachment (Jones et al. 1961;Bender and Raza 1995). The fold geometry within the study area is interpreted as a fault-propagation fold, which is developed in response to the motion along the thrust imbricates and splays from the regional basal decollement, as well as the positive flower structures near the Kingri Fault (Rowlands 1978). The anticlines, such as A1 and A2, are cut by tear faults in the study area. The structures in the study area are divided into plunging upright folds and folds with curvilinear hinges (Fleuty 1964). The Mekhtar Anticlinorium (A7), located toward the northwestern part of the study area, is drawn with a major hinge line that portrays the major trend of this structure, and thus is termed as an anticline for this study (Fig. 2).
Detailed stratigraphic and structural studies for the study area have not yet been reported. Some regional works, undertaken by Jones et al. (1961), Bannert et al. (1992), Malkani (2010), and Jadoon and Zaib (2018), have defined stratigraphic nomenclature and structural elements within the study area. However, the present work undertakes the study at 1:50,000, which provides high resolution information than the previous works. According to Jones et al. (1961) and Malkani (2010), the oldest unit exposed in the study area is the Jurassic Loralai Formation (Jlo), which is composed of the thin to medium bedded Limestone with well-preserved joints. The Cretaceous rock units belonging to the Parh Group (Kpg), Mughalkot Formation (Kmk), and Pab Sandstone (Kpb) are also exposed in the cores of the anticlines and comprise different lithologies (shale, mudstone, marl, and sandstone, respectively). Kmk and Kpb are collectively termed as undifferentiated toward the northwestern part of the study area, as they show similar lithological characteristics (Fig. 2). Paleocene is marked by the onset of the deposition of the Dunghan Formation (Padg). The Dunghan Formation (also locally termed as the Dab Formation -Padb) shows lithological variations (having shales with subordinate limestone and sandstone beds in the eastern part, while limestone becomes dominant southwards) in the study area (Jones et al. 1961). The Ghazij Formation (Eogh), which predominantly consists of shale and sandstone, overlies the Dab Formation (Malkani 2010). The Kirther Group (Eok) is only exposed toward the southeastern parts of the study area and is predominantly composed of limestone (Jones et al. 1961). The description of the Fig. 2 Geological map of the study area incorporating different structural elements (after Jones et al. 1961). The oldest rocks belong to the Jurassic Loralai Formation (Jlo) and are exposed in the cores of anticlines in the northern and northwestern parts while the youngest rocks belonging to the Cenozoic (Eocene) Ghazij Formation (E o gh) are exposed in the cores of synclines as well as the southeastern part of the study area. The crosssection K-K′ depicts the deformational style of the study area and the location of major anticlines (after Jones et al. 1961;Bannert et al. 1992) lithological makeup of each rock unit is given in Table 1. In terms of the structural architecture, the study area exhibits a dominantly folding regime with some minor faulting induced by the Kingri Fault (Jones et al. 1961). The anticlines, located in the eastern part of the study area, comprise Cretaceous strata in their cores, while the anticlines, located in the western part, contain older Jurassic units in their cores. The synclines in the study area have Paleocene to Eocene rocks in their cores. Table 2 summarizes the anticlines that have been used to record the lineament data on the field due to their accessibility.

Methodology
The objective to analyze lineaments on the selected anticlines has been accomplished in four stages (Fig. 3). During the first stage, the joint data from the flanks of the anticlines have been acquired during a field survey of the study area in 2017-2018 (Table 2). Joints can be categorized in the outcrops as systematic and non-systematic sets of definite and random orientations (Suppe 1985), constituted by the mechanical discontinuities and/or mechanical behavior of the rocks (Nelson 1979). For this study, only systematic joints are considered. Joint and bed attitude data has been collected during the field campaign using a Brunton Compass, while the joint spacing data has been acquired using the scanline method (1m line) on different traverses (Priest and Hudson 1981).
The second stage relates to the extraction of lineaments using the Digital elevation models (DEMs) generated from the Shuttle Radar Topography Mission (SRTM). SRTM-DEM is preferred over the other satellite data such as the ASTER-GDEM because of its superior vertical accuracy (Elkhrachy 2018), which helps in the extraction of lineaments in the rugged terrain and has been successfully employed in the lineament studies (Lee and Moon 2002;Solomon and Ghebreab 2006;Masoud and Koike 2011;Das and Pardeshi 2018). SRTM-DEM with 30 m spatial resolution has been acquired from the United States Geological Survey (USGS) Earth Explorer Portal (http://earthexplorer.usgs.gov) and has been subjected to the shaded relief mapping using the directional filtering (Suzen and Toprak 1998), in the ArcGIS 10.7.1, which sharpens the boundary or discontinuity between the adjacent areas and helps in identification of the lineament features, such as the straight valleys, straight streams segments, and rock boundaries. The shaded relief maps have been generated with four azimuth angles (0°, 45°, 90°, 135°) to ensure unbiased extraction of all the lineaments in the study area. LINE module in PCI-Geomatica has been used for the automated extraction of the lineaments from the shaded relief images. LINE module, which is based on edge detection, thresholding, and curve extraction, extracts the lineaments and converts them into vector format based on the input Table 1 Stratigraphic makeup of the SFTB and the study area (after Jones et al. 1961;Kassi et al. 2009 Table 3). These input parameters determine the number and length of the extracted lineaments (Thannoun 2013). Furthermore, the lineaments have been spatially constrained to anticlines by eliminating synclinal lineaments as well as the non-geological lineaments such as roads, canals, and other cultural features using the topographic maps and high-resolution satellite imagery. The finalized lineaments from four azimuth angles have been merged in the ArcGIS to generate a single vector shapefile, which has been further subjected to the duplicate lineaments removal to obtain the final lineament dataset of the study area. The interpretation of these lineaments was done at 1:50,000 scale to statistically correlate these lineaments with the field joints dataset, i.e., JS1=LS1, JS2=LS2, and JS3=LS3. The mean orientation value of the field joint dataset for an individual anticline, along with its standard deviation as a threshold, has been chosen to highlight and extract the lineament set that corresponds to that particular joint set.
In the third stage, the field joint dataset has been correlated with the lineament dataset to quantify the relationship between the observed and extracted values. Spearman's rank correlation coefficient has been adopted to observe if the lineament and joint datasets are statistically similar (Dash and Dash 2017), after the Shapiro-Wilk test of the normality (p < 0.05) disallowed the use of parametric tests. The lineament dataset has been used to generate a lineament density map using the ArcGIS. The lineament density is defined  Fig. 3 Flowchart depicting the methodology adopted in the study. Two major datasets (field observations and remotely sensed observations through the DEM) are correlated in terms of the error and further subjected to the multivariate statistical analysis, i.e., cluster analysis (CA) and principal component analysis (PCA), along with the lineament density analysis to identify the relationships between different lineament sets as the combined length of all lineaments over a unit area (1 km 2 in the present study) (Greenbaum 1985). The histograms and rose diagrams have been used to graphically portray the lineament density and lineament frequency using Grapher 12.0. The interpretation of these lineaments from graphical representations can help recognize and map the surface structural setup of a specific area (Abdullah et al. 2010).
The final part of the study implements two multivariate statistical techniques on the joint and lineament datasets to analyze the lineament architecture in the study area. Cluster analysis (CA) associates and connects the similar observations by normalizing the dataset to compute Euclidean distance using the different algorithms to yield a dendrogram that clusters similar observations together (Davis and Sampson 1986;Spanos et al. 2015;Adhikari and Mal 2019). Principal component analysis (PCA) defines and explains the variance in a large dataset by indicating the associations within a dataset and thus reducing dimensionality by obtaining the eigenvalues and eigenvectors (loadings) from the original observations, which relates observations to certain factors (Ravikumar and Somashekar 2017). All the statistical analysis has been performed using the opensource PAST software package (Hammer et al. 2001).

Field Joints
The 57 sampling points, taken on the nine anticlines, have been measured for the joint orientation, joint dip, joint spacing, bedding strike, bedding dip, and lithology (Table 4). The non-homogenous number of observations (n) on different anticlines depicts the accessibility to different parts of those anticlines. A3 and A7 have the highest number of data points while A1 and A5 have the lowest data density. The field data shows that certain joints on the anticlines have a similar geometrical makeup and have been thus classified into joint sets. In general, 3 joint sets have been observed within the study area anticlines (i.e., JS1, JS2, and JS3); however, only 2 joint sets are continuous on all anticlines (i.e., JS1 and JS2), where each set is characterized based on its relationship with the orientation of the anticlines. The JS1 and JS2 are orthogonal to each other while JS3 is oblique to both JS1 and JS2.
The orientations of the JS1 for all the anticlines vary between 011 and 074 azimuth (NE-SW) with a standard deviation of ± 021, an average dip of 64°, and an average joint spacing of 56 cm. The orientations of the JS2 vary between 279 and 346 (NW-SE) with a standard deviation of ± 017, an average dip of 65°, and an average joint spacing of 96 cm. The Table 4 Quantified average results for field observations and the remotely sensed observations (n). Field observations (n) include the major % of the lithology in which most measurements have been taken, bed attitude, orientation, dip, and spacing (cm) for JS1, JS2, and JS3.
Remotely sensed observations have been classified according to joint data as LS1, LS2, and LS3 and are only analyzed in terms of their orientation. (Sh= shale, Ss= sandstone, Mrl = marl, Lst = limestone) orientations of the JS3 have been measured only for anticlines A1, A2, A3, A4, and A7 and vary between NNE-SSW to WNW-ESE with an average dip of 64°. Based on the relationship between the bedding strike and dip, we can conclude that the JS1 is parallel to the bedding strike, JS2 is perpendicular to the bedding strike, while the JS3 is oblique to the bedding strike. Figure 4 plots the JS1, JS2, and JS3 on rose diagrams and histograms to visually depict the direction of the joint distribution within the study area. The histograms show that the NE quadrant has the highest frequency for the JS1 in the 050-060 bin, while the NW quadrant has the highest frequency for the JS2 in the 320-330 bin. The analysis of the rose diagrams shows that the predominant orientation(s) varies with the location of the data points on the anticlines. Furthermore, due to the curvilinear axis of some anticlines, the joint sets also change their orientation within the respective quadrant. The dips are higher toward the axial part and gentle toward the flanks of the anticlines. However, higher dips have been observed on the flanks of anticlines (such as A1 and A3) which are close to the local faults.

Lineament datasets
The lineaments have been classified as the LS1, LS2, and LS3, where each corresponds to the respective field joint dataset, by defining the thresholds for the lineaments through the mean value of each JS and its one standard deviation value. A total of 7300 lineaments have been obtained after the subsequent processing and classification steps (Table 4). Table 5 represents different thresholds and their effect on predicting the high correlation of the lineaments with the joints. Selecting a threshold of ± 005 only gives an accuracy between 25 and 31% while a threshold of ± 010 increases the accuracy to 72-78%. Therefore, 10°bins have been used throughout this study. The total length of lineaments estimated from the final output is 2863 km, with 3263 lineaments belonging to the LS1, 3037 belonging to LS2, and 218 belonging to the LS3. Figure 5 shows the lineament map for the nine anticlines in the study area and their distribution based on the LS1, LS2, and LS3, while their respective data has been quantified in Table 4. The orientations of the LS1 vary between 031 and 075 azimuth (NE-SW) with a standard deviation of ± 010 while the orientations of the LS2 vary between 301 and 345 azimuths (NE-SW) with a standard deviation of ± 012. The orientations of the LS3 for A1-A4 vary between 009 and 050 (NNE-SSW) and for A7 between 275 and 290 (WNW-ESE). The histograms in Fig. 5 show that the NE quadrant has the highest frequency for the Fig. 4 A map of the data collection points throughout the study area. Bed attitude (strike and dip amount) data along with the joint set data (classified into JS1 as red, JS2 as green, and JS3 as blue) has been collected in and around nine anticlines in the study area. The orientations of the different joint sets have been measured, from 000 to 360, and plotted in the histograms to analyze the distributions of different joint sets with respect to their orientation. The histograms depict the joint sets in two quadrants (NE and NW). These sets are also displayed in the rose diagrams for the respective anticlines. The base map in the background is the shaded relief map with the NE filter which illuminates the morphological cum structural makeup of the area LS1 in the 050-060 bin, while the NW quadrant has the highest frequency for the LS2 in the 290-310 bin. The analysis of the rose diagrams constructed for the lineaments shows that the predominant orientation(s) varies with the orientation of the fold axis of the anticlines. For anticlines, A1, A2, A5, and A6, the dominant orientations for the LS1 and LS2 are NE-SW and NW-SE; for anticline A4, the orientations for the LS1 and the LS2 are ENE-WSW; for anticlines, A3, A8, A9, the orientations vary between NE-SW and ENE-WSW for the LS1, while NW-SE and WNW-ESE for the LS2.

Lineament density mapping
The joint and lineament datasets, plotted on a scatter plot, show a very high correlation coefficient (0.97) with a low p value (< 0.05) signifying a statistically significant correlation (Fig. 6). This validates and justifies the use of 7300 widely distributed lineaments for density mapping along with the spatially limited and non-homogeneously distributed joint dataset. The lineaments, plotted on all nine anticlines, have been subjected to the lineament density mapping using the ArcGIS (Fig. 7). In the study area, the lineament density is Fig. 5 A map of the lineament data extracted from the shaded relief map of the DEM. The lineament data (classified into LS1 as red, LS2 as green, and LS3 as blue) is represented in the histogram for the entire study area and also represented in the rose diagrams for the respective anticlines. The base map in the background is the shaded relief map with the NE filter which illuminates the morphological cum structural makeup of the area Table 5 Average threshold calculation of the lineament sets (LS) with respect to the standard deviation of the joint sets (JS). A lower standard deviation (005) has a lower correlation and thus lower prediction for the LS dataset. A higher standard deviation (010) has a higher correlation and thus higher prediction for the LS dataset. A standard deviation greater than 010 will yield a higher error and thus will reduce the prediction found to be varying from 0 to 8.31 km/km 2 and has been classified into four classes, i.e., very low (< 1 km/km 2 ), low (1 to 2 km/km 2 ), moderate (2 to 5 km/km 2 ), and high (> 5 km/ km 2 ). The lineaments with a moderate density are Fig. 6 A statistically significant correlation on 57 field locations comprising of 122 observations (JS1, JS2, and JS3) and their equivalent remotely sensed observations, depicting the nonnormal distribution of data (low Shapiro-Wilk p value) with a high correlation coefficient of 0.97. The green shaded area depicts the 95% confidence interval of the data set Fig. 7 Lineament density map showing four density classes for the respective anticlines. Respective lineaments, found within each respective zone of the anticlines, have also been plotted on the rose diagrams to analyze the trends for the lineament of each density class predominant with 52% of the area coverage, while the low density lineaments cover about 19% of the study area ( Table 6). The very low and high density classes cover the least area, i.e.,14% and 15% respectively.

Multivariate statistical analysis
Due to the unavailability of the JS3 data on various anticlines, it was removed altogether for the multivariate data analysis. The cluster analysis (CA) shows that the four distinct groups of anticlines can be considered similar based on the field joint datasets with variables such as the orientation, spacing, and dips of the JS1 and JS2 datasets. The clusters are classified as: (1) A2, A6; (2) A5, A8; (3) A1, A3; (4) A4 A7, A9 (Fig. 8a).
The principal component analysis (PCA) for the standardized joint and lineament sets utilized the clusters from CA. The eigenvalues and the subsequent variance illustrate that only the two principal components (PC) axes are responsible for the high variance in the datasets ( Table 7). The PCA for the joint datasets shows the clusters of anticlines and the biplots of the variables that signify the spatial trends (Fig. 8b). Cluster-1 for the joint dataset has high dips for the JS1 and JS2, NNE-SSW dominating orientation due to lower azimuth values in both quadrants and relatively moderate joint spacing. Cluster-2 is categorized by moderate dips, NNE-SSW dominant orientation, and high joint spacing. Cluster-3 is characterized by high joint dips, N-S dominant orientation, and low to moderate joint spacing. Cluster-4 is characterized by the ENE-WSW dominant orientation with low joint dips. Based on the eigenvalues, joint dips have been found as the dominant factor controlling the highest variance in the joint dataset. Based on the geological setup of the area, the PC1 is inferred as proximity to a fault for a particular anticline while the PC2 is inferred as the mechanical compaction/strength of rock units comprising the anticlines.
The PCA for the lineament dataset shows three clusters incorporating nine anticlines and the biplot of variables that signifies the spatial trends (Fig. 8d). The variables for the lineament dataset include the orientation, length (km), and the total area of an anticline (km 2 ) for which the lineaments have been computed. Cluster-1 is found to have relatively a low lineament length and a lower area of the anticlines with a dominant NNE-SSW orientation. Cluster-2 has a dominantly ENE-WSW orientation with a low lineament length and a low anticline area. Cluster-3 is characterized by a high lineament length and anticline area with a dominantly NE-SW orientation. Based on the eigenvalues, the lineament length (km) has been found as the dominant factor controlling the highest variance in the lineament dataset (Table 7). Based on the geological setup of the area, the PC1 is inferred as an interlimb angle for an anticline while the PC2 is inferred as the hinge line/fold axis curvature for the anticline.

Joint characterization and implications on the remotely sensed lineaments
In order to analyze the lineaments and their implications for the study area, it is extremely important to characterize joints based on the spatial and geometrical relationships. Joints, being more widespread and easily measurable on the field, are more likely to manifest themselves as lineaments on the satellite data as compared to the other lineaments. Many researchers such as Cloos (1947), Price (1966, and Stearns (1968) have proposed methods to establish the relationships between the orientation of joint sets and a homogeneous state of stress for the folded layers or fold axis. Figure 9 represents an outcrop with comparable joint sets, as described by Watkins et al. (2018).
Joints that form during a folding episode exhibit orientations related to fold-related stresses while those formed before or after a folding episode exhibit regional stress patterns. To differentiate between the two types requires a rigorous interpretation of the relationship between the timings of folding and jointing (Hancock 1985;Dunne 1986). According to Price and Cosgrove (1990), if the symmetrical relationships cannot be established between the fold geometry and the joint set orientation, the joint sets are not inferred to form during folding. An assumption to this model was that no joints existed before the folding event and had no influence on joint development during the folding phase in a fold-thrust belt. A parallel strike of the joints to the folding geometry testifies to their genetic association, however, does not conclude directly about their timing of formation. Therefore, the joints have been characterized in our study area based on their orientation and genetic associations and are classified into joint sets based on the mean orientation, orientation distribution, spacing, and relative chronology.
In the thrust-related anticlines, similar to our study area, there can be up to four joint sets (Price 1966). There are three Fig. 8 a Cluster Analysis plot for the standardized field data showing anticlines clustered together, i.e., A2, A6; A5-A8; A4, A9, A7; and A1, A3 are different sets identified from the CA and are given respective colors. b Principal component analysis (PCA) plot of the (a) with its biplot showing variables such as the joints dip, orientation, and spacing with two major principal components (PCs) depicting the proximity to a fault and mechanical compaction of rocks. Anticlines conform to the results of CA and show distinct clusters depicting the PC properties. c Cluster analysis plot for the standardized remotely sensed data showing anticlines clustered together, i.e., A2, A6, A1, A5; A4, A7; and A9, A8, A3 are different sets identified from the CA and are given respective colors. d PCA plot of the (c) with its biplot showing variables such as the lineament orientation, length, and surface area of the anticline measured from its surface hinge, with two major principal components (PCs) depicting the fold limb angle and curvature of the fold hinge. Anticlines conform to the results of CA and show distinct clusters depicting the PC properties joint sets (JS1, JS2, and JS3), that have been observed in the field, which form a geometrical relationship with the folding geometry of the anticlines. The overall strike of the JS1 has been observed to be parallel to the fold hinges of the anticlines while the overall dip is observed to be normal to the bedding (Fig. 10a). JS1 may form due to the development of localized tensional stresses in the same orientation as that of regional compression. The overall strike of the JS2 is found to be perpendicular to the fold hinges, while the overall dip is normal to the bedding (Fig. 10c) and may be associated with a localized extension due to elevated curvature on a plunging anticline. JS3 is interpreted as a conjugate joint set, that may form due to a regional compression or localized inner-arc compression associated with a tangential longitudinal strain folding (Ramsay 1980).
The cross-cutting relationship of these three joint sets may allow for the determination of their relative age. The joints formed earliest tend to be long and continuous (JS1, JS2). One field observation (Fig. 9) shows that the JS1 abuts JS2, thus suggesting that the JS2 is an older generation set than the JS1. However, such an abutting relationship between the JS1 and JS2 is not well recorded for other anticlines. With the development of the shorter and younger joints (JS3), the preexisting joints modify the stress orientation, and this results in a poor correlation of the younger joint set with the older, in which the former abuts against the latter. The systematic relationships between the orientations of the joint sets (JS1, JS2, JS3) with the folding geometry and their abutting relationships can define the folding related stresses throughout an anticline (Fischer and Wilkerson 2000), assuming these sets are synfolding. For example, the JS3 is observed to be neither parallel nor perpendicular to the fold axes and thus may be induced by fold-related strains. Furthermore, the orientations of the joints on certain anticlines, such as A3 and A9, have also been influenced by the rotation of the fold hinges. The three joint sets are therefore interpreted as a reflection of the fold growth history. The lineaments, due to their high correlation with the joints, are assumed to have similar geometrical relationships with the folds in the study area, i.e., JS1=LS1, JS2=LS2, and JS3=LS3. It is noteworthy that the different joints form on different positions of an anticline and due to the limitations in data collection and area accessibility, only some selective parts of the anticlines have been used for collecting joint data and therefore this study only focuses on the relevance of joints related to folding episode of SFTB. Therefore, the limited joint data can be confidently spread over the whole anticlines, thus revealing the overall geometry of lineaments with respect to the folding geometry.

Controls on lineament density
The variations in the lineament density can be attributed to the lithological and structural controls (Nelson 1985). Very low to a low density category is found continuously on the outer flanks and hinges of the anticlines, while the cores of the anticlines have high lineament density. Lithological controls, such as mechanical compaction, rock brittleness, mineralogical composition, texture, porosity, and bed thickness, induce a higher degree of jointing (Hugman and Friedman 1979;Marshak and Mitra 1988). Lineament density is also found to be structurally controlled, i.e., a higher curvature area of a fold (forelimb and crust) has a higher lineament density (Barbier et al. 2012;Awdal et al. 2013). In the study area, within high fold curvature regions, the LS1 is dominant and can be correlated with the JS1, which may have formed during the folding episode of SFTB (Table 4). Low curvature regions (low strain domains) of the anticlines, such as the gentle flanks, have longer lineament lengths with variable lineament density (Florez-Niño et al. 2005). These domains, in the study area, are dominated by two orthogonal sets: parallel (LS1) and perpendicular (LS2) to the fold hinge, which form a weak abutting relationship and are thus interpreted to be contemporaneous. Furthermore, the LS3, which is oblique to the dominant lineament sets abuts against LS1 and LS2 through shearing and are thus interpreted to be younger. Within the study area, A1 and A2 show this younger lineament set due to the presence of the nearby Kingri Fault.
According to Watkins et al. (2015), the structural control on the lineaments dominates in the regions with higher strain, while lithological control on the lineaments dominates in the regions with lower strain. This has been confirmed in the lineament density map (Fig. 7), where A3 and A4 have higher lineament density in the core as compared to the other anticlines. The presence of a high lineament density zone within these anticlines is probably due to their proximity from a local fault which makes this part of the fold a high strain domain. Furthermore, in these anticlines, the presence of thin mudstone beds of Parh Group in the cores may also contribute toward the enhancement of the lineament density (Marshak and Mitra 1988). Lithological control is corroborated by the PCA plot where limestone dominant anticlines (A5, A8) appear together based on their high joint spacing and low joint density (Fig. 8b). The structural control is depicted on the PC1 while the lithological control is represented by the PC2. Figure 8d also shows A7 and A4 having a low interlimb angle, which indicates the presence of high strain domains that enhances the lineament density. In general, a high lineament density reflects a high degree of jointing and a higher intensity of deformation (Edet et al. 1998;Hung et al. 2002). Furthermore, some isolated regions of high lineament density areas on the flanks of some anticlines, i.e., A1, A2, and A7, reveal that some local faults could also contribute toward the enhancement of the lineament density.  Fig. 9 A comparison of field joint set with the graphical description (after Watkins et al. 2018). The three joint sets correspond to the folding geometry

Regional implications
Biostratigraphic constraints from northern Pakistan (Beck et al. 1995), the timing of emplacement of ophiolites within the Sulaiman Range (Allemann 1979;Gnos et al. 1997), and Indian Plate motion reconstruction (Copley et al. 2010) puts an age of ∼ 55 Ma for the initial collision of India with Asia in the region of the Sulaiman Range. Since this time, multiple different tectonic phases have intensively affected the rocks in the SFTB. The oblique convergence of the Indian Plate with the Eurasian Plate/Afghan Block has induced a convergence component that has resulted from the change in the angle between the orientation of the plate boundary and the direction of plate motion. This convergence component is also represented as different strike-slip features within the SFTB, which are oblique to the direction of the plate motion. Furthermore, local stress regimes are also responsible for determining the type and degree of deformation. The central and northern parts of the SFTB are dissected by the NW-SE trending faults which result from the uplifting of the Sulaiman Range, dated around the early Miocene with a major period of deformation beginning in the Pliocene (Jones et al. 1961;Humayon et al. 1991;Qayyum et al. 1997;Kassi et al. 2009;Kasi et al. 2012). These events have also affected the lineament makeup in the study area. The three joint and lineament sets, that have been delineated from this study, are represented in a rose diagram to visualize their orientation and statistical distribution (Fig. 11 a and b). The regional stress field, due to the collision at the NW Indian Plate, produces NE-SW oriented lineaments (LS1) in the study area which have the lowest normal stress across them. The major stress direction is toward the SE and acts perpendicular to fold hinges present in the study area. The arrangement of the stress in the NW-SE (310-320) (marked with an arrow showing the maximum principal stress axes σ1) parallels with the orientation of the LS2. Such an arrangement of the LS1 and LS2 dominates the study area and points toward the fact that the lineaments are produced during the folding phase, which was initiated due to the main stress arising from the Indo-Eurasian collision. However, for the anticlines, such as A3 and A9, the LS1 is not homogenous, i.e., the orientation of these lineaments within these folds change from WNW to NNW, implying that the local stresses generating from the local faults also play a role in the deformation of the area. This effect is also evident for the LS3, which orients itself parallel to a local shear, i.e., Kingri Fault. The LS3 is interpreted to be a result of a shear, that could be a consequence of a younger deformational event which might have altered the fold trends and could also have created the LS3 by introducing transpressional structural geometries such as high angle thrusts, blind thrusts, and positive flower structures.
In terms of the validation for this study, the orientation of the maximum stress, based on both the trends of folds and the lineaments, is consistent with the regional thrust fault focalmechanism solutions of earthquakes and the resultant stress tensors (Bernard et al. 2000;Reynolds et al. 2015;Heidbach et al. 2016;Jadoon and Zaib 2018). Fig. 11 a A cumulative rose diagram incorporating the joint data showing JS1 perpendicular to the maximum stress direction while the JS2 parallel to the maximum stress direction. JS3 has both parallel and perpendicular trends. The orientations are scattered due to the nonhomogenous distribution of field data. b A cumulative rose diagram

Conclusions
The orientation of the joint dataset from the field has shown a successful correlation with the remotely sensed lineaments, thus allowing for a regional distribution of the lineaments in the study area. Three lineament sets, that are comparable to three joint sets (LS1=JS1, LS2=JS2, LS3=JS3), have been found to have a symmetrical relationship with the folding geometry. With the fold hinge, the LS1 shows parallel, the LS2 shows perpendicular, while the LS3 shows an oblique geometry, implying a relationship of the lineaments with the regional folding episode as well as a local deformation. Structural and lithological controls on the lineament sets, established through multivariate statistics, indicate that the less competent lithologies, as well as proximity to a local fault, impart higher lineament density. Moreover, the curvature of the fold hinge, depicting a local deformation event, also controls the orientation of the lineaments.
This study establishes lineaments and their relevance to the regional stress regimes of the area by relating the lineaments to the folding episodes and using the folding episodes to delineate the regional stress direction. Orientations of the lineament sets, when correlated to the regional principal stress axis, reveal that the folding event is directly related to the oblique convergence of the Indian Plate with the Afghan block. The study yields NW-SE regional stress direction for the study area which agrees with the folding geometries and the published literature. This is a novel study for this part of SFTB with sparse field joint data due to inaccessibility to the various parts of different anticlines. The pre-folding jointing phenomenon has not been considered toward the interpretation of the joints due to the lack of sufficient and homogenous field observations on all the anticlines. The authors recommend the collection of homogenous field data on all parts of the anticlines for the establishment of additional lineament sets, which will help in improving the overall deformational model of the SFTB.
Acknowledgements The authors would like to acknowledge the Geological Field Party-1 of the Oil and Gas Development Company Limited (OGDCL), Pakistan, without which the field data acquisition in such a remote area would never have been possible. The authors would also like to acknowledge Mr. Shakeel-ur-Rehman from the OGDCL for his thorough review of the structural elements of the study area and Dr. David Peacock (Georg-August-University, Göttingen, Germany) for his valuable advice on the technical nomenclature. The authors would also like to thank Dr. Olaf Lenz (TU Darmstadt, Germany) who provided his valuable comments toward the statistical analysis undertaken in this study. Finally, the authors are grateful to the reviewers who gave detailed feedback which significantly improved this manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.