Groundwater potential zone mapping using GIS and Remote Sensing based models for sustainable groundwater management

Abstract

The present research is conducted in the southern region of Khyber Pakhtunkhwa, Pakistan, to identify groundwater potential zones (GWPZ). We used three models including Weight of Evidence (WOE), Frequency Ratio (FR), and Information Value (IV) with twelve parameters (elevation, slope, aspect, curvature, drainage network, LULC, precipitation, geology, Lineament, NDVI, road, and soil texture, that have been prepared and integrated into ArcGIS 10.8. The reliability of the applied models’ results was validated using Area Under the Receiver Operating Characteristics (AUROC). The GWPZ were reclassified into five classes, i.e. very low, low, medium, high, and very high zone. The area occupied by mentioned classes using WOE are very low (10.14%), low (19.58%), medium (26.75%), high (27.10%), very high (16.40%), while using FR are very low (20.93%), low (32.38%), medium (18.92%), high (13.13%), very high (14.61%) and using IV are very low (14.41%), low (17.17%), medium (29.01%), high (25.85%), and very High (13.53%). The Success Rate Curve of WOE, FR, and IV were 0.86, 0.91, and 0.87, while the Predicted Rate Curve values were 0.89, 0.93, and 0.90, respectively. The results revealed that all applied statistical models performed very well to delineate GWPZ. However, use of the FR technique is strongly encouraged to evaluate the GWPZ, and its findings are especially useful for managing groundwater resources in urban planning. Our approaches for assessing the GWPZ mapping can be applied in any region with similar scenarios and are recommended as a helpful tool for policymakers to manage groundwater.

Keywords:

• AUROC
• bivariate models
• geospatial technology
• groundwater potential zones
• Sentinel-2

1. Introduction

The population of developing countries has increased substantially since the early 1970s and is expected to grow slowly in the ensuing decades (Gerland et al. 2014). Human population growth will presumably necessitate expanded farming operations and as expected, exert a new stress level on groundwater and surface water reservoirs. Excessive water extraction can lead to groundwater level decline, ground subsidence, salinity, decreased well production, enhanced pumping costs, and increased saltwater intrusion (M.Y.A. Khan et al. 2021). In several parts of the world, 60% of significant aquifer’ surfaces were reduced due to all these changes (Khan, ElKashouty, Zaidi, et al. 2023). Moreover, a significant problem of the twenty-first century is climate change, which increases stresses and influences groundwater (Yawar Ali Khan and ElKashouty 2023). Growing demands for water, especially obtaining enough for public demands, is a critical but complicated issue that needs improvement. Emerging water needs, such as the need to obtain enough water to satisfy public expectations, represents another identified as needing improvement (Khan et al. 2022).

Groundwater is the prime source of water storage and provides a means to cleverly help combat water scarcity (Khan, ElKashouty, Subyani, et al. 2023). In many nations, groundwater is a crucial natural resource for drinking water, irrigation, and business supply (Das et al. 2021; Sutradhar et al. 2022). Because of the limited availability of freshwater sources and the overexploitation of most of the world’s groundwater supplies, a severe water crisis is anticipated worldwide by 2025 (Xu et al. 2022; Zhu et al. 2022). Groundwater is more economically significant than surface water because it is easier to access locally, can adjust its fluctuation in response to consumption, provides an extremely consistent source during droughts, and has excellent quality. The dry and semi-arid regions globally, which have seen multiple drought occurrences in recent years due to irregular inadequate precipitation, are especially affected by water scarcity (Venkatesan et al. 2019; Ghosh et al. 2021).

In this water stressed region of Karak, GWPZ mapping is a significant mechanism for agriculture development because rain dependent agriculture is the primary revenue source for the District of Karak. Important crops include wheat, oats, mustard, gramme, and other rainfed varieties. Water scarcity results when there is insufficient water to simultaneously meet the demands of ecosystems and humans. Physical water scarcity essentially creates this basic lack of water. There is also economic water scarcity, or the lack of appropriate infrastructure to enable access to what would otherwise be deemed adequate available water resources. Our current study was conducted in an area deemed water scarce and lacking proper water recharge zonation to recharge the aquifer system. Therefore, the objective of our study was to identify the GWPZ of the region that would ultimately enhance sustainable groundwater management for agriculture and domestic purpose. Particularly, the primary purposes were to recharge an aquifer system to a sustainable level and make the community and help them with long-term water problems.

The GWPZ mapping approach has been used in field investigations in recent decades but can be costly and laborious (Li et al. 2020; Zhuo et al. 2022). Recent technological advances, however, have combined Geographic Information Systems (GIS) and Remote Sensing (RS) to help collect, manipulate, and display various types of data that ultimately generate themed maps (Zhou et al. 2022; Zhou and Liu 2022). Groundwater availability varies among locations depending on a region’s hydrodynamics, rainfall, landscape, geology, soil texture, gradient, etc. (Karimi-Rizvandi et al. 2021). Therefore, these types of parameters are used in a GIS environment to generate the GWPZ.

Scientists worldwide have employed several techniques to understand GWPZs (Islam et al. 2023). Among these, Analytical Hierarchy Process (AHP) (Dong, Zhao, et al. 2023), Logistic Regression (LR) (Dong, Yang, et al. 2023), Frequency Ratio (FR)(Yang et al. 2022), Weight of Evidence (WOE) (Tahmassebipoor et al. 2016), Random Forest (RF) (Prasad et al. 2020), Support Vector Machine (SVM) (Lee et al. 2018), Boosted Regression Tree (BRT) (Naghibi et al. 2016), Linear Discriminant Analysis (LDA) (Lautz et al. 2014), Classification and Regression Tree (CART) (Yin, Wang, T. Li, et al. 2022), Machine learning, Artificial Intelligence and Artificial Neural Network (ANN) (Ahmad et al. 2022). These techniques can generally be divided into qualitative and quantitative variants to provide an extensive GWPZ of a region. Qualitative methods like AHP are utilized frequently because they eliminate the numerical difficulty in decision-making (Abijith et al. 2020). Researchers use the qualitative approach arbitrarily based on their expertise and perspective. A quantitative method is an objective-based approach that generates GWPZ maps by evaluating the dependent and independent variables using mathematical and computational algorithms (Yin, Wang, J. Li, et al. 2022). These GIS-based quantitative methods compute the impacts and association of each class of thematic layers with groundwater (Prasad et al. 2020; Islam et al. 2023).

Herein, we computed the mathematical correlation between groundwater parameters using bivariate statistical techniques including WOE, FR, and IV. The WOE method is the GIS-based quantitative model used for the numerical association of dependent and independent variables to predict GWPZ mapping. This approach has been used for different purposes, including groundwater mapping (Tian et al. 2019). Frequency ratio is a widely used and effective mathematical technique for assessing the probabilistic link between several variables in a geographic manner (Tian et al. 2020), and it also has been used to map GWPZ (Y. Li et al. 2022). The Information Value (IV) technique is a statistical method employed frequently in risk analysis and credit scoring to evaluate the propensity of variables in a dataset (J. Li et al. 2022).

These models are significant for the GWPZ mapping because they compute the numerical-based association between dependent and independent variables. Despite the numerical nature and advantages, these algorithms have yet to be investigated by researchers in our study area. The main aims of this research are to (1) compare the bivariate statistical models to generate GWPZ maps of the study area and (2) provide a scientific means for policymakers and stakeholders to manage the water resources in this region.

2. Materials and Methodology

2.1. Study area

Our study was conducted in District Karak of Khyber Pakhtunkhwa (KPK), Pakistan (Figure 1). Geographically, the Karak is located at latitude and longitude range of 32° 46′ 00′'_ 33° 22′ 00′' North and 70° 43′ 00′' _ 71° 33′ 00′' East (Ahmad et al. 2022). Topographically, the study area is part of the southern Kohat plateau of sub-Himalayas, while climatically, the region is arid to semi-arid (Javed et al. 2019). The climate is dry and non-monsoonal, With cold winters and hot summers. Winter temperatures can reach 0 °C to 15°Cand hot summertime temperatures range from 25 °C to 40 °C. Briefly, in June and July, the temperature may approach 48 °C. The region also receives substantial winter rainfall as clouds move from Afghanistan in a west-to-east direction. The monsoon season produces summertime precipitation, and an annual precipitation ranges from 13 to 37 cm. Rainfall amounts tend to be greatest in summer (July and September) and also parts of winter (December).

Figure 1. Study area map (a) Pakistan association with internal boundaries (b) Study area occurrence in the province of Khyber Pakhtunkhwa (c) Study area with inventory points and elevation in meters.

Geologically, the research area is located in the Southern Kohat Plateau and is influenced by the Himalaya orogenic deformational event (Rehman et al. 2009). The geological age of the regions spans the Plio-Pleistocene to Holocene as indicated by tilting and uplift of the recent terraces along its Southern Topographic front (Bai et al. 2019). Geological records suggest that our study area is bordered to the north by the Karak Fault Zone and to the west by the Bannu Basin; to the east, there are the Shakardara Mountains and the southernmost point is bound by the Shinghar-Surghar ridges (U. Khan et al. 2021). The region’s Siwalik Group is stratigraphically divided into the Lower, Middle, and Upper Siwaliks, respectively (Nazir-Ur-rehman et al. 2020). The Siwalik Group in the region is basically comprised of conglomerates, sandstones, and shales (Nazir-Ur-rehman et al. 2020; Li et al. 2022). The Plateau is the most twisted portion of the Himalayan Foreland Fold and Thrust Belt to the west.

2.2. Datasets

For this research, we used various types of data to produce the thematic layers in GIS environment. The applied data consists of both ground data and satellite imageries. The ground data were acquired from various national research organizations and field visits in the study area. However, the satellite data were acquired from various international open-source platforms to generate parameters for this study. The comprehensive details of both ground and satellite datasets are demonstrated in (Table 1).

Table 1. The comprehensive details of data, its source and purpose.

Data	Resolution	Data Source	Parameter Maps
Sentinel-2	10m	https://dataspace.copernicus.eu/	Land Use Land Cover, NDVI
ALOS DEM	12.5 m	https://search.asf.alaska.edu	Elevation, slope, slope aspect, curvature, and drainage network.
CHIRPS	0.05°	https://chc.ucsb.edu.	To generate precipitation map
Soil	1: 2000000	https://soil.punjab.gov.pk/	To generate soil textural map
Geology	1:650000	https://gsp.gov.pk/	Geological map
Road data		Acquired from KPK highway authority	Road map

2.3. Methodology

We depict details of the approach used in the present study in Figure 2, which consists of the following stages: Initially, the dependent variable (inventory data) map is generated from the ground data (field surveys) and RS data (Sentinel-2, Google Earth). For inventory generation, we first used Google Earth Pro to identify hydrological springs in the study area and exported the KML file of springs to a layer form using KML, and subsequently to a layer tool in GIS platform. In addition to using Google Earth Pro, we used satellite imageries in the Google Earth Engine (GEE) platform using Modified Normalized Difference Water Index (MNDWI). Then, the 12 thematic layers (elevation, slope, aspect, curvature, lineament, drainage, LULC, precipitation, soil, geology, road and NDVI map) were generated using various ground and satellite datasets. Additionally, we computed the associations of both variables. Furthermore, the GWPZ maps were generated using three approaches (WOE, FR and IV). Finally, the accuracy assessment of all algorithms has been performed using the AUROC technique.

Figure 2. The detailed methodological flowchart of the present study.

2.3.1. Inventory map of the study area

Several previous studies have employed the precise positions of springs, wells, and reservoirs as inventory data for the GWPZ mapping (Prasad et al. 2020). In our research here, both ground and satellite data were considered for the analysis to produce GWPZ maps of the study area. One hundred sixteen inventory points across the area of interest were detected using Google Earth Pro, Senitnel-2 data, validated during field visits, and generated the updated inventory of District Karak (Figure 1). Furthermore, the 116 points-based inventory data were divided into two inventory datasets, i.e. training and testing datasets, using a subset tool in ArcGIS 10.8. The training and testing class consisted of 80% and 20% data. For the training, 80% and 20% of the data were used for modelling and applied for validation, respectively.

2.3.2. Thematic layers preparation for GWPZ mapping

It is crucial to employ the appropriate criteria when generating a GWPZ map for a specific region (Naghibi et al. 2018). Based on previous literature (Rahmati et al. 2016; Golkarian et al. 2018; Bui et al. 2019; Chen et al. 2019; Jaafarzadeh et al. 2021) and our field observations, we considered 12 parameters as influential (elevation, slope, aspect, curvature, lineament, drainage, LULC, precipitation, soil, geology, road and NDVI map) to generate a GWPZ map of the study area using a GIS platform. For modelling purposes, each parameter was reclassified and normalized to the spatial resolution of 12.5 m. Details of these parameters used in this research follow.

2.3.2.1. Elevation

Elevation is the potential indicator and it can high influence the approach to GWPZ mapping (Naghibi et al. 2018). Elevation might directly affect rainfall, the amount, frequency, and process of localized precipitation, and thus the recharging of groundwater (Berhanu and Hatiye 2020). The elevation map of the current study was generated from ALOS DEM with a 12.5 m resolution. Additionally, the generated elevation map was reclassified into five classes, i.e. <500 m, 500–700 m, 700–900 m, 900–110 m and >1100 m (Figure 3a).

Figure 3. The thematic layer of the present study.

2.3.2.2. Slope

Slope angle is considered the most significant topographic factor affecting groundwater potential (Rahmati et al. 2015; Lee et al. 2018). Aquifers that store groundwater are replenished more frequently under gentle slopes because they allow rainwater to infiltrate and recharge beneath the surface. Infiltration may be impeded by steeper slopes, increasing surface runoff and decreasing recharge (Maskooni et al. 2020). In normal circumstances, more groundwater percolates into the earth, surface water flow velocity is slower, and slopes are less steep. However, the steeper the slope the greater the surface runoff, which leads to reduced groundwater percolation. The slope map was generated in the GIS platform from ALOS DEM, and the slopes were categorized into five classes using the natural break algorithm (Figure 3b).

2.3.2.3. Aspect

Aspect is a fundamental variable for successful GWPZ mapping. Aspect designates the slope’s orientation, which is determined by the exposure to sunshine, wind, and precipitation (Pande et al. 2021). Aspect also simulates the vegetation, moisture, air, rock bedding’s structure, and moisture retention, all of which can impact the physical components of soil structure, including sloped materials and possible groundwater zones. Moreover, particularly in dry, arid regions devoid of vegetation, the slope aspect substantially influences hydrological parameters such as evapotranspiration, precipitation direction, and weather process (Ercanoglu and Gokceoglu 2002; Benjmel et al. 2022). The aspect map was prepared to identify the groundwater possession at various slope orientations. The aspect map was generated from ALOS DEM and reclassified into nine classes using the natural break algorithm in the ArcGIS platform (Figure 3c).

2.3.2.4. Curvature

Curvature is a significant topographic parameter for GWPZ mapping because it influences surface water and subsurface aquifer systems (Regmi et al. 2015). The curvature of peak slopes, which represents the topography of the area, is the most important feature to measure to observe possible groundwater areas, gauge water hydrology, and estimate terrain instability. Curvatures provide greater precision in describing the dynamics of subsurface hydrology, soil formation, and soil accumulation. The primary functions of curvature are to influence flow convergence and divergence, and flow acceleration and deceleration (Ercanoglu and Gokceoglu 2002; Al-Abadi et al. 2016). The curvature map of the current research is computed from ALOS DEM having 12.5 m resolution and reclassified into three classes, i.e. concave, flat and convex (Figure 3d). Convex terrains, characterized by elevated topography, often create conditions favourable for groundwater occurrence. The collection of rainwater and enhanced infiltration in concave terrains, featuring lower elevations, promotes groundwater recharge. Both terrains contribute to groundwater availability, but convex terrains, with their potential for increased recharge, are particularly conducive to groundwater occurrence.

2.3.2.5. Proximity to drainage network

Drainage networks are directly associated with surface runoff and are inversely associated with groundwater storage (Magesh et al. 2012). Drainage density, an important metric for assessing the infiltration rate of a location, can significantly impact groundwater potential (Arulbalaji et al. 2019). High drainage density values promote runoff and hence signify a zone with little groundwater potential (Muhammad and Khalid 2017). Conversely, low drainage density denotes high infiltration, which can significantly raise the potential of groundwater (Prasad et al. 2008; Arulbalaji et al. 2019). The drainage network of the Karak was extracted from the ALOS DEM and five buffers of <200 m, 200–400 m, 400–600 m, 600–800 m and >800 m were applied on the drainage lines using ArcGIS 10.8 (Figure 3e).

2.3.2.6. Precipitation

Precipitation is considered the most influential parameter for GWPZ mapping. The period and magnitude of precipitation significantly influences the surface flow, recharge, and infiltration (Magesh et al. 2012). Rainfall is the primary source of surface and groundwater in this region, thus the volume of the basin’s recharge is greatly influenced by the amount and distribution of that rainfall (Zahir et al. 2020). The precipitation map of the present study is computed from the CHIRPS data of 2010 to 2023 using a 33-kernel machine learning algorithm in Google Earth Engine (GEE). The rainfall map was finally reclassified into five classes (Figure 3f).

2.3.2.7. LULC

The LULC of any area governs the transportation of groundwater. LULC is significantly responsible for controlling evapotranspiration, surface water runoff, and replenishment of aquifers (Karimi-Rizvandi et al. 2021). Groundwater is crucial for regulating land usage and land cover in an area, and it is also vital for the establishment and accessibility of this resource. It controls several hydrogeological processes, such as surface runoff, evapotranspiration, and infiltration. In addition, groundwater offers valuable information on soil conditions, such as surface and groundwater levels, soil moisture content, and potential indications of groundwater availability (Nguyen et al. 2020). Deeper plant roots loosen the ground and cause rocks to collapse, increasing water penetration and subsequently raising groundwater levels (Fatema et al. 2023). In the urbanized part of our study area, the surface of the region is impervious which leads to surface runoff and reduced infiltration. In the crop cultivated parts of the region, the land is suitable for water infiltration. The LULC map of the study area was generated using machine learning techniques, i.e. Random Forest (RF) in GEE. The prepared map in GEE was exported into ArcGIS 10.8 for further analysis. The LULC map of the present study consists of 7 classes: tree, shrubland, grassland, cropland, built-up, barren, and water (Figure 3g).

2.3.2.8. Geology

Lithology is the general physical characteristics of rocks and its formations affects groundwater and regulates the aquifer’s permeability and porosity because of its conductivity and penetrability (Islam et al. 2023). Geology is the most significant determinant of groundwater potential because penetration is mostly determined by the permeability or conductivity of various rock types (Bera et al. 2020). Information on the subsurface geological formations of our study area is provided by the geological map. These material properties impact groundwater presence, storage, and transportation (Muavhi et al. 2022). The lithological map of the current research was digitized and prepared from the Northern map of Pakistan using ArcGIS 10.8 (Figure 3h).

2.3.2.9. Lineament

Lineaments are tectonic features that are linear, rectilinear, and curvilinear. They are easily visible in satellite imagery and may be used to characterize a variety of attributes, including faults, master joints, cover of vegetation, topographic linearity and formation, roads, bridges, valleys, straight rivers, and lines between different lithological units (Berhanu and Hatiye 2020). Secondary porosity is considered just as important as primary porosity in determining hydrogeological conditions. Lineaments are linear structures of tectonic origin that symbolize secondary porosity and are critical for hydro-geological regulation of groundwater flow and storage (Magesh et al. 2012). Lineaments can be located on the globe’s surface in a linear or curved pattern using DEM, and fieldwork data (Magesh et al. 2012; Lee et al. 2018). The lineament map of the present study was calculated from ALOS DEM hill shades and field survey in the study area, and the buffers were applied on lineament features (Figure 3i).

2.3.2.10. Soil map

Soil mapping is also a significant factor for GWPZ mapping. Soil composition influences water storage because soil characteristics govern the area’s permeability (Karimi-Rizvandi et al. 2021). Generally, soil texture affects the soil’s moisture content, infiltration rate, hydraulic conductivity, soil permeability, grain size, and unique composition, all of which can significantly affect recharge potentiality (Lee et al. 2017). Our study area is comprised of three types of soil groups including loamy clay, loamy, and mainly loamy (Figure 3j).

2.3.2.11. Normalized Difference vegetation Index (NDVI)

The NDVI has a considerable impact on both groundwater retention and availability. A higher NDVI value indicates dense coverage of plants because vegetation both slows runoff and assists in recharging the groundwater reservoir (Hasanuzzaman et al. 2022). The amount of rainfall, temperature, soil adoption, vegetation water content, soil humidity, and evaporation are all related to the change in depth to groundwater. Because evapotranspiration rises during the dry months and precipitation and groundwater recharge in fall, the soil’s fresh water capacity diminishes which may cause related drought stress and restricted growth (Pourtaghi and Pourghasemi 2014). The NDVI map of the study (Figure 3k) was prepared in GEE and exported to ArcGIS for further analysis.

2.3.2.12. Road distance

A network of roads can adversely impact groundwater. During road construction, the surface porosity and permeability within the strata are reduces because of the lithification phenomenon (Islam et al. 2023). Road structure facilitates surface runoff, which in turn are reduced infiltration rate. According to laboratory permeability measurements, loading caused by road embankments greatly decreases the permeability of the alluvial sediments, resulting in an impermeable groundwater system boundary. Numerical modelling in our study demonstrated that groundwater levels rose in the upstream region but decreased downstream when a road embankment served as an impermeable boundary. Furthermore, the intensity of this impact can change with the distance from the road embankment, while the radius of influence of road embankment construction on the upstream area groundwater supplies is determined by the geological and hydraulic conditions (Kim et al. 2023). The road map was prepared from the Khyber Pakhtunkhwa Highway authority data, and five buffers were applied in ArcGIS 10.8 (Figure 3l).

3. GWPZ mapping methods

In the present study, three GIS-based bivariate models (WOE, FR and IV) were used to compute the relationship between the dependent variable (inventory data) and independent parameters (groundwater influencing parameters). These models were used to produce GWPZ mapping for the study area. The overview of the models is as follows.

3.1. Weight of Evidence (WOE)

The WOE is a numerical, ‘data-driven’ Bayes rule-based algorithm that forecasts the likelihood of occurrences (Tahmassebipoor et al. 2016). The predictive ability of an independent variable with respect to the dependent variable is indicated by WOE. This approach previously was used to generate maps and investigate possible gold dispersion (Kim et al. 2018). Since then, the WOE method has been widely used in studies of landslides. Compared to field observation and other deterministic approaches, this method is currently considered dependable and economical (Naghibi and Pourghasemi 2015). We used this approach in our study to determine the prospect of groundwater (Corsini et al. 2009). The WOE approach was initially considered to assess mineral prospective mapping using GIS-based modelling (Cheng et al. 2011). The mathematical formulation of the WOE model is described in depth in Bonham-Carter (Bonham-Carter 1994). This approach determines the weight for each class of groundwater precondition parameter’s existence or absence (M or M*) depending on the presence or absence of the inventory pixels (N or N*) as mentioned in Bonham–Carter Equations (1)(1) $${\mathrm{W}}^{+}=\mathrm{\text{ln}}\frac{\mathrm{P}\left(\frac{\mathrm{M}}{\mathrm{N}}\right)}{\mathrm{P}\left(\frac{\mathrm{M}}{\text{N*}}\right)}$$(1) and (2)(2) $$\mathrm{W}\_=\mathrm{\text{ln}}\frac{\mathrm{h}\left(\frac{\mathrm{M}\mathrm{*}}{\stackrel{⃛}{\mathrm{N}}}\right)}{\mathrm{h}\left(\frac{\overline{\mathrm{M}\mathrm{*}}}{\stackrel{⃛}{\mathrm{N}\mathrm{*}}}\right)}$$(2) (Keller 1995) as follow: (1) $${\mathrm{W}}^{+}=\mathrm{\text{ln}}\frac{\mathrm{P}\left(\frac{\mathrm{M}}{\mathrm{N}}\right)}{\mathrm{P}\left(\frac{\mathrm{M}}{\text{N*}}\right)}$$(1) (2) $$\mathrm{W}\_=\mathrm{\text{ln}}\frac{\mathrm{h}\left(\frac{\mathrm{M}\mathrm{*}}{\stackrel{⃛}{\mathrm{N}}}\right)}{\mathrm{h}\left(\frac{\overline{\mathrm{M}\mathrm{*}}}{\stackrel{⃛}{\mathrm{N}\mathrm{*}}}\right)}$$(2) where P is the probability while M is the presence of groundwater parameters, M* is the absence of groundwater pixels, N is the presence of inventory, and N* is the absence of a well. When a component is present or absent, the weights of the evidence are positive weight (W⁺) and negative weight (W^-), respectively Equation (2)(2) $$\mathrm{W}\_=\mathrm{\text{ln}}\frac{\mathrm{h}\left(\frac{\mathrm{M}\mathrm{*}}{\stackrel{⃛}{\mathrm{N}}}\right)}{\mathrm{h}\left(\frac{\overline{\mathrm{M}\mathrm{*}}}{\stackrel{⃛}{\mathrm{N}\mathrm{*}}}\right)}$$(2) (Corsini et al. 2009).

The magnitude of these comparisons represents the entire spatial relationship between each class of parameter and the groundwater. Positive and negative contrast values imply favourable and unfavourable spatial associations, respectively. The comparison can be computed using Equations (3)(3) $$L=W+‐\mathrm{ }W\_$$(3) and (4)(4) $$\mathit{\text{Wmap}}=\sum L$$(4) (Tahmassebipoor et al. 2016). (3) $$L=W+‐\mathrm{ }W\_$$(3) (4) $$\mathit{\text{Wmap}}=\sum L$$(4)

The Final GWPZ delineation is computed through Equation (5)(5) $$\mathrm{\text{GWPZ}}=\sum \mathit{\text{Wmap}}$$(5) (Getachew and Meten 2021). (5) $$\mathrm{\text{GWPZ}}=\sum \mathit{\text{Wmap}}$$(5)

3.2. Frequency Ratio model (FR)

A FR model is the most frequently used bivariate statistical techniques for research on hazards associated with nature, such as floods, landslides, and droughts (Sarkar et al. 2022). In this study FR, a GIS-based bivariate method, was established to identify the GWPZ by analysing the associations between the influencing parameters (Oh et al. 2011; Naghibi et al. 2015). This algorithm has been used to identify the numerical relationship between dependent variables and independent. The Equation (6)(6) $$\mathrm{\text{FR}}=\frac{\mathrm{R}/\mathrm{\text{TR}}}{\mathrm{\text{AP}}/\mathrm{\text{TA}}}$$(6) (Ozdemir 2011) has been used to compute the FR: (6) $$\mathrm{\text{FR}}=\frac{\mathrm{R}/\mathrm{\text{TR}}}{\mathrm{\text{AP}}/\mathrm{\text{TA}}}$$(6) where FR = Frequency, R is the number of pixels from each thematic map that are well connected, and TR is the total number of pixels in the region of interest. Pixel counts for each thematic map and the region of concern are represented by AP and TA, respectively. Greater ratios typically signify stronger correlations between conditioning factors and dependent factors, and vice versa. The average value of an area’s groundwater pixels throughout the entire area is 1. An FR value >1 indicates high probability of groundwater pixel occurrence. Comparatively, values of <1 indicate little relationship in the probabilities of groundwater pixel occurrence. The raster calculator tool in ArcGIS 10 was used to integrate all the reclassified parameters to estimate the GWPZ in our current research using the Equation (7)(7) $${\mathit{\text{GWPZ}}}_{\mathit{\text{FR}}}={\mathit{\text{FR}}}_E+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_A+{\mathit{\text{FR}}}_L+{\mathit{\text{FR}}}_C+{\mathit{\text{FR}}}_D+{\mathit{\text{FR}}}_R+{\mathit{\text{FR}}}_{\mathit{\text{LULC}}}+{\mathit{\text{FR}}}_{\mathit{\text{linement}}}+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_{\mathit{\text{NDVi}}}+{\mathit{\text{FR}}}_{\mathit{\text{Road}}}$$(7) (Guru et al. 2017). (7) $${\mathit{\text{GWPZ}}}_{\mathit{\text{FR}}}={\mathit{\text{FR}}}_E+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_A+{\mathit{\text{FR}}}_L+{\mathit{\text{FR}}}_C+{\mathit{\text{FR}}}_D+{\mathit{\text{FR}}}_R+{\mathit{\text{FR}}}_{\mathit{\text{LULC}}}+{\mathit{\text{FR}}}_{\mathit{\text{linement}}}+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_{\mathit{\text{NDVi}}}+{\mathit{\text{FR}}}_{\mathit{\text{Road}}}$$(7) where ${\mathit{\text{GWPZ}}}_{\mathit{\text{FR}}}$ is the groundwater potential zone using FR while $E,S,A,L,C,D,R,\mathit{\text{LULC}},\mathit{\text{Linea}},S,\mathit{\text{NDVI}},R$ in Equation (7)(7) $${\mathit{\text{GWPZ}}}_{\mathit{\text{FR}}}={\mathit{\text{FR}}}_E+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_A+{\mathit{\text{FR}}}_L+{\mathit{\text{FR}}}_C+{\mathit{\text{FR}}}_D+{\mathit{\text{FR}}}_R+{\mathit{\text{FR}}}_{\mathit{\text{LULC}}}+{\mathit{\text{FR}}}_{\mathit{\text{linement}}}+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_{\mathit{\text{NDVi}}}+{\mathit{\text{FR}}}_{\mathit{\text{Road}}}$$(7) represent elevation, slope, aspect, curvature, drainage network, rainfall, LULC, lineament, soil, Normalized difference vegetation index and road, respectively. The final map using the FR process is presented in Figure 4b.

Figure 4. GWPZ mapping using (a) WOE, (b) FR, and (c) IV.

3.3. Information value (IV)

The IV model is a bivariate statistical technique used to forecast spatial association between groundwater inventory data and groundwater parameters (Yin and Yan 1988). Yin and Yan designed this approach, and subsequently modified by Sarkar et al. (2013). The estimated IV helps identify each independent variable class’s contribution to the dependent variables’ occurrence. The IV for each class of a factor map in our study were established based on the existence of groundwater in a particular map unit. The same coordinate system was used to convert each factor map into a raster map and then were divided into distinct classes despite having the same pixel size (12.5 m × 12.5 m). Considering their spatial distribution, the GWPZ map was randomly split into training wells and springs and testing data. Using ArcGIS software, the rasterized training GWPZ map was superimposed over the rasterized groundwater factor maps, and the information value model was used to determine the information values for each factor map’s classes.

The log ratio of conditional probability to prior probability was used to determine the IV of a specific parameter class (Wubalem and Meten 2020). The conditional probability was calculated by dividing the groundwater pixels in a single parameter class by the pixel of a subclass of groundwater parameters. However, the prior probability was estimated by dividing the entirety of groundwater pixels of our study area by the total pixels in the entire region (Equation (8)(8) $$\mathrm{\text{IV}}=\mathrm{\text{log}}\frac{\frac{\mathrm{\text{NSLpix}}}{\mathrm{\text{Ncpix}}}}{\mathit{\text{Ntspix}}/\mathit{\text{Ntapix}}}$$(8) ). (8) $$\mathrm{\text{IV}}=\mathrm{\text{log}}\frac{\frac{\mathrm{\text{NSLpix}}}{\mathrm{\text{Ncpix}}}}{\mathit{\text{Ntspix}}/\mathit{\text{Ntapix}}}$$(8)

The NSLpix is the number of groundwater pixels in a given class, $\mathit{\text{Ncpix}}$ is the number of pixels in the class, Nts $\mathit{\text{pix}}$ is the total number of groundwater pixels, and $\mathit{\text{Ntapix}}$ is the total number of pixels in the study area. In this study, all parameters were reclassified, rasterized, normalized to 12.5 m resolution and converted to WGS_1984, UTM Zone, 43 N projected system.

The Equation (9)(9) $$\mathit{\text{GWPZ}}={\sum }_{i=1}^M{X}_{\mathit{\text{ji}}}\mathrm{*}\mathit{\text{Ij}}$$(9) (Farooq and Akram 2021) is used to determine GWPZ by summing each variable category. (9) $$\mathit{\text{GWPZ}}={\sum }_{i=1}^M{X}_{\mathit{\text{ji}}}\mathrm{*}\mathit{\text{Ij}}$$(9) where, M = total number of classes, and $X_{\mathit{\text{ji}}}$= 1 if i class present in parameter j and 0 if i class is absent in variable j (see Figure 4c for the IV generated map). A class’s weight for each factor was determined by dividing the density of springs and wells by the total area’s density of wells and springs; alternatively, the information value could be used to determine the probability of groundwater in each class and overall. The factor classes with the greatest likelihood of groundwater potential are those whose IV value exceed 0.1. Factor classes with negative values, in contrast, suggest the presence of a component that has no appreciable impact on groundwater potential.

3.4. Validation

Validation is the most significant process for validating the accuracy and reliability of the data used in models. Several techniques have been used to evaluate model performance including common approaches of Area Under the Curve Receiving Operating Characteristic (AUCROC), graphic analysis using Success Rate Curve (SRC) and Predicted Rate Curve (PRC). The values associated with these methods vary from 0.5 to 1.0, with values close to 1.0 denoting the greatest degree of accuracy; values close to 0.5 indicate model inaccuracy (Das and Pardeshi 2018). The following classifications distinguish the link between the AUC value and the precision of the predictions: 0.9–1.0 (excellent), 0.8–0.9 (very good), 0.7–0.8 (good), 0.6–0.7 (average) and 0.5–0.6 (poor) (Naghibi et al. 2015). SRC is computed using the training dataset, which has been used to run and train the model for GWPZ. Therefore, SRC alone is not a reliable technique to assess any model’s capacity for prediction (Zabihi et al. 2016). Hence, Predicted Rate Curve is used based on testing data for the prediction. The PRC represents the model’s ability to predict GWPZ as well as its predictor factors (Lee and Pradhan 2007).

4. Results

Understanding how an area’s inventory and GWPZ mapping-effective parameters relate is important. Several predisposing factors can affect groundwater areas’ potentiality (Mukherjee et al. 2021). In the present research, we applied 12 parameters including elevation, slope, aspect, curvature, drainage network, LULC, precipitation, geology, NDVI, road and soil texture, using ArcGIS 10.8, utilizing ground data and satellite imageries. We also developed three GIS-based techniques to produce GWPZ for the research region. The association of parameters with inventory data is included in Table 2.

Table 2. The statistical analysis of dependent and independent variables.

Parameters	Class	No of pixels in class	No of Groundwater pixels in a class	W^+	W^-	WOE	% Pixels in Class	% Groundwater pixels in Class	(FR)	IV = log (A/B)
Elevation	<500	3818355	823	−2.04	0.23	−2.27	22.58	2.93	0.13	−2.04
	500–700	7446408	16709	0.30	−0.32	0.63	44.04	59.52	1.35	0.30
	700- 900	3618070	7914	1.36	−1.55	2.91	21.40	28.19	1.32	2.36
	900–1100	1560886	2125	−0.20	0.02	−0.22	9.23	7.57	0.82	−0.20
	>1100	463101	502	−0.43	0.01	−0.44	2.74	1.79	0.65	−0.43
Slope	<10	7875182	15864	0.20	−0.21	0.41	46.43	56.5098	1.22	0.20
	10–20	6020279	10808	0.08	−0.05	0.13	35.49	38.4996	1.08	0.08
	20–30	2002041	1240	−0.07	0.00	−0.07	11.80	4.41706	0.37	−0.07
	30–40	799850	131	−3.23	0.12	−3.35	4.72	0.46664	0.10	−3.23
	>40	265807	30	−2.69	0.01	−2.70	1.57	0.10686	0.07	−2.69
Aspect	Flat	2125998	3269	−0.07	0.01	−0.08	12.53	11.64	0.93	−0.07
	Northeast	1577145	4147	0.46	−0.06	0.53	9.30	14.77	1.59	0.46
	East	1461113	2551	0.05	−0.01	0.06	8.61	9.09	1.05	0.05
	Southeast	1542041	3467	0.31	−0.04	0.34	9.09	12.35	1.36	0.31
	South	2069026	2946	−0.15	0.02	−0.17	12.20	10.49	0.86	−0.15
	Southwest	1656126	3250	0.17	−0.02	0.19	9.76	11.58	1.19	0.17
	West	1878085	2237	−0.33	0.03	−0.36	11.07	7.97	0.72	−0.33
	Northwest	2046900	3086	−0.09	0.01	−0.11	12.07	10.99	0.91	−0.09
	North	2606725	3120	−0.32	0.05	−0.37	15.37	11.11	0.72	−0.32
Curvature	Concave	3276810	11158	0.7	−0.3	1.01	19.38	39.75	2.05	0.72
	Flat	7232284	12487	0.0	0.0	0.07	42.78	44.48	1.04	0.04
	Convex	6397726	4428	−0.9	0.3	−1.18	37.84	15.77	0.42	−0.88
Distance to Stream	<200	841329	3637	0.96	−0.09	−1.12	4.961	12.956	0.77	−0.34
	200 − 400	818268	3493	0.95	−0.08	0.51	4.825	12.443	1.61	0.47
	400–600	831738	2671	0.66	−0.05	0.71	4.904	9.514	1.94	0.66
	600–800	802851	2134	0.47	−0.03	1.03	4.734	7.602	2.58	0.95
	>800	13665032	16138	−0.34	0.79	1.05	80.57	57.486	2.61	0.96
Precipitation (mm/year)	<700	3669506	507	−2.49	0.23	−2.71	21.7	1.8	0.08	−2.47
	700–800	3603654	1629	−1.30	0.18	−1.48	21.3	5.8	0.27	−1.30
	800–900	4522943	6282	−0.18	0.06	−0.23	26.7	22.4	0.84	−0.18
	900–1000	2424621	7209	0.59	−0.14	0.73	14.3	25.7	1.79	0.58
	>1000	2724869	12446	1.02	−0.41	1.43	16.1	44.3	2.76	1.01
LULC	Shrub	4739024	937	−2.13	0.29	−2.42	27.92	3.34	0.12	−2.12
	Trees	402428	6222	2.25	−0.23	2.48	2.37	22.16	9.35	2.24
	Barren	7043895	12396	0.06	−0.05	0.11	41.49	44.16	1.06	0.06
	Grassland	2233274	2243	−0.50	0.06	−0.56	13.16	7.99	0.61	−0.50
	Crops	2060647	6046	0.57	−0.11	0.69	12.14	21.54	1.77	0.57
	Urban	482231	309	−0.95	0.02	−0.97	2.84	1.10	0.39	−0.95
	Water	14309	11	−0.77	0.00	−0.77	0.08	0.04	0.46	−0.77
Lithology	Siwalics	20656	36	−0.30	0.04	−0.34	12.73	9.47	0.74	−0.30
	Quaternary	7565	131	0.55	−0.20	0.75	20.06	34.47	1.72	0.54
	Rawalpindi	16554	175	0.05	−0.04	0.09	43.89	46.05	1.05	0.05
	Pal	8797	38	−0.85	0.16	−1.01	23.32	10	0.43	−0.85
Lineament Buffer	<500	462160	1927	0.84	−0.04	0.88	2.97	6.86	2.31	0.84
	1500	482907	1597	0.69	−0.03	0.72	2.86	5.69	1.99	0.69
	3000	535801	1341	0.51	−0.02	0.53	2.87	4.78	1.66	0.51
	5000	593403	1005	0.20	−0.01	0.21	2.92	3.58	1.23	0.20
	>5000	14941848	22203	−0.11	0.59	−0.70	88.38	79.09	0.89	−0.11
Soil	Loamy clay	2102769	318	−2.39	0.12	−2.51	12.41	1.132	0.091	−2.39
	Loamy	9321208	27239	0.56	−2.7	0.51	2005	783.63	0.73	0.42
	Mainly Loamy	5518421	516	−2.87	0.37	0.43	1187	14.84	0.59	0.34
NDVI	Low	3954806	1663	−1.37	0.21	−1.58	23.40	14087.58	−0.07	−1.37
	High	12947448	26410	0.21	−1.37	0.22	76.60	46120.64	0.53	0.21
Distance to Road	<200	502160	139	−1.79	0.03	−1.82	2.97	0.50	0.17	−1.79
	200–400	482907	245	−1.19	0.02	−1.21	2.86	0.87	0.31	−1.19
	400–600	485801	283	−1.05	0.02	−1.07	2.87	1.01	0.35	−1.05
	600–800	493403	305	−0.99	0.02	−1.01	2.92	1.09	0.37	−0.99
	>800	14941848	27101	0.09	−1.21	0.30	88.38	96.54	1.09	0.09

4.1. Application of WOE for GWPZ mapping

The GWPZ was generated using WOE technique and the output map was further reclassified into five classes i.e. very low, low, medium, high and very high (Jenks 1967) natural break algorithm in ArcGIS. In this study the classes occupied 262.09 km² (10.14%), 506.00 km² (19.58%), 691.10 km² (26.75%), 700.18 km² (27.10%), 423.75 km² in (16.40%), respectively (Table 3).

Table 3. Showing the results of WOE, FR and IV in area km² (% age).

Models	Classes	km^2	(% age)
WOE	Very Low	262.09	10.14
	Low	506.003	19.58
	Medium	691.10	26.75
	High	700.18	27.10
	Very High	423.75	16.40
FR	Very Low	538.90	20.93
	Low	833.53	32.38
	Medium	487.14	18.92
	High	337.94	13.13
	Very High	376.07	14.61
IV	Very Low	372.46	14.41
	Low	443.67	17.17
	Medium	749.52	29.01
	High	667.97	25.85
	Very High	349.50	13.53

All factors of the WOE technique were computed for each groundwater conditioning variable (Figure 4a and Table 2). The spatial relationships between the classes of each conditioning factor and groundwater occurrence are also shown (Table 2), together with the numeric value of each applicable variable. A positive or negative spatial correlation indicates a positive or negative contrast value for the WOE model.

We determined that the elevations of our study region consisted of five classes, i.e. <500 m, 500–700 m, 700–900 m, 900–1100 m and >1100 m (Table 2). The WOE value for these classes was −2.27, 0.63, 2.91, −0.22 and −0.44, respectively. The results indicated that the greatest groundwater potential zone occurred at 700–900 m, followed by 500–700 m, while the lowest value occurred in elevations <500 m and then those >1100 m. Slopes <10° are the best class for GWPZ, followed by 10–20°, while 30–40° and those >40° demonstrate unreliable associations with inventory data. The WOE for <10° was 0.41, and the association between WOE for >40° was −2.70. These results illustrated a strong positive association between variables using WOE for NE when the value was 0.53 and suggests that NE is the most significant class variable for our study area when using the WOE technique. We also determined that a concave structure is the most feasible class for groundwater potential, while a convex structure did not influence groundwater (Table 2). The results suggest that the <200 m drainage class was the least significant class, whereas the >800 m class was influential when WOE was −1.12 and 1.05, respectively. The precipitation factor was subdivided into five classes i.e. <700 mm/year, 700–800, 800–900, 900–1000 and >1000mm/year. The greatest rainfall class i.e. >1000mm/year, demonstrated the most positive correlation between both variables, while the low precipitation class <700 mm/year, had the least association between variables (Table 2). The WOE for >1000mm/year was 1.43, followed by 900–1000 mm/year, while <700mm/year has the least association of −2.71. The LULC parameter of this study was classified into seven classes having different associations with dependent variables. We found that forest class strongly associated with the dependent variable, as this class had 2.37 for WOE (Table 2). We found that urban and grassland regions had no influential role or association. However, the lithology parameter significantly influenced the GWPZ. The most significant class associated with this factor was Quaternary rocks, with a WOE model value of 0.75, followed by the Rawalpindi group. The lineament predisposing factor had five classes (Table 2), and the WOE <500 m class was 0.88, followed by 1500 m. Simultaneously, the >5000 m class had no significant association of −0.70. The WOE for the Loamy class was 0.51 (Table 2), suggesting that loamy clayey is the least significant class of this parameters in this study. We determined that the NDVI is an essential signal for subsurface water perspective (Table 2). The WOE for the high NDVI class was 0.22. Lastly, the road parameter in this study was categorized into five buffer zones, i.e. <200 m, 200–400 m, 400–600 m, 600–800 m and >800 m classes. Results concluded that road surface restricts water infiltration and adversely impacts GWPZ.

4.2. Application of FR for GWPZ mapping

The FR model was used to generate GWPZ maps and then reclassified into potential zones ranging from very low, low, medium, high, and very high using Natural break algorithm in GIS environment. The zones that we studied contained the following values: very low zone (538.90 km², or 20.93%), the low zone (833.53 km², or 32.38%), the medium (487.14 km², or 18.92%), the high (337.94 km², or13.13%), and the very high zone (376.07 km², or, 14.61%) (Table 3).

Values reported in Table 2 indicate each factor’s FR value calculations using the number of pixels and wells in each contributing variable. The FR is calculated using Equations (6)(6) $$\mathrm{\text{FR}}=\frac{\mathrm{R}/\mathrm{\text{TR}}}{\mathrm{\text{AP}}/\mathrm{\text{TA}}}$$(6) and (7)(7) $${\mathit{\text{GWPZ}}}_{\mathit{\text{FR}}}={\mathit{\text{FR}}}_E+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_A+{\mathit{\text{FR}}}_L+{\mathit{\text{FR}}}_C+{\mathit{\text{FR}}}_D+{\mathit{\text{FR}}}_R+{\mathit{\text{FR}}}_{\mathit{\text{LULC}}}+{\mathit{\text{FR}}}_{\mathit{\text{linement}}}+{\mathit{\text{FR}}}_S+{\mathit{\text{FR}}}_{\mathit{\text{NDVi}}}+{\mathit{\text{FR}}}_{\mathit{\text{Road}}}$$(7) (Table 2). To define GWPZ, 12 parameters that affect groundwater occurrence were combined in the ArcGIS environment. The computed GWPZ using the FR technique was further categorized into five classes, i.e. very low, low, medium, high and very high (Figure 4b). The association of elevation class with inventory data using FR revealed that the greatest potential zone occurred at 500–700 m, followed by 700–900, while the least value occurred in elevations <500 and followed by >1100 m. The FR value for <10° and >40° were 1.22 and 0.07, respectively. The FR for concave and convex structures were 2.05 and 0.42, respectively. The drainage network was categorized into five classes, each with a different association with the dependent variable. The results revealed that the >800 m class was most significant, while <200 m had no influential impact. The WOE for >1000mm/year was 1.43 and for <700mm/year, −2.71. The results demonstrate that <500 m lineament buffer is the most significant class for GWPZ, strongly associating with the dependent variable, while >5000 m class has no significant association with inventory data. Our findings demonstrate that NDVI is an essential sign for groundwater presence (Table 2). The FR for the high NDVI class was 0.53. The results conclude that road surface restricts water infiltration and adversely impacts GWPZ.

4.3. Application of IV for GWPZ mapping

The IV technique is a theory largely applied in credit scoring and risk assessment to examine the ability of various factors to forecast various outcomes. The five classes of final map were generated using Natural Break Algorithm in ArcGIS 10.8. The area of 372.46 Km² (14.41%), 443.67 Km² (17.17%), 749.52 Km² (29.01%), 667.97 Km² (25.85%), and 349.50 Km² (13.53%) were occupied by very low, low, medium, high, and very high zone class of GWPZ map of IV model (Table 3). The statistical results associated with the dependent and independent variables are also provided (Table 2). Results revealed that elevation consisted of five classes, i.e. <500, 500–700, 700–900, 900–1100 and >1100, and the IV value for mentioned classes area were −2.04, 0.30, 2.36, −0.20 and −0.43. Results suggested that the greatest potential zone occurred at 500–700 m, followed by 700–900, while the least value occurred in elevation <500 m and followed by >1100 m. The GWPZ used IV and revealed that slope <10° was the best class for GWPZ, followed by 10–20°, while >40° and 30–40° had deleterious association with inventory data. The IV for <10° are 0.20 and −2.69, respectively. The results illustrated that the strong positive association between variables using IV for NE are 0.46, respectively, followed by SE and SW (Table 2).

The drainage network was categorized into five classes, each with a different association with the dependent variable. The analytical results of >800 m class for IV were 0.96. We illustrate the greatest rainfall class with the most positive correlation between both variables, and vice versa, in Table 2. The IV for >1000mm/year and <700 m was 1.01 − 2.47, respectively. The lithology parameter significantly influenced the GWPZ in this study, where Siwalics, Quaternary, Rawalpindi group and Paleocene occurred in the study area and have a diverse relationship with inventory data. For example, the most significant class of this factor was Quaternary rocks with an IV model values of 0.54, followed by the Rawalpindi group. We found that a < 500 m lineament buffer was the most significant class for GWPZ, strongly associating with the dependent variable. The IV model for <500 m and >5000 m class was 0.84 and −0.11, respectively. Moreover, the loamy sediment had a positive relationship with inventory data, while loamy clayey was an insignificant class (Table 2). The WOE, FR and IV for high NDVI were 0.22, 0.53 and 0.21, respectively. For low NDVI scores, they were −1.58, −0.07, and −1.37 respectively. Our research results suggested that the > 800m class was the most significant class for GWPZ map, while the <200 m class had useless association with the dependent variable (Table 2).

4.4. Validation

We generated GWPZ in this study using three bivariate models. The model’s validation was essential for authenticating GWPZ mapping (Prasad et al. 2020). The ROC curve was frequently utilized in research to statistically verify the models with high predictive values (Golkarian et al. 2018; Chen et al. 2018). We were able to validate the WOE, FR and IV models using the AUROC curve. The AUC values representing the poorest to the best forecasting models for GWPZ vary from lesser (0) to greater (1). According to the AUCROC graphical display, the accuracy of the forecast values can be categorized into five categories: 0.5–0.6 (poor), 0.6–0.7 (average), 0.7–0.8 (good), 0.8–0.9 (very good), and 0.9–1 (excellent) (Yesilnacar 2005). In this study, the SRC for WOE, FR and IV models were 0.86 (86%), 0.91 (91%) and 0.87 (87%), respectively (Figure 5a). The PRC rate for WOE, FR and IV models were 0.89 (89%), 0.93 (93%) and 0.90 (90%), respectively (Figure 5b). The validation results illustrated that all three models perform very well for the current study, but the FR model performed very well compared to WOE and IV models.

Figure 5. Models validation (a) SRC for FR, WOE and IV (b) PRC for WOE, FR and IV.

5. Discussion

The possibility for high groundwater production is often indicated by springs and wells within specific portions of study areas (Oh et al. 2011; Naghibi et al. 2017). However, many scholars have tended to apply statistical and machine learning techniques with positive outcomes to evaluate an entire region’s aquifer potential (Oh et al. 2011). In this research we used three models (WOE, FR, and WOE) to potentially improve on analyses of to produce GWPZ maps of the study area. In this study, a map with 12 parameters served as an independent variable to compute and subsequently understand its association with groundwater inventory and generate the GWPZ. A poor statistical correlation between dependent and independent variables indicates low potential zones, whereas a strong positive correlation indicates highly potential zones for groundwater. Moreover, the geographic distribution of water wells as indicated by the WOE, FR, and IV, enhances the estimation of corresponding weights of every category of precondition parameter.

Because elevations can significantly influence GWPZ mapping (Pham et al. 2021; Zhou et al. 2021; Rui et al. 2023), the fact that our final GWPZ maps contained five groups (very low, low, medium, high and very high zones) aided our understanding of the effect of this variable. For example, the most significant class for GWPZ mapping was 700–900 m elevation range, which strongly correlated with the dependent variable. Because of such a steep gradient, water cannot pool on the ground’s surface for an extended period and penetrate. Additionally, the low slope regions, i.e. <10° showed a strong association with groundwater potential while >40° showed a poor relationship with groundwater (Luo et al. 2022; Cheng et al. 2023). High-slope regions support runoff and consequently affect groundwater recharge. High drainage network in basins show that most rainfall is lost as runoff, while low drainage supports infiltration (Mandal et al. 2016). The current research observed that low drainage networks exhibited a strong association with groundwater potential and vice versa. Rainfall is the most significant parameter for GWPZ and the primary source of groundwater recharge. As a result, the water table rises in the region, mostly because it receives 1200–1400 mm/year rainfall.

In our study, the greatest association class of rainfall occurred in the highest category (>1000mm/year). Presence and flow of groundwater are facilitated by weathered and fractured strata. Movement and existence of groundwater are also influenced by geological factors, including lithology and lineament (Li et al. 2020; Maskooni et al. 2020; Qiu et al. 2023). As a fault or fracture line, a lineament is a large linear surface feature on a planet that promotes secondary porosity and permeability and supports groundwater recharge (Naghibi et al. 2018). In the present research, <500 m buffer zone of lineament showed the greatest correlation with groundwater, while >5000 m buffer was insignificant. Increased NDVI levels raise the probability of finding groundwater, and deep water, for example, is indicated by negative NDVI values (values close to 1). This study observed high associations between groundwater and a high NDVI range. While NDVI can provide valuable insights into vegetation health and its potential relationship with groundwater, we advise that is merely one variable of potentially many influences. Integrating NDVI data with other hydrological and environmental information, such as groundwater monitoring wells, precipitation data, soil properties, and land use patterns, is essential to understand groundwater dynamics and its interactions with the environment comprehensively. Urbanization or impervious surfaces increase surface runoff and decrease infiltration (U. Khan et al. 2021). In this study, urbanization, such as buildings and roads, adversely influenced the GWPZ. The GIS-based statistical models have frequently been used (Rahmati et al. 2015; Tariq, Mumtaz, et al. 2023; Tariq, Hashemi Beni, et al. 2023), and results overwhelmingly support use of the FR model for application in the same study when in similar regions.

5.1. Limitations

The precision and dependability of groundwater potential mapping can be impacted by several data and resolution-related limitations. It is essential to comprehend these constraints to evaluate the outcomes and generate plausible knowledge and recommendations. In the current study, we recognized the typical limitations were the data and resolution in GWPZ mapping i.e. the lack of groundwater data (wells) in various regions, high authentic geological data, and insufficient high advance RS data for monitoring. Moreover, our results were also constrained because of missing geophysical data. As a result, the accuracy of the findings could have been enhanced by incorporating additional thematic layers, such as aquifer thickness, groundwater depth before and after monsoon, recharge rate, and the volume of water used for domestic and agricultural uses. Unfortunately, this study could only consider a small number of thematic layers because of the restricted availability of data at the spatial scale. Furthermore, human activities that were not considered in this study but may have a substantial impact on groundwater availability and recharge rates include pumping groundwater for home or agricultural purposes.

A comprehensive strategy including enhanced data collecting, cutting-edge technologies, interdisciplinary cooperation, and continual modelling technique improvement is needed to overcome these constraints. Improving groundwater potential mapping accuracy and dependability requires frequent changes to datasets and models based on fresh data and research discoveries.

6. Conclusion

The current research attempted to evaluate the groundwater potential in the Karak District of Khyber Pakhtunkhwa, Pakistan, using GIS and RS based modelling. In the present study, groundwater has been investigated using three GIS-based bivariate to generate GWPZ map. According to the literature review and field visit of the area, 12 conditioning parameters including elevation, slope, aspect, curvature, drainage network, LULC, rainfall, lithology, lineament, NDVI, soil and road network were considered for computing their association with groundwater inventory data. These 12 parameters were integrated with inventory data in the GIS platform using WOE, FR and IV models to generate GWPZ map and classify it into five zones. The results of all applied models were validated using the AUROC technique to evaluate the reliability of the techniques. The SRC of WOE, FR and IV were 86%, 91% and 87%, while the PRC rate for WOE, FR and IV models were 89%, 93% and 90%, respectively. The current research findings can be helpful for public and commercial organizations engaged in managing groundwater resources, planning land uses, and safeguarding the environment in the studied region. Additionally, the approach of the current study can be used to examine other watersheds with similar hydro-geologic, topographic, and climatic characteristics. However, some measures can be made to improve the replenishment of groundwater, boost water supply, develop rainwater harvesting projects in urban households, etc. Additionally, advanced innovative geospatial techniques are further required to enhance the management of water resources and control overexploitation and unresearched exploration. Sustainably managed water resources can be supported by more multidisciplinary study in these areas, which can lead to more precise, reliable, and useful groundwater potential zone mapping techniques. In the future, researchers must investigate the impacts of climatic change and urbanization on groundwater recharging and potential. The current research recommends use of advanced RS data such as LIDAR, machine learning algorithms, groundwater-surface water interaction, real time monitoring, capacity building and, integrated approach in future to map GWPZ map for sustainable system. The present investigation emphasizes the need for more research that will include socioeconomic aspects, monitor over an extended period, evaluate the findings in the field, and expand the study to other regions. Future projects like these would help manage water resources sustainably, satisfy the growing demand for water in many industries, and improve the understanding of groundwater potential worldwide.

Author contributions

Abdur Rehman: methodology, software, formal analysis, visualization, data curation, writing—original draft. Fakhrul Islam: methodology, software, formal analysis, visualization, data curation, writing—original draft, investigation, validation, writing—review and editing. Ijaz Ul Islam: investigation, validation, writing review and editing. Aqil Tariq: formal analysis, visualization, data curation, writing—review and editing, Supervision. Davis J. Brain: writing—review and editing. Tehmina Bibi: writing review and editing. Waqar Ahmad: writing—original draft, investigation, validation, writing—review and editing. Liaqat Ali Waseem: writing review and editing. Shankar Karuppannan: writing review and editing. Saad Al-Ahmadi: writing review and editing, Funding. All authors have read and agreed to the published version of the manuscript.

Acknowledgements

The authors extend their appreciation to the Researchers supporting project number (RSP2024R384) King Saud University, Riyadh, Saudi Arabia. We also thank the Forest and Wildlife Research Center, and James C. Kennedy Waterfowl and Wetlands Endowed Chair, Mississippi State University, for support.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability statement

Data available on the reasonable request from the 1^st author of this article.

Additional information

Funding

The authors extend their appreciation to the Researchers supporting project number (RSP2024R384) King Saud University, Riyadh, Saudi Arabia.