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

: The optimization of groundwater conditioning factors (GCFs), the evaluation of groundwater potential (GW pot ), the hydrogeological characterization of aquifer geoelectrical properties and borehole lithological information are of great signiﬁcance in the complex decision-making processes of groundwater resource management (GRM). In this study, the regional GW pot of the Karak watershed in Northern Pakistan was ﬁrst evaluated by means of the multi-inﬂuence factors (MIFs) model of optimized GCFs through geoprocessing tools in geographical information system (GIS). The distribution of petrophysical properties indicated by the measured resistivity ﬂuctuations was then generated to locally verify the GW pot , and to analyze the hydrogeological and geoelectrical characteristics of aquifers. According to the weighted overlay analysis of MIFs, GW pot map was zoned into low, medium, high and very high areas, covering 9.7% (72.3 km 2 ), 52.4% (1307.7 km 2 ), 31.3% (913.4 km 2 ), and 6.6% (44.8 km 2 ) of the study area. The GW pot accuracy sequentially depends on the classiﬁcation criteria, the mean rating score, and the weights assigned to GCFs. The most inﬂuential factors are geology, lineament density, and land use/land cover followed by drainage density, slope, soil type, rainfall, elevation, and groundwater level ﬂuctuations. 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 ﬁve regional hydrological layers (i.e., coarse gravel and sand, silty sand mixed lithology, clayey sand/ﬁne sand, ﬁne sand/gravel, and clayey basement). The preliminary interpretation of ERT results highlights the complexity of the hydrogeological strata and reveals that GW pot is structurally and proximately constrained in the clayey sand and silicate aquifers (sandstone), which is of signiﬁcance 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 reﬂection coefﬁcient, 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 GW pot 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.


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][2][3][4][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][9][10][11][12].
The groundwater potential (GW pot ) 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 GW pot , 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][17][18][19][20][21][22][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 GW pot modeling approach, that can quickly, accurately, cost-effectively, and consequently monitor GW pot [23][24][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 GW pot [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][33][34]. With the proper GW pot 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 GW pot 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

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.

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.

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 mm 3 /year. The alluvial aquifer's average annual recharge is approximately equal to the average annual discharge, which is 2.7 mm 3 /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.

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 char-Figure 2. Regional geology map of the study area illustrates main structural and lithological units.

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 mm 3 /year. The alluvial aquifer's average annual recharge is approximately equal to the average annual discharge, which is 2.7 mm 3 /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.

GCF Analysis and Optimization
The evaluation of groundwater condition factors (GCFs) is essential to effectively determine an accurate groundwater potential (GW pot ) index [50]. GCFs should be considered in terms of regional topographical, geological, hydrological, and land use/land cover characteristics influencing the GW pot [15]. Therefore, the identification of the GW pot spatial distribution was performed by multi-criteria decision-making (MCDM) analysis Water 2021, 13, 1255 6 of 34 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). The drainage density (D d ) 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, D d 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 GW pot . The groundwater favorability is indirectly related to D d , which is related to surface runoff and permeability [52].
where DD represents drainage density, D l is the stream's length, and A is the watershed area (km 2 ). 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 L d 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 L d , the higher the favorability of GW pot .
where LD represent the lineament density, L i 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 inter-polation method. The flat and gentle areas, with less runoff, are more favorable for GW pot 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.
where P is the average precipitation depth, with p1, p2, p3 up to p n 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 GW pot , 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 GW pot . 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 GW pot 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 GW pot 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 GW pot of a given area.

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: 1.
Using RS and GIS toolkits, the database is ready to be input data for the MIF model, after the GCF' analysis and optimization described in Section 3.

2.
Once the GCFs are to be optimized, the weights and ranks of each GCFs are assigned for the multi-criteria decision-making (MCDM) MIF model, and the weighted/ranked GCFs are integrated through the weighted overlay analysis (WOA), based on the principle of superposition in a GIS environment to identify regional GW pot zones of the Karak watershed.

3.
The hydrogeological characteristics of the aquifers are evaluated by the interpretation of electrical resistivity tomography (ERT) and vertical electrical sounding (VES) data. Furthermore, the aquifer potential is further quantified through quantitative analysis of the resistivity mapping, overburden thickness mapping, and reflection coefficient mapping.

4.
To evaluate the accuracy of GW pot 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.
Water 2021, 13, x FOR PEER REVIEW 8 of 33 analysis 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.  The GW pot index is influenced by several hydrological, geological, topographical, environmental, and climatic variables [2]. By means of GCFs analysis and optimization, geology (G EO ), lineament density (L D ), drainage density (D D ), slope (S L ), soil type (S T ), rainfall (R F ), elevation (E L ), 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 GW pot 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 (CF h ) and minor (CF l ) 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 GW pot map is a weighted average. The estimated weight for each conditioning factor was obtained as a percentage using Equation (4).
where, CF h is the major weight of the condition factor, and CF l is the minor weight.
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).
where, CFh is the major weight of the condition factor, and CFl is the minor weight. Total Σ20.5 100   The GW pot 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 GW pot map, as summarized in Equation (5). A GW pot 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 GW pot into four levels of potentiality (i.e., low, medium, high, and very high).
where GW potz is the groundwater potential index, W i is the weight of each condition factor, R i 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 GW pot , 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 GW pot index at the corresponding position of GCFs.

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 GW pot 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 GW pot map with the observed GWT data. The overall accuracy was calculated using the following formula [64].
where, OA is the overall accuracy, C OW L 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].
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.

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). 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.

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).

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 hydro-stratigraphy 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].
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

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 GW pot 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].
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 GW pot and hydro-stratigraphy in the study area.

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 GW pot . 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).
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].
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.

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.

Evaluation of GCFs
The MIF model is an MCDM technique widely used for environmental management and has proven to effectively explain the GW pot 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 (D d ) 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/km 2 ), low (1.61-1.86 km/km 2 ), moderate (1.86-2.11 km/km 2 ), high (2.11-2.38 km/km 2 ), and very high (2.38-3.08 km/km 2 ) (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/km 2 category, indicating high infiltration and low runoff, and the lowest score was assigned to the 2.38-3.08 km/km 2 category (Table 4).
Lineament density (L d ) of the Karak watershed indirectly indicates the GW pot , 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 L d , the higher the probability of GW pot . Therefore, the highest rank was assigned to the 1.46-1.78 km/km 2 category and the lowest was assigned to the 0.17-0.45 km/km 2 category.
Rainfall (R F ) 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 R F , other precipitation characteristics, such as duration and intensity, are also important. For example, a long period of 20 mm R F has a more significant impact on groundwater recharge than a short period of 50 mm R F .
Geology (G EO ) characteristics govern the porosity and permeability of the hydrogeological layer, which in turn influences the formation and distribution of GW pot through physio-mechanical properties that control the water transmitting ability of the hydrogeological layer materials and the rate of groundwater flows. Therefore, the G EO 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).
Soil type (S T ) 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 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.

Assessment of GW pot
Using the weighted overlay analysis in the GIS environment, the GW pot 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 GW pot index, the GW pot 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 km 2 ), 52.4% (1307.7 km 2 ), 31.3% (913.4 km 2 ), and 6.6% (44.8 km 2 ) of the total area, respectively (Figure 7b,c). The spatial distribution of the various GW pot zones typically shows a mirror reflection of key factors. High and very high GW pot zones confirm their excellent capacities as sedimentary groundwater aquifers. The GW pot 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 GW pot . The northwestern, southeastern, and the central part limited regions typically have a low to medium GW pot , accounting for approximately 12.7% of the study area. Water 2021, 13, x FOR PEER REVIEW 16 of 33  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.

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

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 GW pot 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 GW pot ) 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 GW pot 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. 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.

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 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.

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. 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. 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 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 GW pot 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 GW pot 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 GW pot 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 GW pot 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.

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 GW pot 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 (h 1 , h 2 , h 3 , h 4 , and h 5 ), 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.

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

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

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). 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).

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.
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

GW pot 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 GW pot . 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. 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 isoresistivity 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].
(1) High GWpot: the overburden thickness is greater than 13 m with an RC less than 0.8.
(2) Medium GWpot: the overburden thickness is 13-30 m with an RC greater than or equal to 0.8. (3) Low GWpot: the overburden thickness is less than 13 m with an RC greater than or equal to 0.8. 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 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 thick-ness 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 GW pot 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].
(1) High GW pot : the overburden thickness is greater than 13 m with an RC less than 0.8.
(2) Medium GW pot : the overburden thickness is 13-30 m with an RC greater than or equal to 0.8. (3) Low GW pot : the overburden thickness is less than 13 m with an RC greater than or equal to 0.8.
Considering these criteria, the RC and overburden thickness were used to produce the parameters for categorizing VES stations by the GW pot , i.e., VES1, VES5, VES6, VES8, VES9, VES14, VES16, VES17, VES21, VES24, and VES26 have high yield GW pot (Figure 15a), VES3, VES7, VES11, VES13, VES19, VES20, and VES25 have medium yield GW pot (Figure 15b), and VES2, VES4, VES10, VES12, VES15, VES18, and VES22 have low yield GW pot (Figure 15c). Based on these groundwater potentiality variations among the VES stations, a final GW pot contour map of the surveyed area was generated, and it demonstrates that the northern, northeastern and eastern parts have excellent GW pot for future exploitation and development, while the low and medium GW pot 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 GW pot . 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.

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 (m 3 /h), respectively (Table 8).

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 (m 3 /h), respectively (Table 8).

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

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 GW pot 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 GW pot . The MIF model has become a powerful tool for delineating regional GW pot 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 GW pot assessment is usually controlled by geological, morphological, hydrological, and climatic environments [8][9][10][11][12][13][14][15]17]. According to Nampak [75], topographical features (e.g., elevation and slope) have a negative impact on GW pot , 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 GW pot . Hou et al. [76] reported that lithology, altitude, and drainage density have a greater impact on the occurrence of GW pot , while land use and soil type have the least impacts. In this study, a GW pot map was generated based on the MIF model to identify regional GW pot of the Karak watershed. Several GCFs were concluded to have significant impacts on groundwater production. For example, the high GW pot 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 GW pot 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 GW pot , 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 GW pot 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.

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 GW pot 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: (1) According to MCDM-MIF model, approximately 9.7% (72.3 km 2 ), 52.4% (1307.7 km 2 ), 31.3% (913.4 km 2 ), and 6.6% (44.8 km 2 ) areas of the total Karak watershed are classified into the low, medium, high, and very high GW pot , respectively. The southern, southeastern, and the limited northeastern areas have high to medium GW pot due to the distribution of Quaternary alluvial and agricultural land with high infiltration capacity. The final GW pot map will help to manage sustainable groundwater resources in the study area. (2) The predictive performance of MCDM-MIF model is consistent with the groundwater level (GWL) data (as AUC value is 68%, confusion matrix is 68%, and Kappa (K) analysis is 65%). (3) The ERT approach with an optimal compromise between electrode distance and profile length highlights the complexity of hydrogeological layers and reveal that GW pot is structurally controlled and adjacent to clayey sand and silicate aquifers (sandstone). The identified drilling locations on ERT traverses are of great significance for the expansion of drinking water supply and irrigation in the future. The performance analysis between ERT-predicted lithology and well-log lithology indicates suitable matches. (4) Hydro-stratigraphic information followed by apparent resistivity distribution at each VES station shows that the study area is mainly composed of coarse gravel and sand, followed by clayey sand with saline water. According to Olayinka's basic standards, the aquifer geoelectrical characteristics, e.g., reflection coefficient, aquifer overburden thickness and apparent resistivity distribution, were conceptualized. The interpreted potential zones based on VES show satisfactory matches with MIF-based groundwater potential. The drilling protocol and well logs data interpretation of NRLs (short/long) and SPLs reveal that deep sandstone is an effective aquifer type for the groundwater exploitation in the study area.  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.