Next Article in Journal
Optimization and Degradation Studies on Hexahydro-1,3,5-Trinitro-1,3,5-Triazine (RDX) with Selected Indigenous Microbes under Aerobic Conditions
Next Article in Special Issue
Analysis of the Spatiotemporal Annual Rainfall Variability in the Wadi Cheliff Basin (Algeria) over the Period 1970 to 2018
Previous Article in Journal
A Moment-Based Chezy Formula for Bed Shear Stress in Varied Flow
Previous Article in Special Issue
Performance Evaluation of a Two-Parameters Monthly Rainfall-Runoff Model in the Southern Basin of Thailand
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
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

1
MoE Key Laboratory of Metallogenic Prediction of Nonferrous Metals & Geological Environment Monitoring/School of Geosciences & Info-Physics, Central South University, Changsha 410083, China
2
Department of Geology, Bacha Khan University, Charsadda 24600, Pakistan
*
Author to whom correspondence should be addressed.
Water 2021, 13(9), 1255; https://doi.org/10.3390/w13091255
Submission received: 5 April 2021 / Revised: 27 April 2021 / Accepted: 29 April 2021 / Published: 30 April 2021
(This article belongs to the Special Issue Hydrology in Water Resources Management)

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.

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,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 (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), weights-of-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 GWpot modeling approach, that can quickly, accurately, cost-effectively, and consequently monitor GWpot [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 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,33,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.

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.

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).
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].
DD = l = o n D l ( km ) A ( km ) 2 ( km 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.
LD = i = 0 n L i   ( km ) A   ( sq .   km ) ( km 1 )
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 = l = 1 n p l N
where P is the average precipitation depth, with p1, p2, p3 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:
  • 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.
  • 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 GWpot zones of the Karak watershed.
  • 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.
  • 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.

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).
Score = [ ( C F h + C F l )   ( C F h + C F l ) × 100 ]
where, CFh is the major weight of the condition factor, and CFl is the minor weight.

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).
  GW potz = l = 0 n W i × R i = DD c DD w +   LD c LD w + RF c RF w +   SL c SL w +   EL c EL w + GEO c GEO w +   LULC c LULC w +   ST c ST w + GLF c GLF w
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].
OA =   C O W L   O W L
where, OA is the overall accuracy, C O W 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].
K = [   CV   % CAOV %   TC CAOV % ]
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).

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 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].
ρ = K V I
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).
ρ a = π { ( AB 2 ) 2 ( MN 2 ) 2 MN } Ra
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 computer-assisted 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].
RC = { ( ρ n ρ ( n 1 ) ( ρ n + ρ ( n 1 ) }
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/km2 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/km2 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).
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.

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

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

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

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.
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].
(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 in predicting GWpot.

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

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 multi-criteria 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,9,10,11,12,13,14,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:
(1)
According to MCDM-MIF model, approximately 9.7% (72.3 km2), 52.4% (1307.7 km2), 31.3% (913.4 km2), and 6.6% (44.8 km2) areas of the total Karak watershed are classified into the low, medium, high, and very high GWpot, respectively. The southern, southeastern, and the limited northeastern areas have high to medium GWpot due to the distribution of Quaternary alluvial and agricultural land with high infiltration capacity. The final GWpot 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 GWpot 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.

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. Summary results of VES data interpretation demonstrate the inferred lithologies corresponding to resistivity variation and hydrogeologic layers.
VES No.RMS Error (%)No. of LayersResistivity
(Ohm.m)
Thickness
(m)
Depth
(m)
Reflection
Coefficient
Inferred Hydro-Stratigraphic
Lithology
12.9Layer 1923.12.42.40.6723Coarse gravel and sand
Layer 25787.29.8Silty sand mixed lithology
Layer 3685.428.438.2Clayey sand
Layer 45412.650.8Fine sand and gravels
23.2Layer 17481.91.90.9284Coarse gravel and sand
Layer 2598.88.610.5Silty sand and gravels
Layer 357.816.426.9Fine sands and gravels
Layer 434.2InfiniteInfiniteClayey sand and saline water
33.8Layer 1899.32.42.40.8671Coarse gravel and sand
Layer 26738.210.6Silty sand mixed lithology
Layer 3487.418.729.3Silty sand and gravels
Layer 448.314.443.7Fine sand and gravel
Layer 5268.652.3Clayey sand and saline water
44.5Layer 19541.31.30.9752Topsoil
Layer 2637.53.64.9Silty sand and gravels
Layer 3437.412.3Fine sand and gravels
Layer 4683.1InfiniteInfiniteClayey sands
52.8Layer 1854.74.54.50.7841Coarse gravel and sand
Layer 2532.112.817.3Silty sand and gravels
Layer 3678.226.243.5Clayey sands
Layer 443.26.650.1Clayey sand and saline water
61.7Layer 17013.23.20.8911Coarse gravel and sand
Layer 2693.44.67.8Silty sand mixed lithology
Layer 3964.55.112.9Clayey sands
Layer 445.8InfiniteInfiniteClayey sand and saline water
72.6Layer 11093.53.83.80.8497Coarse gravel and sand
Layer 2601.512.416.2Silty sand mixed lithology
Layer 3106518.434.6Clayey sand
Layer 4738.242.8Fine sand and gravels
84.0Layer 11108.54.54.60.7313Topsoil
Layer 2598.114.118.7Silty sand and gravels
Layer 3376.722.441.1Fine sand and gravels
Layer 4985.716.157.2Clayey sands layer
94.7Layer 17352.22.20.6861Alluvium
Layer 23235.67.8Silty sand fine lithology
Layer 3675.812.420.2Silty sand and gravels
Layer 424.71636.2Fine sands and gravels
Layer 5746.442.6Clayey sand and saline water
102.5Layer 1967.81.81.80.9248Coarse gravel and sand
Layer 25604.26Silty sand mixed lithology
Layer 3856.66.812.8Clayey sand
Layer 442.1InfiniteInfiniteFine sand and gravels
113.7Layer 17875.15.10.8410Coarse gravel and sand
Layer 2445.812.217.3Silty sand and gravels
Layer 3943.218.135.4Clayey sands
Layer 4446.842.2Clayey sand and saline water
122.8Layer 1799.81.61.60.9523Coarse gravel and sand
Layer 25884.15.7Silty sand mixed lithology
Layer 31098.34.610.3Clayey sand
Layer 4572.412.7Fine sand and gravels
132.0Layer 1800.31.31.30.8557Coarse gravel and sand
Layer 24548.29.5Silty sand and gravels
Layer 3822.716.225.7Clayey sands
Layer 438.610.636.3Clayey sand and saline water
144.0Layer 18662.62.60.6719Alluvium
Layer 255410.112.7Silty sand fine lithology
Layer 3600.416.829.5Silty sand and gravels
Layer 489.322.842.3Fine sands and gravels
Layer 534.98.250.5Clayey sand and saline water
152.8Layer 1766.12.42.40.9643Coarse gravel and sand
Layer 25436.18.5Silty sand and gravels
Layer 39854.112.6Clayey sands
Layer 427InfiniteInfiniteClayey sand and saline water
163.0Layer 1812.94.24.20.6537Coarse gravel and sand
Layer 2553.48.212.4Silty sand and gravels
Layer 3985.620.833.2Clayey sands
Layer 45824.257.4Clayey sand and saline water
174.5Layer 111093.43.40.5855Topsoil
Layer 26438.712.1Silty sand and gravels
Layer 37615.928Fine sand and gravels
Layer 48322250Clayey sands
182.8Layer 17411.81.80.9778Coarse gravel and sand
Layer 25335.16.9Silty sand and gravels
Layer 39326.012.9Clayey sands
Layer 465.4InfiniteInfiniteClayey sand and saline water
191.6Layer 1979.85.25.20.8923Coarse gravel and sand
Layer 25688.413.6Silty sand and gravels
Layer 3732.616.830.4Clayey sands
Layer 46410.841.2Clayey sand and saline water
204.2Layer 19053.23.20.9211Coarse gravel and sand
Layer 25486.19.3Silty sand mixed lithology
Layer 378819.528.8Clayey sand
Layer 4341543.8Fine sand and gravels
212.4Layer 17454.24.20.7626Alluvium
Layer 262310.414.6Silty sand fine lithology
Layer 3522.318.633.2Silty sand and gravels
Layer 4376.239.4Fine sands and gravels
Layer 584.28.247.6Clayey sand and saline water
223.7Layer 1865.61.61.60.9429Coarse gravel and sand
Layer 25904.56.1Silty sand and gravels
Layer 3955.76.812.9Clayey sands
Layer 467InfiniteInfiniteClayey sand and saline water
233.3Layer 1906.45.65.60.6930Alluvium
Layer 257810.215.8Silty sand fine lithology
Layer 335618.834.6Silty sand and gravels
Layer 46510.445Fine sands and gravels
Layer 545.68.453.4Clayey sand and saline water
242.8Layer 18424.44.40.7889Coarse gravel and sand
Layer 2600.412.416.8Silty sand and gravels
Layer 378922.839.6Clayey sands
Layer 437.214.754.3Clayey sand and saline water
252.5Layer 18902.42.40.8663Coarse gravel and sand
Layer 2566.36.89.2Silty sand and gravels
Layer 374218.327.5Clayey sands
Layer 49816.243.7Clayey sand and saline water
263.0Layer 1975.81.71.70.9184Coarse gravel and sand
Layer 256310.211.9Silty sand and gravels
Layer 3732.11829.9Clayey sands
Layer 463.96.636.5Clayey sand and saline water

References

  1. Wang, L.F.; Chen, Q.; Long, X.P.; Wu, X.B.; Sun, L. Comparative analysis of groundwater fluorine levels and other characteristics in two areas of Laizhou Bay and its explanation on fluorine enrichment. Water Sci. Technol. Water Supply 2015, 15, 384–394. [Google Scholar] [CrossRef]
  2. Pinto, D.; Shrestha, S.; Babel, M.S.; Ninsawat, S. Delineation of groundwater potential zones in the Comoro watershed, Timor Leste using GIS, remote sensing and analytic hierarchy process (AHP) technique. Appl. Water Sci. 2017, 7, 503–519. [Google Scholar] [CrossRef] [Green Version]
  3. Berhanu, K.G.; Hatiye, S.D. Identification of Groundwater Potential Zones Using Proxy Data: Case study of Megech Watershed, Ethiopia. J. Hydrol. Reg. Stud. 2020, 28, 100676. [Google Scholar] [CrossRef]
  4. Chen, Q.; Jia, C.; Wei, J.; Dong, F.; Yang, W.; Hao, D.; Jia, Z.; Ji, Y. Geochemical process of groundwater fluoride evolution along global coastal plains: Evidence from the comparison in seawater intrusion area and soil salinization area. Chem. Geol. 2020, 552, 119779. [Google Scholar] [CrossRef]
  5. Hao, Q.; Xiao, Y.; Chen, K.; Zhu, Y.; Li, J. Comprehensive Understanding of Groundwater Geochemistry and Suitability for Sustainable Drinking Purposes in Confined Aquifers of the Wuyi Region, Central North China Plain. Water 2020, 12, 3052. [Google Scholar] [CrossRef]
  6. Goldman, M.; Kafri, U. Hydrogeophysical applications in coastal aquifers. In Applied Hydrogeophysics; Springer: Dordrecht, The Netherlands, 2006; pp. 233–254. [Google Scholar]
  7. Hasan, M.; Shang, Y.; Jin, W.; Shao, P.; Yi, X.; Akhter, G. Geophysical Assessment of Seawater Intrusion into Coastal Aquifers of Bela Plain, Pakistan. Water 2020, 12, 3408. [Google Scholar] [CrossRef]
  8. Baalousha, H. Using CRD method for quantification of groundwater recharge in the Gaza Strip, Palestine. Environ. Geol. 2005, 48, 889–900. [Google Scholar] [CrossRef]
  9. Shaban, A.; Khawlie, M.; Abdallah, C. Use of remote sensing and GIS to determine recharge potential zones: The case of Occidental Lebanon. Hydrogeol. J. 2005, 14, 433–443. [Google Scholar] [CrossRef]
  10. Baalousha, H. Stochastic water balance model for rainfall recharge quantification in Ruataniwha Basin, New Zealand. Environ. Geol. 2009, 58, 85. [Google Scholar] [CrossRef]
  11. Zghibi, A.; Merzougui, A.; Chenini, I.; Ergaieg, K.; Zouhri, L.; Tarhouni, J. Groundwater vulnerability analysis of Tunisian coastal aquifer: An application of DRASTIC index method in GIS environment. Groundw. Sustain. Dev. 2016, 2-3, 169–181. [Google Scholar] [CrossRef]
  12. Zabihi, M.; Pourghasemi, H.R.; Pourtaghi, Z.S.; Behzadfar, M. GIS-based multivariate adaptive regression spline and random forest models for groundwater potential mapping in Iran. Environ. Earth Sci. 2016, 75, 1–19. [Google Scholar] [CrossRef]
  13. Arshad, A.; Zhang, Z.; Zhang, W.; Dilawar, A. Mapping favorable groundwater potential recharge zones using a GIS-based analytical hierarchical process and probability frequency ratio model: A case study from an agro-urban region of Pakistan. Geosci. Front. 2020, 11, 1805–1819. [Google Scholar] [CrossRef]
  14. Rahmati, O.; Samani, A.N.; Mahdavi, M.; Pourghasemi, H.R.; Zeinivand, H. Groundwater potential mapping at Kurdistan region of Iran using analytic hierarchy process and GIS. Arab. J. Geosci. 2015, 8, 7059–7071. [Google Scholar] [CrossRef]
  15. Cai, Z.; Ofterdinger, U. Analysis of groundwater-level response to rainfall and estimation of annual recharge in fractured hard rock aquifers, NW Ireland. J. Hydrol. 2016, 535, 71–84. [Google Scholar] [CrossRef] [Green Version]
  16. Pourghasemi, H.R.; Pradhan, B.; Gokceoglu, C.; Mohammadi, M.; Moradi, H.R. Application of weights-of-evidence and certainty factor models and their comparison in landslide susceptibility mapping at Haraz watershed, Iran. Arab. J. Geosci. 2013, 6, 2351–2365. [Google Scholar] [CrossRef]
  17. Baalousha, H.M. Using Monte Carlo simulation to estimate natural groundwater recharge in Qatar. Model. Earth Syst. Environ. 2016, 2, 1–7. [Google Scholar] [CrossRef] [Green Version]
  18. Pradhan, B.; Seeni, M.I.; Kalantar, B. Performance evaluation and sensitivity analysis of expert-based, statistical, machine learning, and hybrid models for producing landslide susceptibility maps. In Laser Scanning Applications in Landslide Assessment; Springer: Cham, Germany, 2017; pp. 193–232. [Google Scholar]
  19. Xiao, Y.; Gu, X.; Yin, S.; Pan, X.; Shao, J.; Cui, Y. Investigation of Geochemical Characteristics and Controlling Processes of Groundwater in a Typical Long-Term Reclaimed Water Use Area. Water 2017, 9, 800. [Google Scholar] [CrossRef] [Green Version]
  20. Abrams, W.; Ghoneim, E.; Shew, R.; LaMaskin, T.; Al-Bloushi, K.; Hussein, S.; AbuBakr, M.; Al-Mulla, E.; Al-Awar, M.; El-Baz, F. Delineation of groundwater potential (GWP) in the northern United Arab Emirates and Oman using geospatial technologies in conjunction with Simple Additive Weight (SAW), Analytical Hierarchy Process (AHP), and Probabilistic Frequency Ratio (PFR) techniques. J. Arid. Environ. 2018, 157, 77–96. [Google Scholar] [CrossRef]
  21. Xiao, Y.; Shao, J.; Frape, S.K.; Cui, Y.; Dang, X.; Wang, S.; Ji, Y. Groundwater origin, flow regime and geochemical evolution in arid endorheic watersheds: A case study from the Qaidam Basin, northwestern China. Hydrol. Earth Syst. Sci. 2018, 22, 4381–4400. [Google Scholar] [CrossRef] [Green Version]
  22. Díaz-Alcaide, S.; Martínez-Santos, P. Review: Advances in groundwater potential mapping. Hydrogeol. J. 2019, 27, 2307–2324. [Google Scholar] [CrossRef]
  23. Raju, R.S.; Raju, G.S.; Rajasekhar, M. Identification of groundwater potential zones in Mandavi River basin, Andhra Pradesh, India using remote sensing, GIS and MIF techniques. HydroResearch 2019, 2, 1–11. [Google Scholar] [CrossRef]
  24. Das, S.; Gupta, A.; Ghosh, S. Exploring groundwater potential zones using MIF technique in semi-arid region: A case study of Hingoli district, Maharashtra. Spat. Inf. Res. 2017, 25, 749–756. [Google Scholar] [CrossRef]
  25. Zygouri, V.; Verroios, S.; Anninou, A.; Koukouvelas, I. Fuzzy cognitive mapping vs. statistic approach of geomorphological analysis of active faults: Preliminary results from a case study within the Corinth Graben, Greece. IFAC PapersOnLine 2018, 51, 372–377. [Google Scholar] [CrossRef]
  26. Baalousha, H.M. Development of a groundwater flow model for the highly parameterized Qatar aquifers. Model. Earth Syst. Environ. 2016, 2, 67. [Google Scholar] [CrossRef] [Green Version]
  27. Srivastava, P.K.; Bhattacharya, A.K. Groundwater assessment through an integrated approach using remote sensing, GIS and resistivity techniques: A case study from a hard rock terrain. Int. J. Remote Sens. 2006, 27, 4599–4620. [Google Scholar] [CrossRef]
  28. Farid, A.; Jadoon, K.; Akhter, G.; Iqbal, M.A. Hydrostratigraphy and hydrogeology of the western part of Maira area, Khyber Pakhtunkhwa, Pakistan: A case study by using electrical resistivity. Environ. Monit. Assess. 2013, 185, 2407–2422. [Google Scholar] [CrossRef]
  29. Farid, A.; Khalid, P.; Jadoon, K.Z.; Iqbal, M.A.; Small, J. An application of variogram modelling for electrical resistivity soundings to characterize depositional system and hydrogeology of Bannu Basin, Pakistan. Geosci. J. 2017, 21, 819–839. [Google Scholar] [CrossRef]
  30. Mondal, N.C.; Rao, A.V.; Singh, V.P. Efficacy of electrical resistivity and induced polarization methods for revealing fluoride contaminated groundwater in granite terrain. Environ. Monit. Assess. 2010, 168, 103–114. [Google Scholar] [CrossRef]
  31. Son, Y. Assessment of concentration in contaminated soil by potentially toxic elements using electrical properties. Environ. Monit. Assess. 2010, 176, 1–11. [Google Scholar] [CrossRef]
  32. Devi, S.P.; Srinivasulu, S.; Raju, K.K. Delineation of groundwater potential zones and electrical resistivity studies for groundwater exploration. Environ. Earth Sci. 2001, 40, 1252–1264. [Google Scholar] [CrossRef]
  33. Anechana, R.; Noye, R.M.; Menyeh, A.; Manu, E.; Okrah, C. Electromagnetic Method and Vertical Electrical Sounding for Groundwater Potential Assessment of Kintampo North Municipality of Ghana. J. Environ. Earth Sci. 2015, 5, 9–16. [Google Scholar]
  34. Gupta, G.; Patil, S.N.; Padmane, S.T.; Erram, V.C.; Mahajan, S. Geoelectric investigation to delineate groundwater potential and recharge zones in Suki river basin, north Maharashtra. J. Earth Syst. Sci. 2015, 124, 1487–1501. [Google Scholar] [CrossRef] [Green Version]
  35. Olorunfemi, M.; Fasuyi, S. Aquifer types and the geoelectric/hydrogeologic characteristics of part of the central basement terrain of Nigeria (Niger State). J. Afr. Earth Sci. 1993, 16, 309–317. [Google Scholar] [CrossRef]
  36. Frid, V.; Sharabi, I.; Frid, M.; Averbakh, A. Leachate detection via statistical analysis of electrical resistivity and induced polarization data at a waste disposal site (Northern Israel). Environ. Earth Sci. 2017, 76, 233. [Google Scholar] [CrossRef]
  37. Kumar, D.; Mondal, S.; Nandan, M.J.; Harini, P.; Sekhar, B.M.V.S.; Sen, M.K. Two-dimensional electrical resistivity tomography (ERT) and time-domain-induced polarization (TDIP) study in hard rock for groundwater investigation: A case study at Choutuppal Telangana, India. Arab. J. Geosci. 2016, 9, 355. [Google Scholar] [CrossRef]
  38. Hamzah, U.; Samsudin, A.R.; Malim, E.P. Groundwater investigation in Kuala Selangor using vertical electrical sounding (VES) surveys. Environ. Earth Sci. 2006, 51, 1349–1359. [Google Scholar] [CrossRef]
  39. Olasehinde, P.I.; Bayewu, O.O. Evaluation of electrical resistivity anisotropy in geological mapping: A case study of Odo Ara, West Central Nigeria. Afr. J. Environ. Sci. Technol. 2011, 5, 553–566. [Google Scholar]
  40. Webb, S.J.; Ngobeni, D.; Jones, M.; Abiye, T.; Devkurran, N.; Goba, R.; Lee, M.; Burrows, D.; Pellerin, L. Hydrogeophysical investigation for groundwater at the Dayspring Children’s Village, South Africa. Lead. Edge 2011, 30, 434–440. [Google Scholar]
  41. Manu, E.; Agyekum, W.A.; Duah, A.A.; Mainoo, P.A.; Okrah, C.; Asare, V.S. Improving Access to Potable Water Supply using Integrated Geophysical Approach in a Rural Setting of Eastern Ghana. Elixir Environ. For. 2016, 95, 40714–40719. [Google Scholar]
  42. Kruseman, G.P.; Naqavi, S.A.H. Hydrogeology and groundwater resources of the North-West Frontier Province Pakistan. WAPDA Hydrogeol. Dir. Peshawar 1988. [Google Scholar]
  43. Gee, E.; Gee, D. Overview of the geology and structure of the Salt Range, with observations on related areas of northern Pakistan. Geol. Soc. Am. Spec. Paper 1989, 232, 95–112. [Google Scholar]
  44. Ali, F.; Khan, M.I.; Ahmad, S.; Rehman, G.; Rehman, I.; Ali, T.H. Range front structural style: An example from Surghar Range, North Pakistan. J. Himal. Earth Sci. 2014, 47, 193. [Google Scholar]
  45. Nazir-ur-Rehman, A.S.; Fayaz, A.; Iftikhar, A.; Amin, S. Joints/fractures analyses of Shinawah area, district Karak, Khyber Pakhtunkhwa, Pakistan. J. Himal. Earth Sci. Vol. 2017, 50, 93–113. [Google Scholar]
  46. Jaswal, T.M.; Lillie, R.J.; Lawrence, R.D. Structure and evolution of the northern Potwar deformed zone, Pakistan. AAPG Bull. 1997, 81, 308–328. [Google Scholar]
  47. Akhter, G.; Farid, A.; Ahmad, Z. Determining the depositional pattern by resistivity–seismic inversion for the aquifer system of Maira area, Pakistan. Environ. Monit. Assess. 2012, 184, 161–170. [Google Scholar] [CrossRef]
  48. Ullah, K.; Arif, M.; Shah, M.T. Petrography of sandstones from the Kamlial and Chinji formations, southwestern Kohat plateau, NW Pakistan: Implications for source lithology and paleoclimate. J. Himal. Earth Sci. 2006, 39, 1–13. [Google Scholar]
  49. Shahzad, F.; Mahmood, S.A.; Gloaguen, R. Drainage network and lineament analysis: An approach for Potwar Plateau (Northern Pakistan). J. Mt. Sci. 2009, 6, 14–24. [Google Scholar] [CrossRef]
  50. Kalantar, B.; Al-Najjar, H.A.H.; Pradhan, B.; Saeidi, V.; Halin, A.A.; Ueda, N.; Naghibi, S.A.; Najjar, A.-. Ueda Optimized Conditioning Factors Using Machine Learning Techniques for Groundwater Potential Mapping. Water 2019, 11, 1909. [Google Scholar] [CrossRef] [Green Version]
  51. Magesh, N.; Chandrasekar, N.; Soundranayagam, J.P. Delineation of groundwater potential zones in Theni district, Tamil Nadu, using remote sensing, GIS and MIF techniques. Geosci. Front. 2012, 3, 189–196. [Google Scholar] [CrossRef] [Green Version]
  52. Sener, E.; Davraz, A.; Ozcelik, M. An integration of GIS and remote sensing in groundwater investigations: A case study in Burdur, Turkey. Hydrogeol. J. 2005, 13, 826–834. [Google Scholar] [CrossRef]
  53. Edet, A.E.; Okereke, C.S.; Teme, S.C.; Esu, E.O. Application of remote-sensing data to groundwater exploration: A case study of the Cross River State, southeastern Nigeria. Hydrogeol. J. 1998, 6, 394–404. [Google Scholar] [CrossRef]
  54. Pavelic, P.; Giordano, M.; Keraita, B.N.; Ramesh, V.; Rao, T. Groundwater Availability and Use in Sub-Saharan Africa: A Review of 15 Countries; International Water Management Institute: Colombo, Sri Lanka, 2012. [Google Scholar]
  55. Ndlovu, S.; Mpofu, V.; Manatsa, D.; Muchuweni, E. Mapping groundwater aquifers using dowsing, slingram electromagnetic survey method and vertical electrical sounding jointly in the granite rock formation: A case of Matshetshe rural area in Zimbabwe. J. Sustain. Dev. Afr. 2010, 12, 199–208. [Google Scholar]
  56. Steelman, C.M.; Kennedy, C.S.; Capes, D.C.; Parker, B.L. Electrical resistivity dynamics beneath a fractured sedimentary bedrock riverbed in response to temperature and groundwater–surface water exchange. Hydrol. Earth Syst. Sci. 2017, 21, 3105–3123. [Google Scholar] [CrossRef] [Green Version]
  57. Stephenson, D.J. Rockfill in Hydraulic Engineering; Elsevier: Amsterdam, The Netherlands, 1979. [Google Scholar]
  58. Hutti, B.; Nijagunappa, R. Identification of groundwater potential zone using geoinformatics in Ghataprabha basin, North Karnataka, India. Int. J. Geomat. Geosci. 2011, 2, 91. [Google Scholar]
  59. Zghibi, A.; Mirchi, A.; Msaddek, M.H.; Merzougui, A.; Zouhri, L.; Taupin, J.-D.; Chekirbane, A.; Chenini, I.; Tarhouni, J. Using Analytical Hierarchy Process and Multi-Influencing Factors to Map Groundwater Recharge Zones in a Semi-arid Mediterranean Coastal Aquifer. Water 2020, 12, 2525. [Google Scholar] [CrossRef]
  60. Malczewski, J. On the Use of Weighted Linear Combination Method in GIS: Common and Best Practice Approaches. Trans. GIS 2000, 4, 5–22. [Google Scholar] [CrossRef]
  61. Chen, W.; Panahi, M.; Khosravi, K.; Pourghasemi, H.R.; Rezaie, F.; Parvinnezhad, D. Spatial prediction of groundwater potentiality using ANFIS ensembled with teaching-learning-based and biogeography-based optimization. J. Hydrol. 2019, 572, 435–448. [Google Scholar] [CrossRef]
  62. Rajaveni, S.P.; Brindha, K.; Elango, L. Geological and geomorphological controls on groundwater occurrence in a hard rock region. Appl. Water Sci. 2017, 7, 1377–1389. [Google Scholar] [CrossRef] [Green Version]
  63. Chepchumba, M.C.; Raude, J.M.; Sang, J.K. Geospatial delineation and mapping of groundwater potential in Embu County, Kenya. Acque Sotter. Ital. J. Groundw. 2019, 8. [Google Scholar] [CrossRef]
  64. Jensen, J.R. Introductory Digital Image Processing: A Remote Sensing Perspective; Prentice-Hall Inc.: Upper Sadder River, NJ, USA, 1996. [Google Scholar]
  65. Usman, M.; Liedl, R.; Shahid, M.A.; Abbas, A. Land use/land cover classification and its change detection using multi-temporal MODIS NDVI data. J. Geogr. Sci. 2015, 25, 1479–1506. [Google Scholar] [CrossRef]
  66. Muchingami, I.; Hlatywayo, D.; Nel, J.; Chuma, C. Electrical resistivity survey for groundwater investigations and shallow subsurface evaluation of the basaltic-greenstone formation of the urban Bulawayo aquifer. Phys. Chem. Earth, Parts A/B/C 2012, 50, 44–51. [Google Scholar] [CrossRef]
  67. Sasaki, Y. Resolution of resistivity tomography inferred from numerical simulation 1. Geophys. Prospect. 1992, 40, 453–463. [Google Scholar] [CrossRef]
  68. Loke, M.H.; Barker, R.D. Rapid least-squares inversion of apparent resistivity pseudosections by a quasi-Newton method1. Geophys. Prospect. 1996, 44, 131–152. [Google Scholar] [CrossRef]
  69. Ernstson, K.; Kirsch, R. Geoelectrical methods. In Groundwater Geophysics: A Tool for Hydrogeology; Spinger: Berlin, Germany, 2006; pp. 84–117. [Google Scholar]
  70. Lenkey, L.; Hámori, Z.; Mihálffy, P. Investigating the hydrogeology of a water-supply area using direct-current vertical electrical soundings. Geophysics 2005, 70, H11–H19. [Google Scholar] [CrossRef]
  71. Xiao, Y.; Yin, S.; Hao, Q.; Gu, X.; Pei, Q.; Zhang, Y. Hydrogeochemical appraisal of groundwater quality and health risk in a near-suburb area of North China. J. Water Supply Res. Technol. 2020, 69, 55–69. [Google Scholar] [CrossRef]
  72. Anomohanran, O. Investigation of Groundwater Potential in some selected towns in Delta North District of Nigeria. Int. J. Appl. Sci. Technol. 2013, 3, 61–66. [Google Scholar]
  73. Olayinka, A. Non uniqueness in the interpretation of bedrock resistivity from sounding curves and its hydrological implications. Water Resour. J. NAH 1996, 7, 55–60. [Google Scholar]
  74. Bayewu, O.O.; Oloruntola, M.O.; Mosuro, G.O.; Laniyan, T.A.; Ariyo, S.O.; Fatoba, J.O. Geophysical evaluation of groundwater potential in part of southwestern Basement Complex terrain of Nigeria. Appl. Water Sci. 2017, 7, 4615–4632. [Google Scholar] [CrossRef] [Green Version]
  75. Nampak, H.; Pradhan, B.; Manap, M.A. Application of GIS based data driven evidential belief function model to predict groundwater potential zonation. J. Hydrol. 2014, 513, 283–300. [Google Scholar] [CrossRef]
  76. Hou, E.; Wang, J.; Chen, W. A comparative study on groundwater spring potential analysis based on statistical index, index of entropy and certainty factors models. Geocarto Int. 2018, 33, 754–769. [Google Scholar] [CrossRef]
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.
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.
Water 13 01255 g001
Figure 2. Regional geology map of the study area illustrates main structural and lithological units.
Figure 2. Regional geology map of the study area illustrates main structural and lithological units.
Water 13 01255 g002
Figure 3. Framework to delineate groundwater potential and to identify hydrogeological characteristics.
Figure 3. Framework to delineate groundwater potential and to identify hydrogeological characteristics.
Water 13 01255 g003
Figure 4. Interrelationship between the GCFs concerning the GWpot index.
Figure 4. Interrelationship between the GCFs concerning the GWpot index.
Water 13 01255 g004
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.
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.
Water 13 01255 g005
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 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.
Water 13 01255 g006
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.
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.
Water 13 01255 g007
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.
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.
Water 13 01255 g008
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.
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.
Water 13 01255 g009
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 10. VES data interpretation result-based partial curve matching (PCM) along 26 VES stations: (a) hydrologic layers depth variation; (b) hydrologic layers resistivity variation.
Water 13 01255 g010
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.
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.
Water 13 01255 g011
Figure 12. Correlation of VES data interpretation results with borehole lithological information.
Figure 12. Correlation of VES data interpretation results with borehole lithological information.
Water 13 01255 g012
Figure 13. (a) Reflection coefficient (RC); (b) overburden thickness along 26 VES stations in the surveyed area.
Figure 13. (a) Reflection coefficient (RC); (b) overburden thickness along 26 VES stations in the surveyed area.
Water 13 01255 g013
Figure 14. (a) Reflection coefficient map; (b) overburden thickness map; (c) apparent resistivity map based on the interpretation of VES data.
Figure 14. (a) Reflection coefficient map; (b) overburden thickness map; (c) apparent resistivity map based on the interpretation of VES data.
Water 13 01255 g014
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.
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.
Water 13 01255 g015
Figure 16. Groundwater potential map of the vertical electrical sounding (VES) surveyed area.
Figure 16. Groundwater potential map of the vertical electrical sounding (VES) surveyed area.
Water 13 01255 g016
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.
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.
Water 13 01255 g017
Table 1. Lithological characteristics of the Karak watershed.
Table 1. Lithological characteristics of the Karak watershed.
ProductFormation NamesLithology Characteristics
(C Fm)Chinji formationSandstones and shales (abundant quartz with subordinate feldspars)
(DS)Darzinda shaleDark-brown to gray claystone and subordinate fossiliferous marl beds
(K Fm)Kamlial formationMainly composed of sandstone (subordinate feldspars, lithic grains, micas)
(J/T)Jurassic or Triassic rocksSandstone, siltstone, shale and dolomite
DP FmDhok Pathan formationEqual amount of sandstone and clay
(K Fm)Kohat formationMainly composed of limestone and divided into three members
(N Fm)Nagri formationPrimary sandstone and minor number of clays
(Q)Quaternary alluviumMainly composed of sand, gravel, silt and clay
Table 2. GCFs used for mapping groundwater potential of the Karak watershed.
Table 2. GCFs used for mapping groundwater potential of the Karak watershed.
GCFsData sourcesFormatProduct
Drainage densityDigital elevation model (DEM) (ASTER 30 m)Digital(Dd)
SlopeDigital elevation model (DEM) (30 m spatial resolution)Digital(SL)
ElevationShuttle Radar Topography Mission (SRTM) data from United States Geological Survey (USGS), resolution: 30 mDigital(EL)
RainfallAnnual rainfall data from Pakistan Meteorological Department (PMD)numbers(RF)
Land use/coverForest Management Center Peshawar (FCMP), KPK, PakistanDigital(LULC)
GeologyGeological map from National Centre of Excellence in Geology (NCEG), University of PeshawarDigital(GEO)
Lineament densityLandsat 8 OLI imagery and Shuttle Radar Topography Mission (SRTM)Digital(LD)
Soil type
GW fluctuation
Directorate General Soil and Water Conservation (DGSC), Khyber Pakhtunkhwa (KPK), Pakistan
Pre-monsoon and post-monsoon groundwater table data (onsite survey)
Digital
Points
(ST)
(GLF)
Table 3. Effect of GCFs, relative weight and score for each GCFs.
Table 3. Effect of GCFs, relative weight and score for each GCFs.
Groundwater Conditioning Factors (GCFs)Major Effect (GCFh)Minor Effect (GCFl)Relative Weights (GCFe + GCFm)Proposed Score of GCFs
Rainfall10.51.506
LU&LC10.5 + 0.5208
Geology1 + 1 + 10.5 + 0.5524
Lineament density1 + 10.52.510
Drainage density1 + 10.52.510
Slope10.5 + 0.5 + 0.52.510
Soil type10.5 + 0.5208
Elevation1 + 10208
GWL fluctuation1 + 10.5 + 0.5316
Total Σ20.5100
Table 4. Classification of weight and ranks of GCFs.
Table 4. Classification of weight and ranks of GCFs.
Groundwater Conditioning Factor (GCF)Categories within the GCFQualitative RankRanksWeight of GCF
Rainfall629–664Very high0606
604–628High04
577–603Moderate03
281–576Low02
13–280Very low 01
Land use/
land cover
Agriculture, Rivers/streamVery high0808
Barren landHigh06
Dam/pondModerate05
Shrub landLow04
Built-upLow–Very low03
ForestVery low02
Range landVery low01
GeologyQAVery high2424
K Fm, N FmHigh18
DP Fm, K FmModerate12
J/T rocksLow08
C Fm, DSVery low04
Lineament density
(km/km2)
1.47–1.78Very high1010
1.15–1.46High08
0.82–1.14Moderate06
0.50–0.81Low04
0.07–0.49Very low 02
Drainage density
(km/km2)
1.08–1.61Very high1010
1.62–1.86High08
1.87–2.11Moderate06
2.12–2.38Low04
2.39–3.08Very high 02
Slope
(degree)
0.0–5.78Flat1010
5.79–13.5Gentle08
13.6–23.4Moderate06
23.5–35.3Steep04
35.4–81.7Very steep02
Soil typeLoamyHigh0608
Loamy clayModerate04
Mainly loamyLow02
Elevation
(meter)
1419High0608
706Moderate04
303Low02
Groundwater level fluctuation (meter)1.57–5.29High1416
5.3–1.08Moderate10
1.09–14.6Low06
14.7–19.3Very low02
Table 5. Error matrix of the GWpot zone-based confusion matrix and Kappa (K) analysis.
Table 5. Error matrix of the GWpot zone-based confusion matrix and Kappa (K) analysis.
S. NoGWpot zonesVery HighHighModerateLowTotalCS 1
1very high0000111
2high120406022416
3moderate040010164
4low100011
Total170407043222
1 CS refer to the correct sample.
Table 6. Average inferred hydro-stratigraphy corresponding to resistivity in the study area. The detailed VES interpretation results are shown in Appendix A.
Table 6. Average inferred hydro-stratigraphy corresponding to resistivity in the study area. The detailed VES interpretation results are shown in Appendix A.
Inferred Hydro-Stratigraphic LithologyInferred Resistivity (ρ) VariationReflection Coefficient VariationThickness Variation (m)
Topsoil954–11090.6723–0.78893.3–4.5
Coarse gravel and sand748–923.10.7313–0.76261.3–4.7
Silty sand mixed lithology323–6730.7841–0.86634.6–8.2
Clayey sand685–1098.30.8911–0.95236.8–28.4
Fine sand and gravels34–98.80.9643–0.975212.6–22.8
Clayey sand with saline water26–84.20.5885–0.743410.8–32.5
Table 7. The following screen schedule is proposed for conversion.
Table 7. The following screen schedule is proposed for conversion.
NoDepth (m) Screen (m)Slot Size (m)
01106.6–113.97.31/12.1–1/15.2
02117.6–128.610.91/12.1–1/15.2
03139.5–150.510.91/12.1–1/15.2
Table 8. Derived borehole lithology-based Normal resistivity logs (NRLs) (short/long) and spontaneous potential logs (SPLs).
Table 8. Derived borehole lithology-based Normal resistivity logs (NRLs) (short/long) and spontaneous potential logs (SPLs).
NoDepth (m) Classified LithologyThickness (m)
010–134Gravel-boulder134
02134–196Gravel-boulder-sandstone62
03196–269Gravel-boulder73
04269–238Sandstone hard59
05238–344Sandstone fine grained16
06344–377Sandstone hard33
07377–383Clay06
08383–429Sandstone hard46
09429–442Sandstone fine grained13
10442–488Sandstone hard46
11488–503Clay15
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

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

AMA Style

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(9):1255. https://doi.org/10.3390/w13091255

Chicago/Turabian Style

Khan, Umair, Haris Faheem, Zhengwen Jiang, Muhammad Wajid, Muhammad Younas, and Baoyi Zhang. 2021. "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 13, no. 9: 1255. https://doi.org/10.3390/w13091255

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop