*Article* **Integrating a GIS-Based Multi-Influence Factors Model with Hydro-Geophysical Exploration for Groundwater Potential and Hydrogeological Assessment: A Case Study in the Karak Watershed, Northern Pakistan**

**Umair Khan 1,2, Haris Faheem 1, Zhengwen Jiang 1, Muhammad Wajid 1, Muhammad Younas <sup>2</sup> and Baoyi Zhang 1,\***


**Citation:** Khan, U.; Faheem, H.; Jiang, Z.; Wajid, M.; Younas, M.; Zhang, B. Integrating a GIS-Based Multi-Influence Factors Model with Hydro-Geophysical Exploration for Groundwater Potential and Hydrogeological Assessment: A Case Study in the Karak Watershed, Northern Pakistan. *Water* **2021**, *13*, 1255. https://doi.org/10.3390/ w13091255

Academic Editors: Andrzej Wał ˛ega and Tamara Tokarczyk

Received: 5 April 2021 Accepted: 29 April 2021 Published: 30 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** The optimization of groundwater conditioning factors (GCFs), the evaluation of groundwater potential (GWpot), the hydrogeological characterization of aquifer geoelectrical properties and borehole lithological information are of great significance in the complex decision-making processes of groundwater resource management (GRM). In this study, the regional GWpot of the Karak watershed in Northern Pakistan was first evaluated by means of the multi-influence factors (MIFs) model of optimized GCFs through geoprocessing tools in geographical information system (GIS). The distribution of petrophysical properties indicated by the measured resistivity fluctuations was then generated to locally verify the GWpot, and to analyze the hydrogeological and geoelectrical characteristics of aquifers. According to the weighted overlay analysis of MIFs, GWpot map was zoned into low, medium, high and very high areas, covering 9.7% (72.3 km2), 52.4% (1307.7 km2), 31.3% (913.4 km2), and 6.6% (44.8 km2) of the study area. The GWpot accuracy sequentially depends on the classification criteria, the mean rating score, and the weights assigned to GCFs. The most influential factors are geology, lineament density, and land use/land cover followed by drainage density, slope, soil type, rainfall, elevation, and groundwater level fluctuations. The receiver operating characteristic (ROC) curve, the confusion matrix, and Kappa (K) analysis show satisfactory and consistent results and expected performances (the area under the curve value 68%, confusion matrix 68%, Kappa (K) analysis 65%). The electrical resistivity tomography (ERT) and vertical electrical sounding (VES) data interpretations reveals five regional hydrological layers (i.e., coarse gravel and sand, silty sand mixed lithology, clayey sand/fine sand, fine sand/gravel, and clayey basement). The preliminary interpretation of ERT results highlights the complexity of the hydrogeological strata and reveals that GWpot is structurally and proximately constrained in the clayey sand and silicate aquifers (sandstone), which is of significance for the determination of drilling sites, expansion of drinking water supply and irrigation in the future. Moreover, quantifying the spatial distribution of aquifer hydrogeological characteristics (such as reflection coefficient, isopach, and resistivity mapping) based on Olayinka's basic standards, indirectly and locally verify the performance of the MIF model and ultimately determine new locations for groundwater exploitation. The combined methods of regional GWpot mapping and hydrogeological characterization, through the geospatial MIFs model and aquifer geoelectrical interpretation, respectively, facilitate decision-makers for sustainable GRM not only in the Karak watershed but also in other similar areas worldwide.

**Keywords:** multi-influencing factors (MIF); vertical electrical sounding (VES); electrical resistivity tomography (ERT); groundwater resource management (GRM); hydro-stratigraphy; well logs

#### **1. Introduction**

Increasing anthropogenic repression, climate change, and environmental problems are affecting the supply and demand of domestic and irrigation water. The efficient and innovative use of geospatial and geophysical datasets for understanding groundwater management and hydrological processes in various climatic and vegetation regimes under topographical, geological, hydrological, and land-covered influence has become an important challenge, which offers a wide range of research opportunities [1–5]. There are several conventional geological, geophysical, and hydrogeological methods, and the most commonly used methods are geophysical, but they are time-consuming and mainly applicable on a small scale [6,7]. However, remote sensing (RS) and geographical information system (GIS) provide spatial, temporal, and spectral data availability that can cover large and inaccessible areas within a short period and serve as a useful tool for assessing and managing groundwater resources [8–12].

The groundwater potential (GWpot) is influenced by multiple geological, hydrological, and land-covered processes [10,12,13]. Usually, the occurrence and movement of surface water and groundwater could be assessed by optimized groundwater conditioning factors (GCFs), i.e., rainfall, lineament density, slope, soil types, drainage density, land use/land cover, lithology, elevation, and groundwater level fluctuation. [14,15]. GIS and RS analysis are useful for large-scale estimates of surface water and groundwater. Several methods have been employed to monitor GWpot, such as cumulative rainfall departure (CRD), Monte Carlo (MC) simulation, frequency ratio (FR), certainty factor (CF), weightsof-evidence (WoE), fuzzy logic index models, logistic regression (LR) model, analytical hierarchy process (AHP), and multi-influence factors (MIFs) [8,16–23]. The CRD is a water balance method which defines groundwater level fluctuations in shallow aquifers as a function of rainfall. The statistical methods (e.g., FR, LR, WoE) estimates the coefficient for each GCF by defining the relationship between the dependent variable and independent variables, while the AHP assigns a score to each conditioning factor based on expert's opinion. The MC simulation is considered to be the main tool to quantify the uncertainty in groundwater predictions. To reduce the mathematical complexity by incorporating a decision-making reasoning process based on expert system judgment, the MIF technique has become a useful GWpot modeling approach, that can quickly, accurately, cost-effectively, and consequently monitor GWpot [23–25]. MIFs constitute a GIS-based multi-criteria decision-making (MCDM) technique that enumerates the spatial relationships between dependent and independent variables according to scores assigned based on major and minor GCFs influencing GWpot [24,26]. This method is economical as it relatively simple and useful for practical applications before starting an expensive field survey [3,9,20]. It helps in narrowing down the potential areas for conducting detailed hydrogeological and geophysical surveys and ultimately locating the drilling sites [7,27].

Hydro-stratigraphy and hydrogeology are essential for characterizing aquifer potentiality and developing hydrological models to predict groundwater resources for future availability [28,29]. For geoscientists, finding and locating the source and availability of the groundwater in a complex area with multiple hydrogeological features is a vital task. Although surface geophysical measurements can provide effective spatial coverage services [30,31], these measurements depend on the area extent to be investigated, cost, geological condition, and the acquired data readability. They contribute information on groundwater levels, hydrogeological behaviors, and corresponding lithology, ensuring a higher positioning accuracy for groundwater resources [32–34]. With the proper GWpot and hydrogeological evaluation, geophysical techniques can be combined to improve efficiency. Specifically, the electrical resistivity techniques are well established and commonly used to solve numerous geological and environmental problems [35,36], which are considered as the most effective geophysical methods for the characterization of GWpot and hydrologic stratigraphy. These methods are widely used to scrutinize high-resistance and low-resistance layers, and are, therefore, valuable tools for studying aquifer vulnerability [32,37]. The quantification of the aquifer potential analyzed by VES-based reflection

coefficient, isopach, and resistivity mapping can directly verify the predicted result of the MIF model and its performance. The combination of vertical electrical sounding (VES) and electrical resistivity tomography (ERT) methods produces a high ratio of 90% compared to 82% for the VES method and 85% for the ERT method [33,38]. The spatial distribution of aquifer hydrologic characteristics, such as resistivity, reflection coefficient, overburden thickness, hydraulic conductivity, and specific productivity, plays an essential role in assessing and managing GWpot [39,40]. Apparent resistivity and reflection coefficients are the most critical hydrogeological data needed to manage groundwater resources [41]. These parameters also outline variances in the hydrological strata that help to explain aquifer models for GWpot modeling. In addition, geophysical well logging also generates useful information about the geological structure and the formations' lithology [22]. The feasibility study of resistivity surveys through boreholes has been used worldwide and is supported by general hydrogeological studies. Drilling (machinery type deployed subsurface soil/rock conditions) and electrical logs record the true location of the aquifer and corresponding lithology.

The phenomenon of surface water resource depletion and irregular spatial-temporal distribution of precipitation have made groundwater a vital natural resource for the reliable and economic provision of potable water supply in low- and mid-income regions of the Karak watershed, Northern Pakistan. In this context, this study addresses the applicability of comprehensive MCDM-MIFs model with optimized GCFs for GWpot assessment and hydro-geophysical investigation for hydrogeological characterization. As the GWpot mapping depends on the suitable GCFs and the weights assigned to them, various GCFs, such as geology, lineament density, land use/land cover, drainage density, slope, soil type, rainfall, elevation, and groundwater level fluctuations, were processed and optimized through geospatial analysis in GIS environment. The predicted GWpot results using the MIF model were then analyzed by the receiver operating characteristic (ROC) curve and the confusion matrix, and Kappa (K) analysis. However, groundwater is an invisible resource that is difficult to measure or quantify directly. Therefore, the interpretation of VES and ERT data was employed to predict hydrogeological properties, aquifer hydraulic characteristics and GWpot zones for future exploitation and installing tube wells for its utilization. Moreover, our methodology not only improves the reliability of the integrated geospatial and geoelectrical modeling and bridges the gap of GWpot evaluation and hydrogeological characterization in the Karak watershed, but also provides an optional solution of groundwater assessment in other similar areas worldwide.

#### **2. Study Area**

#### *2.1. Physical Geographical Background*

The study area is located at geodetic coordinates between the latitudes of 32◦46 and 33◦22 N and between the longitudes of 70◦43 and 71◦33 E, covering an area of approximately 2372 km2 (Figure 1b). A 123 km road from Peshawar on the Indus Highway leads to Karachi and is easily accessible from various parts of the country via mettled roads (Figure 1a). Geographically, the Karak watershed is located in the southern part of the Kohat Plateau of the upper Indus basin Pakistan. The Kohat Plateau itself lies between 70◦–74◦ E and 32◦–34◦ N, covering an area of approximately 10,000 km2. Most of the region's climate is semi-arid, with two major seasons, i.e., the rainy season and the dry season. Precipitation is the primary source of groundwater replenishment where the average precipitation is 450 mm/year, and the minimum and maximum average temperatures in the Karak (at an altitude of 706 m) are 10.3 ◦C and 43.5 ◦C, respectively, varying by altitude. The harvest depends on the amount of precipitation or pipeline well supply. Annual precipitation in the northeast ranges from 500 to 750 mm. In the study area, rainfall from June to November is 68% and is 32% from December to May. During the short rainy season, rainfall is scarce, unstable, and concentrated, and it is relatively or absolutely dry for the rest of the time. High temperature and rainfall intensity cause large amounts of precipitation loss due to evaporation and runoff, respectively [42]. The highest elevation

area of the Karak watershed is in the eastern Surghar Shinghar ranges (Figure 1b), where elevations typically exceed 1415 m above sea level. The lowest elevation area is associated with the Bannu boundary, where the river level is below 305 m.

**Figure 1.** (**a**) Generalized physical geographical features of the Potwar region; and (**b**) Location map of the Karak watershed with the surveyed boreholes and vertical electrical sounding (VES) and electrical resistivity tomography (ERT) measurements.

#### *2.2. Geological Background*

A regional geological map of the study area was prepared to plot major geological structures and lithological units (Figure 2). The Karak watershed is part of a large intermontane basin where sedimentation has taken place from weathering and erosion of the surrounding Bannu mountain belts [43,44]. The Bannu basin is located in a depression behind the Trans-Indus current uplift boundary, which leads to the formation of the Bhittani, Khisor/Marwat, and Shinghar mountains. The basin is formed by the uplift boundary from the Kohat mountain range to the Bhittanni and Marwat/Khisor mountain ranges [43], as shown in Figure 2. In the Potwar Plateau and the adjacent Kohat Plateau, the exposed sedimentary formations are Eocene limestone, evaporite, and red beds [45]. Subsurface deposits of the area widely vary from very coarse sediments (such as gravel and boulders) to very fine sediments (such as silt and clay). There are three types of sediments in this region, including alluvial fans, floodplains, and basin-filled sediments [46,47]. An alluvial fan is composed of various proportions of boulders, gravel, sand, silt, and clay. The sediments in the floodplain are mainly clay and silt, with minor amount of sand. Sandy sediments were primarily formed in the Marwat range, mainly due to erosion [48]. The ages of the exposed strata in the study area range from the Precambrian to the Quaternary. The lithological distributions of the Karak watershed are illustrated in Table 1.

**Table 1.** Lithological characteristics of the Karak watershed.


**Figure 2.** Regional geology map of the study area illustrates main structural and lithological units.

#### *2.3. Hydrological Background*

In the study area, the estimated thickness of semi-confined aquifers ranges from 10 to 30 m. The groundwater quality in the northeastern part of the northwest catchment is inferior [42]. This situation occurs due to the salt rock in the northern mountainous region, which is dissolved by runoff water and polluted groundwater due to deep infiltration. Under diving conditions, groundwater flows through weathered layers and fault zones. The alluvial filling is very uneven and contains high level of silt and clay. Locally, sand and gravel beds were encountered in boreholes. The flow rate through the open wells is calculated to be 0.035 mm3/year. The alluvial aquifer's average annual recharge is approximately equal to the average annual discharge, which is 2.7 mm3/year. The groundwater level is between 29.03 and 238.66 m. This indicates that a fuzzy groundwater boundary exists corresponding to a surface water boundary [49]. A small dam (Chambia dam) was constructed to maintain the groundwater level in the Karak watershed. The soil texture of the study area is predominantly medium clay, pure sand, cultivable soil and crops.

#### **3. GCF Analysis and Optimization**

The evaluation of groundwater condition factors (GCFs) is essential to effectively determine an accurate groundwater potential (GWpot) index [50]. GCFs should be considered in terms of regional topographical, geological, hydrological, and land use/land cover characteristics influencing the GWpot [15]. Therefore, the identification of the GWpot spatial distribution was performed by multi-criteria decision-making (MCDM) analysis

of nine factors, i.e., drainage density, geology, lineament density, slope, soil type, rainfall, elevation, land use/land cover, and groundwater level fluctuations. These GCFs were extracted independently from appropriate remote sensing, geological, and conventional map datasets (Table 2).

**Table 2.** GCFs used for mapping groundwater potential of the Karak watershed.


The drainage density (Dd) is a measure of the total length of all streams per unit area, regardless of the stream networks [51]. The hydrology toolkit in ArcGIS 10.4 was used to extract stream networks from a digital elevation model (DEM). Accordingly, Dd was calculated as the stream's total length divided by the total drainage using Equation (1) [14]. Subsequently, the drainage frequency was classified into five categories using a natural break classification scheme [16]. High drainage frequency is associated with high permeable lithology and accordingly high GWpot. The groundwater favorability is indirectly related to Dd, which is related to surface runoff and permeability [52].

$$\text{DD} = \sum\_{l=o}^{n} \frac{D\_l(\text{km})}{A(\text{km})^2} \left(\text{km}^{-1}\right) \tag{1}$$

where DD represents drainage density, *Dl* is the stream's length, and *A* is the watershed area (km2).

Lineaments are surface manifestations of linear or curvilinear features, such as joints, straight streams, and regional vegetation placement, reflecting potential topographical or geological structure [15]. The seven bands of the Landsat 8 image were stacked using ENVI 4.8 (Harris Geospatial, Broomfield, CO, USA), and principal component analysis (PCA) was performed on the stacked image in QGIS (Open Source Geospatial Foundation, Bern and Chur, Switzerland). The thematic layer for Ld can be defined as the total length of all recorded lineaments divided by the catchment area under consideration, as shown in Equation (2) [53]. The higher the Ld, the higher the favorability of GWpot.

$$\text{LD} = \sum\_{i=0}^{n} \frac{L\_i \text{ (km)}}{A \text{ (sq. km)}} \left(\text{km}^{-1}\right) \tag{2}$$

where LD represent the lineament density, *Li* is the lineament's length in km, and *A* is the grid area in square kilometer.

Data from 17 metrological stations were processed using simple arithmetic mean, isometric, and Thiessen polygon interpolation methods to obtain sufficient uniform precipitation in the catchment area. After these three interpolation methods were used for comparison, isometric interpolation (Equation (3)) was considered the best technical interpolation method. The flat and gentle areas, with less runoff, are more favorable for GWpot than steep slopes [54]. In addition to rainfall quantity, other precipitation characteristics (such as duration and intensity) are also important. For example, a 20 mm rainfall in a long period may have a more significant impact on groundwater recharge than a 50 mm rainfall in a short period.

$$P = \frac{\sum\_{l=1}^{n} p\_l}{N} \tag{3}$$

where P is the average precipitation depth, with *p*1, *p*2, *p*3 up to *pn* being the rainfall records of measurement stations 1, 2, 3, up to *n*, respectively.

The slope is an important factor that directly controls the infiltration of surface water. A 30-m resolution DEM was processed to generate a slope map in the ArcGIS 10.4 spatial analyst toolkit. The slope gradient was reclassified into five classes using the quantile classification scheme presented by [18]. A higher slope is more conducive to runoff but has a smaller impact on groundwater recharge. Elevation or altitude can have an indirectly inverse effect on GWpot, which relates primarily to the occurrence of rainfall and the resulting recharging. However, high altitudes favor more recharge and ensure groundwater availability in low land areas in a watershed. Mountainous regions are often favorable for the recharge of deep-seated confined aquifers situated at low land areas [55].

Stratum lithology influences the porosity and permeability of aquifers and directly affects the GWpot. The porosity of rocks, alluvial/sedimentary layers, sand, silt, and clay beds determine water infiltration and percolation [56]. Therefore, the lithology factor was also considered concerning groundwater characteristics. The lithology map was extracted, digitized, and reclassified from the geological map of Northern Pakistan. Accordingly, different weights were assigned to rock units depending on the infiltration capacity and GWpot as per multiple influence factor criterion.

Vegetation cover areas, such as forests and agriculture traps, retain water by the roots of plants. In contrast, the built-up and rocky land cover decreases groundwater recharge by increasing the runoff during rainfall [24]. Therefore, to conduct GWpot studies, it is necessary to investigate the land use land cover (LU/LC) characteristics of the study area. Therefore, the LU/LC map from the Forest Management Center Peshawar (FCMP) was reclassified with different score values assigned to several subclasses.

The water retention capacity of an area depends on the type of soil and its permeability. Permeability is directly related to the soil effective porosity which is greatly influenced by the particle shape, size, adsorbed water, porosity, saturation, and the presence of impurities in the soil [57]. The soil type map was primarily derived from the Directorate General Soil and Water Conservation (DGSC), KPK, and updated through onsite inspections. Soil mainly influences infiltration and percolation processes that eventually affect the groundwater recharge and then the GWpot of a given area.

#### **4. Methods**

In this study, the application of remote sensing (RS) and geographic information system (GIS)-based spatial data and geoelectric data assisted hydrogeological assessment to distinguish the sediments and rock units of groundwater significance. The flowchart developed in this study is shown in Figure 3, which contains four steps:


of the resistivity mapping, overburden thickness mapping, and reflection coefficient mapping.

4. To evaluate the accuracy of GWpot mapping, the performance of the MIF model is assessed based on groundwater level (GWL) data through a confusion matrix, Kappa (K) analysis, and a Receiver Operating Characteristics (ROC) curve. In addition, the quantitative aquifer potential interpreted by VES data indirectly verifies the MIF model's predictive performance. Meanwhile, the hydrologic stratigraphic prediction derived from ERT and VES numerical models is correlated with known boreholes lithological information.

**Figure 3.** Framework to delineate groundwater potential and to identify hydrogeological characteristics.

#### *4.1. Multi-Influence Factors (MIF) Model*

4.1.1. Assigning of Weights and Ranks

The GWpot index is influenced by several hydrological, geological, topographical, environmental, and climatic variables [2]. By means of GCFs analysis and optimization, geology (GEO), lineament density (LD), drainage density (DD), slope (SL), soil type (ST), rainfall (RF), elevation (EL), land use/land cover (LULC), and groundwater level fluctuations (GLF) were identified as the input data of the MIF model. The MIF model involves drawing a graph with correlations between conditioning factors and assigning weights based on the strength of the interrelationships (Figure 4) [2]. In Figure 4, a continuous arrow shows a major influence, and a dashed arrow indicates a minor influence on the other GCFs. The weights and ranks were assigned to each GCFs and different classes based on their relative contribution to GWpot using the heuristic approaches/knowledge-driven method [11,58,59]. Weights of 1.0 and 0.5 were allocated to each major and minor effective variable, respectively. The combined weights of both major (*CFh*) and minor (*CFl*) were considered for calculating the comparative ranks (Table 3). Since the estimated weight of each GCF is equally distributed and applied to each GCF' category, the final GWpot map is a weighted average. The estimated weight for each conditioning factor was obtained as a percentage using Equation (4).

$$\text{Score} = \left[ \frac{(\text{CF}\_{l} + \text{CF}\_{l})}{\sum \left( \text{CF}\_{h} + \text{CF}\_{l} \right)} \times 100 \right] \tag{4}$$

where, *CFh* is the major weight of the condition factor, and *CFl* is the minor weight.

**Figure 4.** Interrelationship between the GCFs concerning the GWpot index.


#### 4.1.2. Weighted Overlay Analysis (WOA)

The GWpot index quality is influenced by the quality and quantity of the input data and the predictive models used [2]. Weighted overlay analysis [60,61] in ArcGIS 10.4 (Environmental System Research Institute, Redlands, California, United States) was used to outline the spatial distribution of the groundwater potential index based on nine GCFs' superimposition and their corresponding percentage effects on the groundwater potential. This work was done by multiplying each factor's category cell value by the factor's weight and summing the resulting cell values to generate a GWpot map, as summarized in Equation (5). A GWpot index is a calculated dimensionless number considering the weight assigned for each GCF and its categories [3]. After the WOA analysis had been completed, the natural break method was used to categorized GWpot into four levels of potentiality (i.e., low, medium, high, and very high).

$$\begin{aligned} \text{GW}\_{\text{potz}} &= \sum\_{i=0}^{n} \text{W}\_{i} \times \text{R}\_{i} \\ &= \text{DD}\_{\text{\text{\textdegree}}} \text{DD}\_{\text{W}} + \text{LD}\_{\text{\textdegree}} \text{LD}\_{\text{W}} + \text{RF}\_{\text{\textdegree}} \text{RF}\_{\text{W}} + \text{SL}\_{\text{\textdegree}} \text{SL}\_{\text{W}} + \text{EL}\_{\text{\textdegree}} \text{EL}\_{\text{W}} + \text{GEO}\_{\text{\textdegree}} \text{GEO}\_{\text{W}} + \\ &\quad \text{LLUL}\_{\text{\textdegree}} \text{LLLC}\_{\text{W}} + \text{ST}\_{\text{\textdegree}} \text{ST}\_{\text{W}} + \text{GLF}\_{\text{\textdegree}} \text{GLF}\_{\text{W}} \end{aligned}$$

where GWpotz is the groundwater potential index, Wi is the weight of each condition factor, Ri is the rank of each GCF's category, DD is the drainage density, LD is the lineament density, RF is the rainfall, SL is the slope variation, EL is the elevation, GEO is the lithology, LULC is land-use/land-cover, ST is the slope type, and GLF is the groundwater level fluctuation. The subscripts c and w indicate a category of a GCF's thematic layer and its corresponding percent influence on GWpot, respectively. This overlay analysis was done by multiplying the rank of each GCF's category (each individual category has a rank) with the weight of each condition factor (each GCF has a unique weight) to obtain the GWpot index at the corresponding position of GCFs.

#### *4.2. Accuracy Assessment of the MIF Model*

The pre-monsoon and post-monsoon groundwater table (GWT) data from 32 observed boreholes with global positioning system (GPS) positions were collected for validation purposes. The area under the curve (AUC) based receiver operating characteristic (ROC) curve, the confusion matrix, and Kappa (K) analysis were used to test the performance of the MIF model. The ROC is a mathematical technique developed to explain the efficiency of probabilistic deterministic detection and prediction systems [62,63]. In this study, ROC was used to assess the spatial consistency between real events and to predict the model probability. In the validation phase, pre-monsoon and post-monsoon GWT data of 32 observed boreholes/tube wells were compared with the GWpot result obtained by the MIF model. The ROC curve provides a quantitative evaluation that can determine the uncertainty of modeling and evaluate the spatial model effectiveness. The confusion matrix and Kappa (K) analysis [26] were also used for accuracy evaluation by correlating the GWpot map with the observed GWT data. The overall accuracy was calculated using the following formula [64].

$$\text{OA} = \frac{\sum \text{C}\_{OWL}}{\sum \text{OWL}} \tag{6}$$

where, OA is the overall accuracy, *COWL* represent the number of correct observation boreholes/well's locations and *OWL* is the number of observation boreholes/well's locations.

The Kappa (K) analysis is a multivariate approach for MIF accuracy evaluation. It was calculated by the following formula [65].

$$\mathbf{K} = \begin{bmatrix} \frac{\sum \text{ CV } \%-\text{CAOV}\%}{\sum \text{ TC} - \text{CAOV}\%} \end{bmatrix} \tag{7}$$

where, CV % is the percentage of the correct values, CAOV% is the percentage of the correct agreement to observed values, TC is the total number of class.

#### *4.3. Interpretation of Geoelectrical Data*

The geophysical techniques have typically been used to assess hydrogeological structures, hydro-stratigraphic characteristics and the spatial distribution of aquifers [34]. Fundamentally, an electrical current is injected into the ground by two current electrodes and measures the potential difference between the other two pairs of electrodes. In this study, two-dimensional electrical resistivity tomography (ERT) based on dipole–dipole configuration and vertical electrical sounding (VES) based on Schlumberger configuration measurements were performed using essential field equipment (Terameter SAS 100 and SAS 1000 Lund imaging systems and their accessories, ABEM, Sundbyberg, Sweden) (Figure 5).

**Figure 5.** Schematic diagram of (**a**) the Schlumberger array configuration for vertical electrical sounding (VES), and (**b**) dipole–dipole array configuration for electrical resistivity tomography (ERT) techniques.

#### 4.3.1. Electrical Resistivity Tomography (ERT)

The ERT technique was effectively applied in the surveyed area to provide information about subsurface hydrogeological characteristics to fully understand the GWpot and hydrostratigraphy through vertical and horizontal two-dimensional sections capable of reaching lengths and depths up to 176 m and 30.2 m, respectively. A multi-electrode 2D device (Terameter SAS 100) along a dipole–dipole configuration including electrodes connected to a transmitter/receiver system via a multi-core cable was used to acquire data (Figure 5b). The dipole–dipole configuration exhibits an excellent vertical and horizontal resolution of subsurface geological features, which has great horizontal coverage and penetration depth [66]. The apparent resistivity was calculated for every electrode quadrupole by Equation (8) [34].

$$
\rho = K \frac{V}{I} \tag{8}
$$

where *V* is the voltage, *I* is the current, and *K* is a geometric factor.

The dipole–dipole configuration data were concatenated to obtain combined apparent resistivity pseudo-sections. The degree of consistency between resistivity and actual subsurface resistivity distribution depends on the combination of acquisition parameters and inversion strategy. The smoothness constrained least-squares technique in the RES2DINV (Landviser, League, Texas, United States) program was used to process the apparent resistivity data [67,68]. This process automatically creates 2D models in a rectangular block by selecting the optimal data inversion parameters (e.g., the damping coefficient, and the vertical and horizontal flatness filter ratio, convergence limit, number of iterations). We used the finite difference method to calculate the module's apparent resistivity and compared it to the measured data. Iteratively, we adjusted the resistivity of the model block until the calculated apparent resistivity value of the model matched the actual measurement. Finally, the program produces a pseudo section (a qualitative method for measuring or calculating resistivity changes) and an inverse model section (slice depth and resistivity tomography image) [68]. As a follow-up to the observation results of ERT lines L1, L2, L3, L4, L5, and L6 were acquired in the E-W, S-W, N-E, E-W, E-W, and E-W directions, respectively. In this study, the ERT technique estimated the spatial subsurface resistivity

caused by the lateral and longitudinal inhomogeneities of petrophysical properties. The distribution of petrophysical properties indicated by the measured resistivity fluctuations were generated to guide GWpot and hydro-stratigraphy in the study area.

#### 4.3.2. Vertical Electrical Sounding (VES)

The VES method was used in the surveyed area to evaluate the hydro-stratigraphic structure of the sedimental layer (i.e., the structure of the subsurface sediments), aquifer characteristics (e.g., thickness, resistivity (ρ), overburden thickness, and reflection coefficient), and GWpot. The VES technique is one of the most commonly used conventional resistivity methods to determine the vertical variation of subsurface resistivity parameters [34]. In the surveyed area, 26 VES measurement stations were operated at different positions using the Schlumberger electrode configuration with half-current electrode spacing (AB / 2) ranging from 1.5 to 1000 m in each successive electrode probe to determine the depth to the sediments and apparent resistivity (ρa). Meanwhile, using the Schlumberger array (Figure 5a), the adequate penetration depth is typically 20–40% of the external electrode spacing (AB), depending on the subsurface resistivity structure [69]. In this study, we first plotted all resistivity data collected to confirm qualitive and qualitative characteristics. The statistical apparent resistivity (ρa) values of the Schlumberger array for each sounding were calculated using Equation (9).

$$\rho\_{\rm a} = \pi \left\{ \frac{\left(\frac{\rm AB}{2}\right)^2 - \left(\frac{\rm MN}{2}\right)^2}{\rm MN} \right\} \text{Ra} \tag{9}$$

where, AB represent the distance between two current electrodes, MN is the distance between the potential electrode, and Ra is the apparent electrical resistance.

The preliminary interpretation was performed using Partial Curve Matching (PCM) and auxiliary tools to summarize VES values, i.e., the relationship between the apparent resistivity and corresponding half current electrode spacing (AB/2) on the double logarithmic graph. The results obtained from the exercises were used as an input model for computerassisted iterations using the WinResist™ (Geotomo Software, Gelugor, Penang, Malaysia) program. The preliminary interpretation of VES data was quantitative, determining the thickness (h) and resistivity (ρ) of different layers, and qualitative inferring lithology was based on the resistivity and reflection coefficient (RC) values of each sounding station. For better depiction, six VES measurements were performed in the two boreholes' immediate vicinity (BH06/BH09) and correlated with known lithological information. The Schlumberger configuration was characterized by tracking and tracing each VES subsurface layer, the vertical changes, and the geoelectric profile with a known borehole/well lithology to horizontally correlate the measured VES to perceive a unified layer model applicable to all field curves. Moreover, geological information of known borehole/wells can improve interpretations that lead to lithological results from VES data, while software analysis can only provide resistivity distinction by depth. The statistical apparent resistivity values of each VES measurements were outlined to create an iso-resistivity map. The RC values for the surveyed area were calculated using the following expression [70].

$$\text{RC} = \left\{ \frac{(\rho n - \rho(n-1))}{(\rho n + \rho(n-1))} \right\} \tag{10}$$

where RC represents the reflection coefficient, *ρn* is the resistivity of the *n*-th layer, and *ρ*(*n* − 1) is the resistivity overlying the *n*-th layer.

#### *4.4. Geophysical Well Logging*

Hydrogeological characterization of aquifers using geophysical well/borehole logs has been emphasized in many studies [5,71]. Effective groundwater exploration and well/borehole lithology evaluation require a complete understanding of aquifer hydrogeological characteristics and well/borehole design. In the study area, the drilling sites were selected based on the experience of the MIF model to determine prerequisites for the successful construction of the tube well and evaluate the availability of groundwater supply that can meet the demand for domestic and irrigation water. The GeoLog International (GLI) groundwater and engineering services with reference to Ms. Manahil Engineering & Cons conducted St. Rotary (SR) drilling and geophysical logging in Marwatan Banda, Karak. The borehole's logging survey was conducted using multi-parameter methods, i.e., normal resistivity logs (NRLs) (short and long configuration) and spontaneous potential logs (SPLs). The Geo logger 3030/Mark-2 3433 (GLI, Peshawar, Pakistan) was used for petrophysical property measurements. Through these significant hydrogeological properties, e.g., the formations' lithology, depth, thickness, groundwater water table level, and groundwater quality in total dissolved solids (TDS) were evaluated.

#### **5. Results**

#### *5.1. Evaluation of GCFs*

The MIF model is an MCDM technique widely used for environmental management and has proven to effectively explain the GWpot influential factors. It can effectively determine GCF weights. Table 4 illustrates the weights and qualitative ranks assigned to each influencing factor described below.

Drainage density (Dd) is a measure of the total length of all streams per unit area, regardless of the stream networks [51]. Subsequently, the drainage frequency was classified into five categories, i.e., very low (1.08–1.61 km/km2), low (1.61–1.86 km/km2), moderate (1.86–2.11 km/km2), high (2.11–2.38 km/km2), and very high (2.38–3.08 km/km2) (Figure 6a), according to a natural break classification scheme. The groundwater favorability is indirectly related to drainage density, as are surface runoff and permeability. Therefore, the highest score was assigned to the 1.08–1.61 km/km2 category, indicating high infiltration and low runoff, and the lowest score was assigned to the 2.38–3.08 km/km<sup>2</sup> category (Table 4).

Lineament density (Ld) of the Karak watershed indirectly indicates the GWpot, as the presence of lineaments usually means a porous zone. The lineaments are spatial distributed in the study area aligned in the directions of E-SW, NNE-SSW, NW-SE, and E-W, and their density was classified into five frequency categories (Figure 6b). The higher the Ld, the higher the probability of GWpot. Therefore, the highest rank was assigned to the 1.46–1.78 km/km2 category and the lowest was assigned to the 0.17–0.45 km/km<sup>2</sup> category.

Rainfall (RF) interpolated data were reclassified into five categories, i.e., very low (13–281 mm), low (282–577 mm), moderate (578–604 mm), high (605–629 mm), and very high (630–663 mm) (Figure 6c). In addition to the quantity of RF, other precipitation characteristics, such as duration and intensity, are also important. For example, a long period of 20 mm RF has a more significant impact on groundwater recharge than a short period of 50 mm RF.

The slope (SL) map was reclassified into five categories, i.e., flat (0–5.78◦), gentle (5.78◦–13.5◦), moderate (13.5◦–23.1◦), steep (23.1◦–35.0◦), and very steep (35.0◦–81.9◦) using the quantile classification scheme presented in [18]. The flat and gentle areas are more suitable for GWpot than steep slopes, as a gentle and flat slope allows for less runoff, and a steep slope is more conducive to runoff [54]. The highest rank was assigned to flat area (0–5◦.780◦), and the lowest was assigned to a very steep area (35.0◦–81.9◦), which has a smaller impact on recharge in the study area (Figure 6d).


**Table 4.** Classification of weight and ranks of GCFs.

The elevation (EL) map in Figure 6e shows three elevation categories, i.e., high (707–1419 m), moderate (304–706 m), and low (0–303 m).

Geology (GEO) characteristics govern the porosity and permeability of the hydrogeological layer, which in turn influences the formation and distribution of GWpot through physio-mechanical properties that control the water transmitting ability of the hydrogeological layer materials and the rate of groundwater flows. Therefore, the GEO factor was also considered concerning groundwater characteristics. The study area consisted of eight lithological units of formation types and geological ages. The confirmed lithology outcrops are the Quaternary alluvium (Q), Dhok Pathan formation (DP Fm), Chinji formation (C Fm), Jurassic-Triassic rocks (J/T), Kohat formation (K Fm), Nagri formation (N Fm), Kamlial formation (K Fm), and Drazinda shale (DS) (Figure 6f).

Land use/land cover (LULC) greatly influences groundwater occurrence and exploitation. The major portion of the study area is agriculture (62%; 1345 km2), followed by forest area (15%; 576 km2), barren land (12%; 292km2), rangeland (4%; 58.3 km2), shrubland (3%;

40.7 km2), built-up (2%; 30 km2), river/stream (1%; 22 km2), and dam/pond (1%; 8 km2) (Figure 6g).

Soil type (ST) and its permeability decides the water retention capacity of an area. The soil types of the study area include loamy soil, loamy clay, and mainly loamy (Figure 6h). The dominant soil type in the study area is loamy soil. The coverage of the other two soil types (i.e., loam and mainly loam) are relatively low. According to composition and soil water holding capacity, the loam is regarded as the highest grade, and mainly loam is regarded as the lowest grade.

Groundwater level fluctuation (GLF) is of significance in the successful management of GWpot. Pre-monsoon and post-monsoon groundwater levels (GWLs) indicate the degree of saturation and the extent of recharge aquifers. In this study, hydrogeological data of 32 boreholes/wells over 10 years 2009–2019 (from Pakistan Water and Power Development Authority (WAPDA)) was collected through onsite investigation. During the period 2009–2019, the pre-monsoon and post-monsoon water level varies from 5.9 to 15.4 mbgl and from 7.3 to 32.6 mbgl, respectively (Figure 6i). The groundwater fluctuation levels were calculated for the period of 2009–2019, with a minimum of 1.57 m and a maximum of 19.3 m. In the study area, the aquifer is partially saturated due to the inadequate precipitation and other influencing factors. In the northern region, slight fluctuations of groundwater level (about 6 m) were observed, which may have been due to groundwater recharge by surface irrigation. However, groundwater levels fluctuated significantly in the southern and central regions, which may have been caused by topographical influence and the excessive exploitation of groundwater.

#### *5.2. Assessment of GWpot*

Using the weighted overlay analysis in the GIS environment, the GWpot zones were evaluated by integrating several conditioning factors (i.e., rainfall, slope, geology, soil type, drainage density, lineament density, land use/cover, elevation and groundwater fluctuation). Based on natural breaks in the histogram of the GWpot index, the GWpot map was categorized into four levels of potentiality, i.e., low, medium, high, and very high (Figure 7a), with the distribution ranges of 9.7% (72.3 km2), 52.4% (1307.7 km2), 31.3% (913.4 km2), and 6.6% (44.8 km2) of the total area, respectively (Figure 7b,c). The spatial distribution of the various GWpot zones typically shows a mirror reflection of key factors. High and very high GWpot zones confirm their excellent capacities as sedimentary groundwater aquifers. The GWpot map demonstrates that the excellent groundwater is concentrated due to the distribution of Quaternary alluvial and agricultural land with high infiltration ability. Moreover, high drainage densities and low slope gradients can increase groundwater infiltration capacity, which may be related to the evaluated high GWpot. The northwestern, southeastern, and the central part limited regions typically have a low to medium GWpot, accounting for approximately 12.7% of the study area.

**Figure 6.** GCFs considered in this study: (**a**) drainage density; (**b**) lineament density; (**c**) rainfall; (**d**) slope; (**e**) soil type; (**f**) land use land cover; (**g**) geology/lithology; (**h**) elevation; (**i**) groundwater level fluctuations.

**Figure 7.** (**a**) GWpot zones and groundwater level depths of boreholes/wells; (**b**) groundwater potentiality in square kilometers and (**c**) in percentage of the Karak watershed.

#### *5.3. MIF Model's Performance*

The ROC curve, the confusion matrix, and Kappa (K) analysis were used to evaluate the accuracy of the assessment result and the performance of MIF the model.

ROC graphs are useful tools for visualizing a classifier's performance and for determining the area under the curve (AUC) value to evaluate an algorithm [62]. The ROC curves were implemented in the present study as a goodness of fit, and the success rate can be distinctly visualized. In this study, the predicted GWpot map was examined and compared with 32 pre-monsoon and post-monsoon groundwater level (GWL) fluctuations to evaluate the spatial coincidence between the favorability values (from GWpot) and the actual GWL fluctuation events (Figure 8a). The GWL fluctuations range from 1.57 to 19.3 m (Figure 8b). Since a larger area under the ROC curve indicates that the spatial GWpot mapping is more effective, an AUC value of 1 shows a perfect prediction of the model and indicates that the highest ranked probabilities coincide with the groundwater fluctuation [63]. The result of the ROC chart analysis shows that the AUC value of the presented MIF performance is 68% (Figure 8c) which is consistent with GWL fluctuation.

**Figure 8.** (**a**) The pre-monsoon and post-monsoon groundwater level (mbgl) fluctuations; (**b**) average groundwater level fluctuation (m) of the Karak watershed; (**c**) receiver operating characteristics (ROC) curve of the MIF model.

The confusion matrix and Kappa (K) analysis were performed using the 32 actual groundwater depths from boreholes/wells. The groundwater depth in the study area is between 6.7 and 190 m. These 32 depths were divided into four categories, i.e., 6.7–36 m, 36–80.76 m, 80.76–130 m, and 130–190 m. The groundwater depth data were used to calculate classification accuracy by the confusion matrix and Kappa (K) analysis. Overlay analysis shows that most of the boreholes/wells with higher groundwater levels are located in areas with demarcated higher groundwater potential. The performance evaluation of the MIF model shows that the overall accuracy is 68%, and the Kappa coefficient is 0.65 or 65% (Table 5), which indicates that the estimated potential of groundwater is consistent with the investigated groundwater depths in the study area.


**Table 5.** Error matrix of the GWpot zone-based confusion matrix and Kappa (K) analysis.

<sup>1</sup> CS refer to the correct sample.

#### *5.4. ERT Interpretation*

In this study, the ERT approach with an optimal compromise between the electrode distance and profile length produced a deep characterization of the hydro-stratigraphical layers and groundwater potentiality. The smoothness constrained least-squares outputs by the RES2DINV software show an apparent lateral homogeneity with a gradual increase in resistivity, with depth caused by lateral and longitudinal inhomogeneities of rock physical properties (Figure 9a). Each inversed resistivity section obtained a distribution of petrophysical properties of resistivity variability and possible resistivity anomalies (which may be water-bearing zones). The final depth of the inversed sections ranges from 5 to 30.2 m.

**Figure 9.** (**a**) Inverse model resistivity section of ERT survey lines (i.e., L1, L2, L3, L4, L5, and L6) containing existing borehole lithological information on L3; (**b**) correlation of acquired ERT hydrologic stratigraphy with (**c**) The existing borehole lithological information; and (**d**) ERT measurements line alignment in the study area, in which the red dot shows the location of the existing borehole, and the yellow dots indicate the proposed locations of wells for groundwater exploitation.

Generally, the root means square (RMS) error at the end of eight iterations of almost every ERT section is less than 8%. The interpretation of ERT sections is based on a standard resistivity range of values. The recommended GWpot zones were based on an understanding of the subsurface sediment/ rock lithology of the study area. Meanwhile, the subsurface lithology related to the resistivity range was derived from the existing standard resistivity chart, which considers other local factors that may cause the resistivity deviation.

In the study area, the L1, L2, L3, L4, L5, and L6 ERT inversed resistivity values area 12.8–189 Ωm, 12.8–189 Ωm, 3.62–792 Ωm, 12.8–189 Ωm, 3.62–792 Ωm, and 3.62–792 Ωm, respectively (Figure 9a). The inverse resistivity models using dipole–dipole configurations on L2, L4, and L6 ultimately revealed the vertical and lateral distribution of subsurface resistivity. According to the predicated GWpot on diffusion and array configuration, the

groundwater prospect resistivity values are 12.8–48.3 Ωm, 14–76.4 Ωm, and 3.62–16.3 Ωm on L2, L4, and L6, respectively. Variation of resistivity characteristics within the primary lithological unit ultimately indicates the GWpot prospect adjacent to clayey sand and silicate aquifers (sandstone) (Figure 9a). This result is consistent with the Karak watershed regional geology, which is mainly composed of interlayers of fine sand, sandstone, clay, and gravel. Since the GWpot is structurally controlled, it also needs to locate potential fracture zones, e.g., fractured sandstone, which are considered good aquifer sources. The ERT techniques should be applied with a proper understanding of the hydrogeological background. Therefore, five lithological sequences (i.e., topsoil with coarse gravel and sand, silty sand mixed lithology, clayey sand/fine sand, fine sand/gravel, and clayey basement) of the drilled borehole on L3 at final depths of 45 m were normalized with the ERT model by mean of quantitative quota (Figure 9b). The ERT-predicted hydro-stratigraphy and borehole lithological log signature (Figure 9c) performance analysis shows suitable matches. The marked yellow points on the L3, L4, and L6 sections are considered future prospects for groundwater exploitation (Figure 9d). These high groundwater potential zones will play a vital role in the future expansion of drinking water and irrigation development in the surveyed area.

#### *5.5. VES Interpretation*

#### 5.5.1. Hydrogeological Characteristics

The VES technique has been proven efficient in evaluating hydrogeology, aquifer properties, and aquifer potential. In this study, aquifer characteristics (such as thickness, lithology, and resistivity, reflection coefficient, and isopach) were determined, which is an essential factor in hydro-stratigraphic inheritance and GWpot assessment. The apparent resistivity data obtained from the VES positions were plotted against half of the current electrode spacing (AB/2), and a curve matching technique was used to interpret resistivity sounding curves (Figure 10). This technique involves matching small segments of the field curve against the trendline curve to determine the thickness of a particular layer in half-space and the apparent resistivity. As far as the evaluation of the statistical apparent resistivity is concerned, the qualitative interpretation results indicating that the curves, stratification properties, and RMS errors are in complete agreement (Figure 11) (Appendix A). Depending on the shape of the VES curve, the resistivity distributions of various hydro-stratigraphy can be classified into H, K, A, and Q types, which can be mutually combined to generate HA, HK, KH, and QH types [72]. In this study, the type of curves observed include 3-layer H-type (26%), 4-layer HA-type (9%) and KH (52%), and 5-layer HKH-type (13%). Qualitative hydrological inferences can usually be based on the type of curve.

The geoelectrical interpretation based on curve matching reveals hydrologic resistance and depth variation (Figure 10). According to the corresponding resistivity values (ρ1, ρ2, ρ3, ρ4, and ρ5) and thicknesses (h1, h2, h3, h4, and h5), the geoelectric units indicate four to six sequences of lithologies, i.e., topsoil (coarse gravel and sand), alluvial layer, silty sand, clayey sand, fine sand and gravels, and clayey sand with saline water. Table 6 summarizes the VES interpretation, including the number of hydrologic layers and their corresponding resistivity values and the inferred lithology information. Appendix A presents the detailed explanation of geoelectrical stratification for all the VES surveys carried out in the Karak watershed and the resistivity variation.

**Figure 10.** VES data interpretation result-based partial curve matching (PCM) along 26 VES stations: (**a**) hydrologic layers depth variation; (**b**) hydrologic layers resistivity variation.

**Figure 11.** Root mean square (RMS) error of 26 VES stations (**left**) and the reflection coefficient variation of VES stations (**right**) in the study area.

**Table 6.** Average inferred hydro-stratigraphy corresponding to resistivity in the study area. The detailed VES interpretation results are shown in Appendix A.


5.5.2. VES Correlation with Boreholes

For better delineation of the hydro-stratigraphy, six VES results adjacent to two boreholes (BH06/BH09) were correlated with known lithological information. Performance analysis shows that VES1 yields five lithological units (coarse gravel/sand, silty sand mixed lithology, clayey sand, and fine sand/gravel) (Figure 12). The zone of interest with water saturation lies at a depth of 29.8 m. VES3 penetrates up to 48.1 m where water is predominantly saline, with freshwater saturation having a lithology of coarse gravel/sand, silty sand mixed lithology, silty sand/gravels, fine sand/gravel, and clayey sand/saline water. VES2 yields three lithological units, where the zone of interest lies at a shallow depth. Furthermore, VES9, VES17, and VES8 correlated with borehole BH09 show suitable matches, where salinity and freshwater saturation are encountered at a shallow depth (25 to 35 m) due to capillary action. However, the VES8 upper portion is mainly composed of unconsolidated alluvium, and the freshwater zone is at a shallow depth due to elevation. The main lithological characteristics of the topsoil at each VES station are predominantly alluvium. The VES and borehole log signature performance analysis show suitable matches between them (Figure 12).

**Figure 12.** Correlation of VES data interpretation results with borehole lithological information.

#### 5.5.3. GWpot Based on VES

Aiming at monitoring aquifer potential, a preliminary conceptualization of geoelectrical properties governing the reflection coefficient, the aquifer's overburden thickness, and resistivity is needed during VES measurements. These basic and essential interpretative criteria are described below.

The reflection coefficient (RC) is an essential geoelectric factor, as it helps to identify the permeable hydrologic layers carrying the GWpot. The RC values of the VES positions in the surveyed area were calculated using Loke's method [72]. Figure 13a shows the changes in RC values detected by each VES station. Differences in subsurface resistivity and lithology cause the RC fluctuations. The calculated RC values were contoured in Surfer 15 software, and an RC map shows a value range of 0.50–0.95 (Figure 14a). Olayinka [73] observed that the subsurface topography usually shows a good aquifer when the overburden is relatively thick and/or the reflection coefficient is low (<0.8). RC mapping has been found to be useful in investigating the hydrogeological aquifer because it reveals whether a permeable aquifer is filled with water. Therefore, an anisotropy coefficient for this parameter was considered in this study.

**Figure 13.** (**a**) Reflection coefficient (RC); (**b**) overburden thickness along 26 VES stations in the surveyed area.

**Figure 14.** (**a**) Reflection coefficient map; (**b**) overburden thickness map; (**c**) apparent resistivity map based on the interpretation of VES data.

An overburden thickness/isopach map was plotted and contoured according to the interpreted depths to the sedimentary rock (Figure 14b). The isopach map illustrates the thickness variation in a hydro-stratigraphic layer, a tabular unit, or a stratum [29]. Isopach mapping is essential in the hydrogeological investigation because it shows the number of hydrogeologic layers above the aquifer, and where groundwater can be observed in areas considering the overburden thickness. The overburden thickness variation of the aquifer along VES can be seen in Figure 13b. The overburden thickness in the surveyed area varies between 6.3 and 65.6 m. The isopach map shows that the overburden thickness in the northern, eastern, and southern parts of the surveyed area ranges from 20 to 50 m (Figure 14b). In contrast, the relatively thin overburden thickness of 5–15 m is virtually around the central and western parts of the surveyed area. The overburden thickness is shallow in most probing stations, indicating that the basement is close to the surface. Therefore, groundwater in these areas is highly dependent on the occurrence of fractures [29].

The apparent resistivity values of all VES stations were contoured to produce an iso-resistivity map (Figure 14c), indicating that the apparent resistivity increases radially outward from the center of the region and the resistivity values are 10–1150 Ωm. The resistivity of the bedrock represents the resistivity of the deepest hydrological layer in the surveyed area. It has been found that the resistivity of the bedrock is of significance in many aspects of hydrogeological and hydro-geophysical measurements because it plays a vital role in assessing the potential of groundwater. After all, the resistivity of the bedrock has the potential to reveal fractured aquifers.

The lower RC and relatively high overburden thickness can increase a well's groundwater productivity [74]. In this study, the considered GWpot geoelectrical factors includes reflection coefficient, overburden thickness, and iso-resistivity obtained from the interpretation of VES data. This quantification of aquifer potential indirectly verified the accuracy of the MIF model and its predictive performance. The VES stations in the surveyed area were divided into high yield, medium yield, and low yield groundwater by employing Olayinka's basic criteria [73].


Considering these criteria, the RC and overburden thickness were used to produce the parameters for categorizing VES stations by the GWpot, i.e., VES1, VES5, VES6, VES8, VES9, VES14, VES16, VES17, VES21, VES24, and VES26 have high yield GWpot (Figure 15a), VES3, VES7, VES11, VES13, VES19, VES20, and VES25 have medium yield GWpot (Figure 15b), and VES2, VES4, VES10, VES12, VES15, VES18, and VES22 have low yield GWpot (Figure 15c). Based on these groundwater potentiality variations among the VES stations, a final GWpot contour map of the surveyed area was generated, and it demonstrates that the northern, northeastern and eastern parts have excellent GWpot for future exploitation and development, while the low and medium GWpot regions are located in the western and central parts of the surveyed area (Figure 16). The VES-based groundwater potential map was compared with the groundwater potential map obtained by the RS and GIS-based MIF method. This indicated that the MIF method is accurate and consistent in predicting GWpot.

**Figure 15.** Groundwater potential VES distribution corresponding overburden thickness and reflection coefficient (RC): (**a**) high yield GWpot VES stations, (**b**) medium yield GWpot VES stations, and (**c**) low yield GWpot VES stations. 42

**Figure 16.** Groundwater potential map of the vertical electrical sounding (VES) surveyed area.

#### *5.6. Geophysical Well Logs Interpretation*

Information obtained from technical reports of SPLs and NRLs (short and long) and drilling protocols show that the slightly denser thick and deep sandstone is an effective aquifer type for groundwater exploitation in the study area (Figure 17). The geophysical well logs approach has great significance in determining the exact location (depth) of any permeable aquifers and impermeable aquitards (Table 7). In this study, NRLs (short and long) were appropriately calibrated and quantitatively interpreted. Moreover, log measurements were converted to the apparent resistivity and adjusted for mud resistivity, bed thickness, borehole diameter, mud cake, and invasion to arrive at true resistivity (Figure 17). SPL interpretation can be complex, particularly in freshwater aquifers. This complexity commences to the perversion of groundwater and misinterpretations of spontaneous potential (SP) logging. SPLs record the potential or voltage caused by contact between a shale/clay layer and an aquifer. The natural flow of current and the SP curve were offered under the salinity conditions. The NRLs (short/long), SPLs, and drilling protocol at a depth of 152.4 m showed that the major lithology's units are clay, gravel-boulders, and sandstone (Table 8). The quality of groundwater measured by TDS is fresh. The static water level depth is about 88.3 m (Figure 17). The proposed slot opening, and the estimated discharge volume, are 1/40–1/50 and 11.35–13.24 cubic meters per hour (m3/h), respectively (Table 8).

**Figure 17.** Spontaneous potential (SP), short normal resistivity (SNR), and long normal resistivity (LNR) log curves obtained in Well -1 of the experimental site of the Marwatan Banda, Karak.



**Table 8.** Derived borehole lithology-based Normal resistivity logs (NRLs) (short/long) and spontaneous potential logs (SPLs).


#### **6. Discussion**

The Karak watershed, located in Northern Pakistan, has experienced significant economic development associated with hydrology and groundwater exploitation. The superficial resource depletion, the irregular spatial-temporal distribution of precipitation, and the deformation of the Indian and Eurasian tectonic plate environment, which affect the occurrence and movement of groundwater, together with widespread salt in the northern mountainous catchments, which is dissolved by runoff water and polluted groundwater due to deep infiltration, have made groundwater a key resource in the study area. However, the collaboration of remote sensing observations, aquifer geoelectrical properties and accurate hydrogeological measurements, and the optimization of groundwater influential factors are major challenges. Therefore, the GWpot mapping are essential for planning artificial recharge programs to mitigate groundwater decline [6]. The multicriteria decision-making (MCDM)-based multi-influence factor (MIF) model approach can be useful for groundwater resource management (GRM) and monitoring purposes, which is an efficient bivariate statistical technique mainly used to calculate the degree to which each dependent or independent conditioning factor influences the GWpot. The MIF model has become a powerful tool for delineating regional GWpot and narrowing down the target areas for conducting detailed hydrogeological and hydro-geophysical surveys in the scattered areas. However, in the MIF method, weights and ranks are subjectively assigned according to expert knowledge and literatures. In a comprehensive analysis, it is important to determine the weight of each category because the output result depends on the correct weight distribution. It is used to depict groundwater prediction zones taking into account various surface and subsurface hydrological influential factors. However, several studies report that the importance and predictive power of GCFs employed in GWpot assessment is usually controlled by geological, morphological, hydrological, and climatic environments [8–15,17]. According to Nampak [75], topographical features (e.g., elevation and slope) have a negative impact on GWpot, while lineament density and drainage density have positive impacts. Similar research reports that topographical, soil cover, structural and hydrogeological characteristics affect precipitation runoff and permeability, thereby affecting the occurrence of GWpot. Hou et al. [76] reported that lithology, altitude, and drainage density have a greater impact on the occurrence of GWpot, while land use and soil type have the least impacts. In this study, a GWpot map was generated based on the MIF model to identify regional GWpot of the Karak watershed. Several GCFs were concluded to have significant impacts on groundwater production. For example, the high GWpot zones on the final map are closely correlated to lineament density and drainage density. Usually, the lineaments indicate the areas of faults and fractures, leading to increased secondary porosity and permeability. This factor is of great significance in hydrogeology because it provides a pathway for groundwater infiltration. However, the lineament density is only an indirect indicator of the GWpot in the Karak watershed, because the lineaments usually show a permeable area. In the study area, a larger slope produces a smaller recharge, because surface water will quickly flow over the steep slope during rainfall, so there is not enough time for water seeping into the ground and recharge the unsaturated zone. However, the distribution of LU/LC usually depends on the subsurface soil and geological conditions, thereby increasing the groundwater recharge on the surfaces covered by vegetation (such as agricultural plants and forests).

The hydrogeological interpretation of the 2D high-resolution resistivity tomography dataset of six traverses revealed the prospect of groundwater at different depths with variation in the resistivities in the aquifer zone. The high resistivity of the subsurface geological sediments was well delineated, which shows a large resistivity contrast within the complex geological background in the study area. This phenomenon is suggested to be caused by different degrees of weathering, fracturing and saturated weathered/fractured part of the sediments in the Karak region. In future, four to five boreholes/wells will be drilled in potential areas identified by ERT and VES to check the availability of groundwater and the performance of geoelectric surveys. The analyzed regional GWpot, hydrogeological

and aquifer geoelectrical information provides a beneficial prospect for the development of GRM in the study area. However, the geoelectrical exploration methods can only locally verify the result of GWpot mapping, and they are too costly and time-consuming to cover the whole study area. The acquired results are expected to help practitioners to drill boreholes/wells in order to supply domestic water and irrigation in the Karak watershed of Northern Pakistan. Moreover, combined geospatial and geoelectrical methods through the MIFs model and Olayinka's basic criteria will help to assess groundwater resources in other similar areas worldwide.

#### **7. Conclusions**

This study addresses the applicability of the comprehensive MCDM-MIF model and hydro-geophysical investigation in GRM in the Karak watershed. The GIS-based MIF model facilitates the regional GWpot assessment using the topographical, geological, hydrological, and land-cover GCFs, meanwhile, the geophysical exploration and data interpretation reveals the hydrogeological structure and aquifer geoelectrical characteristics. The main findings are as follows:


**Author Contributions:** Conceptualization, U.K. and B.Z.; methodology, U.K. and H.F.; software, Z.J. and M.W.; validation, M.Y. and B.Z.; data curation, H.F. and M.Y.; writing—original draft preparation, U.K. and B.Z.; writing—review and editing, U.K. and B.Z.; funding acquisition, B.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (grant number 42072326 and 41772348) and the National Key Research and Development Program of China (grant number 2019YFC1805905 and 2017YFC0601503).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** No conflict of interest exists in the submission of this manuscript, and the manuscript was approved by all authors for publication.

#### **Appendix A**

**Table A1.** Summary results of VES data interpretation demonstrate the inferred lithologies corresponding to resistivity variation and hydrogeologic layers.



**Table A1.** *Cont.*


**Table A1.** *Cont.*

#### **References**

