Modeling Groundwater Potential Zone in a Semi ‐ Arid Region of Aseer Using Fuzzy ‐ AHP and Geoinformation Techniques

: Saudi Arabia’s arid and semi ‐ arid regions suffer from water scarcity because of climatic constraints and rapid growth of domestic and industrial water uses. The growing demand for high ‐ quality water supplies and to reduce the dependency on desalination creates an urgent need to explore groundwater resources as an alternative. The weighted overlay analysis method using the fuzzy ‐ analytical hierarchy process (FAHP) multi ‐ criteria decision making (MCDM) techniques combined with geoinformation technology was used in this study to explore the groundwater potential zones in the Itwad ‐ Khamis watershed of Saudi Arabia. Twelve thematic layers were prepared and processed in a GIS setting to produce the groundwater potential zone map (GPZM). Subsequently, potential groundwater areas were delineated and drawn into five classes: very good potential, good potential, moderate potential, poor potential, and very poor potential. The estimated GWPZ (groundwater potential zones) was validated by analyzing the existing open wells distribution and the yield data of selected wells within the studied watershed. With this quality ‐ based zoning, it was found that 82% of existing wells were located in a very good and good potential area. The statistical analysis showed that 14.6% and 28.8% of the total area were under very good and good, while 27.3% and 20.2% were accounted for the moderate and poor potential zone, respectively. To achieve sustainable groundwater management in the Aseer region, Saudi Arabia, this research provided a primary estimate and significant insights for local water managers and authorities by providing groundwater potential zone map.


Introduction
Groundwater is one of nature's most essential resources that occurs in pore spaces and fractures in rock and sediment below the surface of the earth [1]. It plays a significant role in human well-being, ecological balance, and economic growth [2]. About 30% of the world's freshwater is concealed and reserved as groundwater, whereas surface water accounts for only 0.3% in the form of lakes, marshes, reservoirs, and rivers [3]. The rainwater and snowmelt water are the primary source of groundwater, which flows in the aquifer through the soil pores. Groundwater consumption is more reliable and fresh compared to surface water due to being more convenient and less vulnerable to pollution. Since the last century, the demand for freshwater has been a growing concern worldwide due to rapid industrialization and population growth [4]. Consequently, groundwater extraction has become a significant component of water management and planning, particularly in rural regions [5]. The most vital part of groundwater resource management is to delineate the water occurrence. The groundwater potential zone delineation is crucial for sustainable resource utilization and further irrigation system development [6]. In the current water scarcity scenario, groundwater potential should be identified and evaluated through groundwater exploration. Concerning groundwater exploration, the term potential groundwater could be described as the probability of groundwater occurring in the area.
Saudi Arabia is located in one of the driest regions in the world, with average annual precipitation ranging from 80 to 140 mm [7]. Some areas at higher elevations in the west and south-western of the Kingdom, however, receive up to 500 mm of rainfall annually [8] due to topographically-driven convective events [9]. The spatial and temporal variability of rainfall is very high, along with occurrences of extreme climatic events [10]. The United Nations has already classified Saudi Arabia as water-scarce nations [11]. Future climatic predictions indicate that Saudi Arabia would experience a decrease in the water supply and rainfall, leading to water scarcity in an already water-stressed region [12]. To meet the increasing demand for water in the light of fresh surface water's limited resources, Saudi Arabia is dependent on desalination and has overlooked its high operational costs and potential effects on the environment. Therefore, research into groundwater resources is crucial in Saudi Arabia and could contribute effectively to cover the water requirements and sustainable groundwater management. However, identifying prospective groundwater sites with abundant water supply is not a simple task. Conventional approaches, such as drilling, hydrogeological, geological, field surveys, and geophysical methods, involve extensive labor in exploration activities that are costly in terms of money, time, and resources and also demand the participation of experts [13].
Groundwater occurrence in any place is the result of interactions among geology, hydrology, air, and biology factors [14]. As a consequence, lithological variation, topographical status, slope, geology, rainfall pattern, soil texture, and many more identify groundwater movement through soil pores [14][15][16]. GIS and remote sensing (RS) integration of groundwater exploration in groundwater studies is now a paradigm shift in groundwater research that helps assess, monitor, and conserve groundwater resources [17]. Satellite-based remote sensing (RS) techniques allow for greater ground surface observations than is possible through terrestrial observations [18]. It provides extensive, accurate, neutral, and readily available information about the location, scope, and dynamics of changes. It assists in the analysis process related to the effects of geographical phenomenon factors, processes of change, and magnitudes of change. The set of tools, to this extent, i.e., RS and GIS, is beneficial for the study of water resources management, seasonal environmental conditions, climate change, and ecosystem management. Remote sensing data collection facilitates the synoptic analysis of function, patterns, and changes in the Earth system over time at local, regional, and global levels. Also, this information provides an important relationship between comprehensive local water resources research and conservation and management of regional, national, and international biodiversity. Recent developments in geoinformation technology (RS and GIS) provide an important tool for groundwater mapping and exploration [18][19][20][21][22][23]. Delineating prospective areas of groundwater using RS and GIS is also a powerful tool and is very valuable in determining groundwater potential zones in arid to semi-arid climatic regions [20,[24][25][26][27][28]. However, feature classes used to delineate the groundwater potential zones (GWPZ) differ across studies, and the attribute layer selection is arbitrary. A literature review (Table 1) of expert opinion was conducted to identify key parameters affecting groundwater availability.
Integrated multi-criteria decision analysis (MCDA), geoinformation techniques were used [16,18] to identify prospective of the groundwater and artificial recharge sites. Several researchers [29][30][31][32][33] have recognized MCDA as an efficient tool for water resources management and the environment. Some studies have shown that the GIS-based multi-criteria analysis is also useful in identifying zones for groundwater recharge [34].
Based on expert opinion and the location-specific conditions, weights are assigned to different thematic layers and their feature classes. Several researchers have used the analytical hierarchy process (AHP) for the computation of the thematic parameter's relative importance [35][36][37]. AHP method for multi-criteria decision making (MCDM) has been widely used and effectively implemented in various fields of management of natural resources, environmental impact analysis, and regional planning [38][39][40][41]. AHP's graphical users (IGUs) interfaces, automatic priority, variable calculation, and sensitivity analyses have become even more significant in this area [42]. Despite its popularity, the imprecision and inherent uncertainties linked to the geoscientific visualization of a crisp number are not efficiently dealt with decision-makers [41,42]. Many researchers have also addressed the theoretical validity and empirical effectiveness of AHP [43], and this debate has concentrated on four important key fields: the axiomatic foundation, the correct meaning of priorities, the 1-9 measurement scale, and the rank reversal problem [44]. Indeed, most of the critics have been partly resolved in these areas, mainly three-level hierarchic structures [45]. The decision-making problems are hierarchically organized at different levels in the AHP method, a finite number of elements at each level; however, in many circumstances, the preferential model of the decision-maker is uncertain and fuzzy, and crisp numerical values of comparable proportions to be given by subjective perception are extremely difficult [42,45]. Because of incomplete information or knowledge and uncertainty regarding decision-making, it may be subjective and ambiguous for decision-makers about their level of preference. In the analysis of uncertainty, AHP can be combined with fuzzy logic methods [46][47][48] and setting up a framework for the assessment and reliability of the criteria through the use of fuzzy membership functions (FMFs) [49]. Assigning to each object a membership or non-membership function of each criterion, the thematic criterion in the MCDM framework is standardized using fuzzy sets [50]. The combination of an AHP with a fuzzy set theory allows for more flexibility in the analysis of results, and subsequent decision-making [49]. Table 1. Literature review of the theme used to identify groundwater potential zones (GWPZ). AHP method for MCDM has been widely used in recent decades and has been successfully applied in groundwater potential zone mapping [38][39][40][41]. However, despite the popularity of AHP, the method is sometimes criticized for being unable to adequately deal with the inherent uncertainties and imprecision associated with mapping the perception of a decision-maker to crisp numbers [42]. Further, AHP method for MCDM has been used in arid and semi-arid region for groundwater potential mapping [28,30,41], and none of the studies performed using the integration of fuzzy set theory with MCDA and, in particular, with FAHP (fuzzy-analytical hierarchy process) would lead to an improvement in the accuracy of groundwater potential maps due to the flexibility of fuzzy membership functions. Based on this context, the present paper applied the integrated approach of RS and GIS with Fuzzy-AHP to the development of thematic data layers for the delineation of GWPZ in the Itwad-Khamis watershed of Saudi Arabia. The main objectives were as follows:

Literature Review G GM SW ST DEM SL R DD LD TWI LULC VC WT
(a) propose a methodology for the identification and delineation of groundwater potential zones in Itwad-Khamis watershed using integrated techniques of RS, GIS with Fuzzy-AHP and demonstrate this by a case study; (b) carry out a sensitivity analysis and recognize the most sensitive factors that influence the identification of potential groundwater zones; (c) demonstrate geoinformation technology capabilities in groundwater mapping.
The present study, a Fuzzy-AHP MCDA method, was used to develop potential groundwater zones that could serve decision-makers, policymakers, and water resource planners to make effective and sustainable use of groundwater resources.

Study Area
The Itwad-Khamis watershed located in the Aseer region of Saudi Arabia (Figure 1) covers an area of 1208 km 2 . The study site lies between the latitude of 17°59'33.519" and 18°21'33.586" North and longitude of 42°34'4.228" and 43°6'17.364" East, and the elevation ranges between 1972 m and 2596 m in mean sea level (MSL) with the mean of 2172 m. According to the Saudi Geological Survey, the watershed area is underlain with sedimentary soft, hard silt, and clay soil. The watershed has a heterogeneous landscape in terms of terrain complexity. The climatic condition of the watershed is cold semi-arid climates, and it is influenced by high elevation. For the last 53 years (1965-2017), the average annual rainfall is 242 mm, with the majority of the rainfall occurring between Feb and May, and the average min. and max. air temperature temperatures are 9.3 °C and 30.6 °C, respectively. The region is susceptible to heavy rainfall (high intensity), and some of its villages and rural areas, within or surrounding, experience flash flooding during the winter.
The Food and Agriculture Organization (FAO) and AQUASTAT [69] reported that Saudi Arabia accounted for 88% and 9% of freshwater withdrawal by agricultural use (793 m 3 /p/year) and domestic use (81 m 3 /p/year), respectively. With the anticipated future limited desalination capacity and reuse of wastewater [70], further exploration of groundwater resources would have to meet this requirement. Therefore, socio-economic activities are related to the use of natural resources that need immediate concern for the conservation and development of water resources [71]. The prospective mapping of groundwater in the studied watershed, which is one of the most relevant regions for Eastern Afromontane biodiversity hotspots [72], will have a significant impact on the Aseer area [73]. Since the significant part of drainage contributes to Lake Tindaha and Upper Bisha Dam, the mapping groundwater potential would improve the sustainable management of the Aseer region and the kingdom. These characteristics were appropriate for this study to demonstrate the applicability of the methodology.

Data Used and Methodology
The current research demonstrated the use of different information sources for potential groundwater assessment in the Itwad-Khamis watershed. Groundwater recharge is the primary method by which water reaches an aquifer. The availability of groundwater from an aquifer at a specific location is based on the rates of withdrawal and replenishment (recharge) [74]. Areas with high recharge rates indicate a high infiltration capacity of the soil and helping to identify potentially good groundwater zones [52]. Bhattacharjee [75] estimated the annual groundwater storage based on fluctuations in water levels and specific yield approach. However, in the current study, the groundwater level data, the amount of groundwater data, and the groundwater quality data were limited due to the economic constraint and regional situation of the study area [51,53,56,60,61,65]. Therefore, in this study, an estimation of the values in regions without incorporating the recharge and groundwater head data (i.e., pumping well tests) was carried out using an MCDM model based on hydrogeological factors.
The relative influence of the individual themes on groundwater potential was determined based on the experts' opinion, knowledge of influencing factors, and practical experience. A total of eight international experts (hydrogeologists) were interviewed through a questionnaire in the field of hydrogeology and water resources engineering explicitly designed to obtain their opinions on the relative importance of hydrological/hydrogeological variables affecting groundwater occurrence (Supplementary S1). Besides, local hydrogeologists were also consulted for their views regarding the relative importance of theme and their feature classes. Twelve thematic maps for this purpose were chosen as per the expert opinion and extensive literature review (Table 1) in context with semi-arid climatic conditions, such as geology, slope, elevation, rainfall, drainage density, lineament density, land use/land cover, proximity to surface water bodies, vegetation cover, soil texture, hydraulic conductivity, and topographic wetness index (TWI). All these themes were generated using remote sensing (RS) and conventional data, with the help of ArcGIS (ESRI) 10.3 and ERDAS 9.2 imagine software (Hexagon Geospatial). The dataset was geometrically rectified to a common projection (Universal Transverse Mercator and datum WGS84) based on the topographical base map and global positioning system points, and it was resampled using the nearest-neighbor algorithm in ERDAS Imagine software. The description of the data collection process, the existing open wells map preparation, including yield data, groundwater influencing parameters, and development of spatial modeling and validation have been defined as below.

3.1.Data and Material Used
Data collection included data sources and data types used in the present study. Geological data and the geological maps from the 1:250,000 Abha quadrangle GM-75 were obtained from the geological survey of Saudi Arabia. Meteorological data and rainfall data were acquired from the PME (Presidency of Meteorology and Environment) of Saudi Arabia for the last 53 years (1965-2017) from the four (4) rain gauges stations, networked within and nearby the watershed area. Satellite images: Cloud-free Landsat 8 OLI (Operational Land Imager) Level 1T (terrain corrected) data with 30 m spatial resolution were acquired on 16 September 2018, for path 167 and row 047 with a UTM Zone 38 N coordinate system that was obtained from the archives of the United States Geological Survey website [76]. ALOS PALSAR RTC (radiometrically terrain corrected) digital elevation model (DEM) [77] with a resolution of 12.5 m was obtained from Earth Science Data Systems NASA. Field survey and the reconnaissance survey were carried out during the 16-19 September 2018; during 15 January-25 February 2019 (cloud-free date), the different land use/land cover (LULC) classes, soil samples collection, and the location of existing wells were identified from the various locations. Experimental data and the collected soil samples were analyzed in the laboratory for grain size distribution and permeability. The equipment used to accomplish this study was the global positioning system (GPS-MONTANA 680 Garmin, Garmin USA), for the positioning of the soil samples and identification of various LULC classes; for the laboratory analysis, hydrometers for the sieve analysis (grain size distribution), and falling head permeability test, BS 1377-5:1990 and ASTM D2435-04 were utilized. The data processing was done by ArcGIS 10.3 (incl. ArcHydro), ERDAS, ENVI, Microsoft office software, and SPSS software.

Spatial Data Processing for Groundwater Potential Zone Mapping
For the evaluation of the groundwater potential of the Itwad-Khamis watershed, twelve thematic maps, namely, geology, slope, rainfall, elevation, Normalized difference vegetation index (NDVI), drainage density LULC, lineament density, topographic wetness index (TWI), proximity to surface water bodies, hydraulic conductivity (K), and soil texture, were generated using remote sensing and conventional data with the help of ERDAS and ArcGIS software. Out of these thematic maps, topographic elevation, slope, drainage density, and TWI maps were generated from ALOS PALSAR RTC (radiometrically terrain corrected) DEM data, whereas the remaining maps were generated using Landsat-8 satellite and conventional data, such as geology, rainfall, surface water body, permeability (K), NDVI, LULC maps, and existing well and specific yield data. All the datasets were resampled using the nearest-neighbor algorithm to a spatial resolution of 30 m.
Geology determines the groundwater storage aquifers. The geological rock type exposed to the surface has a significant impact on groundwater recharge [78]. In the present study, the geological map was obtained from the geological survey of Saudi Arabia. The map was first scanned, rectified, and then digitized using ArcGIS software to prepare thematic layers of geology. Figure 2 shows the produced geology layer, consisting of the following classes: dioritic and gabbroic rocks, tonalite suite, Jiddah group, granodiorite and granite suite, granite suite, Jiddah and Baha groups undivided, sedimentary, volcanic, and metamorphic rocks. Table 2 shows the main rock types (group), such as biotite granodiorite and monzogranite-foliated uniform body (gdn); mixture formed of irregular layers and bodies in amphibolite (gdv) under the granodiorite and granite suite and hornblende-biotite tonalite-gneissic tonalite (gdn); hornblende diorite associated with hornblende-biotite under the tonalite suite covering a significant part of the watershed. A small coverage of sedimentary sandstone that is highly permeable and unstable is located in central-southern (Wajid Sandstone "O€w") and eastern (Qal and Tb) part of the study area. Due to the hardness and low fractures (cracks), the movement of groundwater in diorite gabbro and diorite rocks is difficult. Hence, it may result in poor groundwater potential [28]. The topographical elevation is likely to affect the occurrence of groundwater prospects and regulated by various geomorphological and hydrogeological processes (i.e., geology, meteorological conditions, land degradation, etc.) [31]. The topographical elevation map for the study area was developed by the digital elevation model (DEM) extracted from ALOS PALSAR RTC data [77]. The DEM of the study area acquired from the (ESDS) NASA website was unfilled and unprocessed, and therefore, the DEM was cleaned up by filling sinks to cover-up the local depressions. Figure 3 shows the elevation map ranges between 1972 and 2596 m. The higher elevation was found to the south and south-east, while the lower elevation was found to the north of the map. While assigning the weight to delineate the GWPZ, the high weight potential was assigned for lower elevation and low weight for higher elevation [52]. The existence and flow of groundwater are strongly controlled by the slope [59]. The slope was derived directly from the elevation layer, as shown in Figure 3. In the study area, the slope angle ranged between 0 and 67.7. The high slope was found towards central, central south, and eastern parts of the study area, while the west and northern parts consisted of a low slope. From the theoretical and conceptual perspective drawn from the comprehensive literature review, it has been shown that the slope (0-10°) with overlay rainfall extracts the area highly potential to groundwater [56,59]. The gentle slopes have been designated for groundwater in the "good" category, as they give more residence time for rainwater to percolate than the steep slopes [52]. The lineaments can be defined as the area of weakness surface that shows some linear to curvilinear features in the geological structure, such as fracture, joint, fault, etc. Visual interpretation methods [79] using Landsat-8 satellite data (band 4, 7, and 1) and geological maps were used to extract the lineaments map. Studies have shown that groundwater intensity increases with higher lineament density [51,59]. The low weight was assigned for low lineament density and high weight for high lineament density ( Figure 4). The drainage density represents the sum-up of all streams length (wadis) in a watershed divided by the total watershed area [59]. A drainage network's structural analysis helps to assess the characteristics of a groundwater recharge zone [59], and the drainage system quality relies on the geology, which offers a significant index of percolation rate. Greenbaum [80] defined drainage-length density as an essential feature for groundwater potential zone modeling. For groundwater potential zones, the high weight was attributed to high drainage density, and low weight was attributed to low drainage density. The Euclidean distance tool in the ArcGIS environment was used to prepare proximity to surface water bodies map. It has been shown that groundwater prospect is decreasing with increasing distance from the surface water bodies [55] ( Figure 4). As per a study carried out by Beven and Kirkby [81], TWI relates upslope areas as a measure of water flowing to a specific point (i.e., to the local slope), which is a measure of the subsurface lateral transmissivity. TWI was widely used to quantify topographical control of hydrological processes and represent the potential exfiltration of groundwater caused by topographical effects; as a result, higher TWI value described higher groundwater potential value [53]. Equation (1) given below was used for the calculation of TWI [81]. "α = Upslope contributing area, β = Topographic gradient (Slope)". The TWI varied from 1.89 to 25.21. The high weights were allocated to the high value of TWI and vice versa. Figure 4 shows the TWI map of the studied watershed.
In the delineation of groundwater potential zone, one of the essential weather variables was rainfall. It was, therefore, necessary to understanding the spatial and temporal phenomena of rains that were necessary for the study of water resources. Rainfall data were acquired from the PME (Presidency of Meteorology and Environment) of Saudi Arabia for the last 53 years (1965-2017) from the four (4) rain gauges stations (Figure 4). The networks were distributed within and nearby the watershed area. The average annual rainfall was calculated based on the daily rainfall data [82]. The inverse distance weighted (IDW) method [83] in the ArcGIS platform was interpolated to analyze rainfall variation. The results revealed that the rainfall ranged between 128.5 mm to 147.6 mm, with an annual mean of 138 mm (Figure 4). The spatial distribution and occurrence of rainfall depend mostly on climatic and topographic factors [84]. High weights were allocated for high rainfall and vice versa.
Globally, anthropogenic activities (directly or indirectly) are responsible for deforestation, soil erosion, groundwater depletion, loss of soil quality, causing damage to the environment [16]. LULC, due to natural surface change, is considered as one of the crucial factors for the groundwater occurrence. The maximum likelihood classification (MLC) algorithm (supervised classification) was performed to prepare the LULC map [85]. Confusion or error matrix was used to validate the classification result [86], and a good agreement was found (overall accuracy and Kappa coefficient of 89.38 and 0.8711, respectively). Figure 3 shows the LULC map of the study area. The most dominant class was the scrubland (42%) followed by exposed rock (34.7%), fallow land (8.37%), built-up area (6%), and sparse vegetation (2.61%). The high weight was allotted to the dense vegetation, sparse vegetation, agricultural cropland, and water bodies. The low weight was consigned for the built-up class, bare soil, and the exposed rocky area.
Vegetation has an important role in the stability of slopes by hydrological and mechanical processes. The high vegetation has good soil properties, for example, high organic matter and a greater accumulation, resulting in higher infiltration capacity [87]. Besides, the type and amount of vegetation in the region may also play an important role in groundwater occurrence [88]. Groundwater is an essential source of water for the preservation, abundance, and development of various ecosystems, especially in arid and semi-arid regions [88,89]. Groundwater depletion has been related to stress, losses, and a lack of a new generation of groundwater-dependent vegetation in many parts of the world [90,91]. NDVI is the conventional vegetation index that has been used for most hydro-environmental studies to derive an abundance of vegetation from remotely sensed data [92]. Aguilar et al. [93] investigated a strong linear relationship between maximum and mean NDVI and groundwater levels in dry periods; however, in wet periods, relations are much weaker. In practical terms, the algorithm isolates and normalizes the significant increase in reflection over the visible red to near-infrared wavelengths by dividing each pixel's overall brightness into those wavelengths, as shown in Equation (2). (2) where values are converted from raw digital number (DN) values to solar electromagnetic radiation reflectance in either band. This equation results in a single band data; NDVI values range from -1 to +1, and the values close to + 1 mean more vegetation cover. Figure 3 shows the highest vegetation patch in the area of the south, southwest, and wadies, as well as a good number of linear patches that were scattered over the central part. For the groundwater potential, the high value of NDVI was assigned as high weight, whereas the lower values were assigned as low weight.
Samples of the soil were obtained from various locations in the studied watershed ( Figure 5). With the use of a GPS navigator, a total of 64 soil samples were collected from the study area of approximately 1 kg of aggregate stability (0-30 cm depth). Soil sampling was performed using a stratified composite process, and the field was sub-grouped into areas with different elevations, LULC, and soil moisture. The site was then surveyed separately, and two replicates were obtained at each survey site, two to three meters apart. Every specimen was carefully weighed and sieved through a 2 mm and, subsequently, analyzed in the laboratory for its soil texture and organic matter, following the standard procedure defined in Carter [94]. The soil grain sizes (texture analysis) were determined using the hydrometer method (Stokes law). Figure 5 shows the distribution of soil types in the watershed. The five soil classes, such as sandy loam, loam, loamy sand, silt loam, and sandy clay loam, were recognized in the watershed. Sandy loam had a high rate of infiltration when wet and was, therefore, categorized as highly drained soil. Silty loam and sandy clay loam were considered as low drained soil due to low infiltration rates. Hydraulic conductivity has been studied along with other soil parameters for assessing groundwater recharge methods [95]. The groundwater stock is directly affected by surface soil characteristics. High hydraulic conductivity allows for increased infiltration and percolation through soil defining groundwater recharge zones [96]. Soil samples were collected from the field survey; after that, the hydraulic conductivity was determined using a falling head method conforming to BS 1377-5:1990 and ASTM D2435-04. Finally, the point data were interpolated using ordinary kriging (OK) method in the GIS platform to analyze the hydraulic conductivity variation in the Itwad-Khamis watershed (Figure 4). The selection of ordinary kriging (OK) was based on the study conducted by Bhunia et al. [97]. In that study, the authors used five interpolation methods, such as inverse distance weighting (IDW), local polynomial interpolation (LPI), radial basis function (RBF), ordinary kriging (OK), and empirical Bayes kriging (EBK), to generate a spatial distribution of soil properties. The results indicated that OK was a superior method with the least RMSE and highest R 2 value for interpolation of soil properties spatial distribution. The high hydraulic conductivity was located in the south-western, central, and south-eastern parts of the watershed, whereas it had a lower value towards the north-eastern region. For the groundwater intensity, high weights were considered for the high value of hydraulic conductivity and vice versa.

MCDM: Fuzzy Set Theory
Fuzzy set theory in MCDM, proposed by Zadeh [98], is a modeling technique, stimulating complicated systems that are difficult to describe by crisp numbers. Fuzzy provides a remarkably easy way to draw specific conclusions from ambiguity, vagueness, and imprecision information [99]. For spatial planning, fuzzy logic is normally included when making decisions to execute the spatial object on a map as members of a set. In classical set theory, often referred to as "crisp," an element is either part of a set or is not part of a set. Though through the fuzzy set theory, a feature object can be used as membership values between 0 and 1, this represents the degree of membership function [98]. A triangular fuzzy number (TFN) M ̃ is shown in Figure 6. TFNs are represented simply by (l/m, m/u) or (l, m, u), which is the lowest possible value, the most likely value, and the highest possible value, respectively. In the term of its membership, the TFN has a linear representation on its right and left sides that can be referred to as (Equation (3)): Below equation shows the fuzzy number of each degree of membership based on its respective left and right representation [100]: where l(y) and r(y) represent the left side and the right side of a fuzzy number, respectively.

Fuzzy Membership Function (FMF)
The primary function of the fuzzy set theory and FMF can signify ambiguous data. Math and programming functions are also allowed to use in a fuzzy environment. A fuzzy set is an object class defined through a membership function that assigns a membership value from 0 to 1 to each object and vice versa [98]. For groundwater prospective zoning mapping, a fuzzy set theory allowed the concept of partial location membership considered for multiple classes. In this conceptual background, FMFs were assigned to the analysis of spatial variance, and their pattern led to the development of fuzzy boundaries for each potential zone. FMFs were assigned to the variance, and their trend resulted in fuzzy boundaries for each potential zone being developed. The change from 0 to 1 could be determined by applying each FMF's shape.

Feature Data Standardization Using FMFs
The featured layers for the mapping of potential groundwater areas have been created in different units and also measurement levels. Four scales of estimation have been described as ordinal, nominal, interval, and ratio [101], which require a standardization of data. For this reason, the assessment process involves the incorporation of all influence factors for groundwater potential thematic layers into one output. Consequently, the fuzzy membership method has considered standardization methods. The usage of fuzzy set theory in the mapping of potential groundwater was seen as a better result [102]. Consequently, all groundwater influence factors have been in the range from low potential (0) to high potential (1). Following the objectives and hypothesis, we utilized two membership functions for groundwater potential, viz. linear FMFs. (linear increase or decrease in membership between two inputs: linearized sigmoid shape) and categorical FMFs (the expert has allocated membership value for each named class) ( Figure 6). In many fuzzy logic applications, the first two sigmoidal membership features are frequently used and provide for a progressive change from 0 (non-member) to 1 (full membership) [101], though selecting user-defined or categorical membership functions is sometimes inevitable. Figure 7 shows membership functions relied on Fuzzy-AHP, such as (Type I) linear FMFs for proximity to surface water bodies, slope, elevation, drainage density, vegetation, lineament density, rainfall, hydraulic conductivity TWI, and (Type II) categorical FMFs for geology, LULC, and soil texture. All the criteria used for the FMF and subsequent Figures 7 and 8 showed the following output raster maps.

Weights Assignments and Normalization
A total of eight international experts (hydrogeologists) were interviewed through a questionnaire in the field of hydrogeology and water resources engineering explicitly designed to obtain their opinions on the relative importance of hydrological/hydrogeological variables affecting groundwater occurrence (Supplementary S1). Besides, local hydrogeologists were also consulted for their views regarding the relative importance of theme and their feature classes. The experts' opinions and the existing literature review (Table1) were investigated for the assignment and standardization of appropriate weights. Saaty [103] suggested weight assignment, but in earlier studies, it was not considered significant. The AHP has been widely used by the MCDA to procure appropriate weights for various criteria [104]. An AHP estimates the necessary weights related to the appropriate thematic layers to support the selected matrix by comparing and analyzing all the criteria (thematic layers) identified [41,82]. The weights combined with the criterion values for each decision variant (e.g., for each location) are used to obtain a single scalar value signifying the relative strength of the variant. The objective of the AHP is to take expert knowledge into account, and as the conventional AHP cannot accurately represent the human choice based on the quantitative articulation of preferences, a fuzzy upgrade of AHP (called Fuzzy-AHP) has been developed to address the fuzzy hierarchical issues. In this research, we used the Fuzzy-AHP method to fuzzify hierarchical analysis by permitting fuzzy numbers to be used for pairwise comparisons to assess fuzzy weights. To evaluate the weights of the evaluation criteria using a FAHP, the following stages were considered [105].
Stage I: Pairwise comparison matrixes were developed using all criteria in the dimensions of the hierarchy system. The linguistic principles were applied to the pairwise assessments as follows, and in each case, one of the two parameters was more relevant.
Stage II: Buckley's geometric mean method was used to determine the criterion's fuzzy geometric mean and fuzzy weighting [107]: where is fuzzy comparison value of criterion i to criterion n; therefore, ̃ is geometric mean of fuzzy comparison value of criterion i to each criterion. In , i is the fuzzy weight of the ith criterion, and it can be shown by a TFN. = (lwi, mwi, uwi), where lwi, mwi, uwi stand for the lower, middle, and upper values of the fuzzy weight of the ith criterion, respectively.

Groundwater Potential Map Development
The GWPZ of the watershed was obtained through the integration of feature maps into the ArcGIS platform. The weighted linear combination (WLC) aggregation method was used to calculate the GWPZ [108], as shown below in Equation (7). (7) where GWPI = groundwater potential index, = FMF thematic maps, and = normalized weight of the theme. m = total number of themes, and n = total number of classes in a theme.

Sensitivity Analysis
The sensitivity test was performed to evaluate the effects on the model output performance of the input criteria and also verify the impact on the process of modifying parameters or factor conditions [109]. Data on the influence of the weighted values and weights assigned to each variable were tested by sensitivity analysis [110]. The effective weights of each criterion were, therefore, compared to the weight assigned. Equation (8) below was used to calculate the effective weight [110].
where w and s are the weight and scaled value of the theme (theme) assigned to each grid, respectively, and GWPI is the groundwater potential index, as calculated from Equation (7).

Relative Operating Characteristics (ROC)
The ROC curve analysis was used for the qualitative validation of groundwater potential datasets developed by Fuzzy-AHP-MCDA [111][112][113]. ROC curve analysis is a standard technique to assess the accuracy of a diagnostic test. The ROC curve plots, on the X-axis, the false-positive value, and, on the Y-axis, the true-positive value [114]. The area under the curve (AUC) of ROC demonstrates the accuracy of a prediction process by explaining the ability of the system to expect the accurate occurrence or non-occurrence of pre-defined "events". A curve with the largest AUC is the best method. Naghibi et al. [1] summarized the relationship between AUC and prediction accuracy as follows: excellent (0.9-1), very good (0.8-0.9), good (0.7-0.8), average (0.6-0.7), and poor (0.5-0.6). A total of 240 existing wells were identified and used to prepare the ROC curve to validate the predicted potential map of groundwater.

Results and Discussion
The relationships between the location of existing wells and twelve groundwater potential factors were identified by using GIS-based statistical models, including Fuzzy-AHP and WLC. As per the selection of criteria, expert's opinions, local hydrologist views, and literature review were considered. The theme standardization of groundwater prospect criteria was performed, and after that, the GWPZ was created based on a GIS-based groundwater potential mapping technique. The reliability of predictive models ("What is likely to be expected?") has been a significant concern in the majority of environmental modeling applications, including water resources assessments [30,31,39,49,56]. The details of the results and discussions have been discussed below.

Sensitivity Analysis
The statistical analysis is shown in Table 5 (derived from Equation (8)), confirming significant GWPI variations. For each groundwater potential parameter, the effective weights with the theoretical weight are shown in Table 5. As shown in Table 5, the significant variation of GWPI was verified by the statistical analysis (derived from Equation (8)). This was anticipated after the removal of the geological criteria, as the layer had the highest indices of variation (31.33%). It is primarily owing to the existence of the Wajid sedimentary sandstone and tertiary laterite deposited in the Palaeozoic era, where the terrain is most appropriate for groundwater (gentle slope, high porosity provides large quantities of water storage, unconsolidated). The Paleozoic Wajid sandstone generally demonstrates very high porosity (20%-31%) and permeability (500-7000 mD) [115]. A negative correlation existed between the slope and the infiltration rate. The exclusion of slope also significantly influenced the variations in the GWPI assessment (20.56). Because of the almost flat terrain and the relatively high infiltration rate, areas with a 0-6° slope fell into the 'very good' category. In the present study area, 66% of the land was covered by a 0-6° slope. The gentle slope percolated more rainwater than the steep slopes. The GWPI also seemed to be sensitive to rainfall, elevation, and lineament density (with an average variation of 10.94%, 11.20%, and 5.90%, respectively). This could be due to the high theoretical weight of all these criteria (11.02%, 10.24%, and 6.45%). The removal of the drainage density, LULC, proximity to surface water bodies, vegetation cover, TWI, hydraulic conductivity, and soil texture also led to the sensitivity value variations in the study area; their mean value was 5.48%, 4.99%, 4.29%, 3.00%, 2.69%, 2.13%, and 1.75%, respectively.

Overlay Method
A total of 240 existing wells were identified in the watershed, all of which were evaluated with the results of GWPZ. Only 13 of these 240 wells were found to exist within the very poor and poor GWPZ due to various reasons. These wells were located near dense settlements or areas of intensive agriculture. Whereas 196 wells found within the good and very good groundwater potential zones (Table 4) and remaining 31 were located within the moderate. Therefore, it was concluded that the model was valid as it indicated that 82% of wells existed in good and very good groundwater potential zones. From the study, it could be found that the GIS and Fuzzy-AHP-based techniques for delineating potential groundwater zones adopted herein are a useful method that could be applied in the watershed-based planning and development of semi-arid regions with different geo-environmental settings.

Relative Operating Characteristics (ROC)
The validation using ROC depicted that Fuzzy-AHP had a very good prediction accuracy of 0.815 with a standard error of 0.023. The results of potential groundwater maps obtained using the GIS-based Fuzzy-AHP technique were, therefore, good predictions.
The GWPZ map was also verified by the yield data collected by the Ministry of Environment, Water, and Agriculture from 13 wells. On the derived groundwater potential map, the positions of the wells were overlaid. The validation checked that the classes for higher water yields ranged from 120 to 200 m 3 /h in very good and good GWPZs. In moderate potential areas, the yield ranged from 80 to 120 m 3 /h. The very poor to poor areas of groundwater are made up of hard-rock regions, mainly consisting of basalt and andesite-pillow lava, subordinate chert, siltstone, slate and conglomerate, minor interbedded basalt, dacite, and andesite. Most of the high-yield wells lie on the sedimentary Wajid sandstone and tertiary laterite and have good aquifer recharging potential. A linear regression coefficient of 0.785 was observed between GWPI and groundwater yield data ( Figure 9).

Analysis of GWPZ Classification Map
In the Itwad-Khamis watershed (Saudi Arabia), there was a need to evaluate the groundwater potential study due to their flash-flood events and agricultural development [73]. The results of the study summarized the weighted overlay analysis approach using the GIS-based Fuzzy-AHP technique, which is one of the effective technologies for mapping potential groundwater areas. The twelve (12) thematic layers, namely, geology, slope, elevation, rainfall, NDVI, drainage density, LULC, lineament density, topographic wetness index (TWI), proximity to surface water bodies, hydraulic conductivity (K), and soil texture, classified the study area into various groundwater potential zones. Figure 10 illustrates the GWPZ map that was quantitatively developed for analysis using the groundwater potential index (GWPI). The groundwater potential index's mean value was observed to be 50 with a standard deviation of 0.097, and the min. and max. GWPI values were 0.17 and 0.80, respectively. The groundwater potential index was categorized in zones according to the histogram profile. The histogram profile showed the GWPI value, indicating the data distribution frequency. The histogram ( Figure 11) indicates that the distribution values were unevenly distributed. Therefore, the natural break classification scheme was chosen for zoning mapping. Subsequently, potential groundwater areas were identified and drawn into five classes: very good potential, good potential, moderate potential, poor potential, and very poor potential ( Figure 10). Table 6 shows the statistics of integrated groundwater categories for groundwater potential zone. According to the statistical analysis, 14.59% and 28.78% of the total area were under good and very good GWPZ, while 27.28% and 20.18% accounted for the poor and moderate potential zone, respectively. As quality-based zoning, the 82% (195 wells out of 240 wells) of the existing open wells sites were located over the very good and good potential area covered by 43% of the watershed, as indicated by the presence of sedimentary Wajid sandstone and tertiary laterite. Concerning spatial distribution and quantitative analysis, the spatial variation of good and very good potential zones of groundwater was found in the eastern, central, and southern parts of the watershed.  Furthermore, a closer analysis of the GWPZ map revealed that the distribution was more or less a consequence of the slope and rainfall pattern, in addition to the geological influence. Additionally, in the areas underlain by basalt, diorite gabbro, and diorite rocks, the movement of groundwater was difficult, especially in the south-eastern and north-western regions, characterized by relatively low annual rainfall and scarce vegetation cover. Similar observations were also noted by Rahmati et al. [28]. On the other side, due to the presence of a high density of lineament and deep weathering, areas underlain by Wajid sandstone and tertiary laterite and quartzite and quartz-schist exhibit high groundwater potential [116], while areas underlain by biotite granodiorite and monzogranite-foliated rocks have moderate groundwater potential. Moreover, high slope areas and high elevation areas could contribute to the identified poor groundwater potential in the southern part of the study area. However, high rainfall, high drainage density, and low slope that could increase water infiltration into the groundwater system could be attributed to the high groundwater potential of the western part of the study area. Vegetation plays a vital role in the stability of slopes by hydrological and mechanical processes. According to Aguilar et al. [93], maximum and mean NDVI had a good relationship with groundwater levels; in this study also, the higher NDVI value had a significant role in the existence of high groundwater potential. Das [117] studied that the wetlands, vegetation, and agricultural croplands had good groundwater potential; in the present study, it was revealed that these types of land cover had a significant role in the existence of moderate to very high groundwater potential. The methodology adopted, therefore, showed promising outcomes for predicting groundwater potential, irrespective of the opinion of the expert.
In order to take advantage of both the fuzzy technique and the AHP, Van Laarhoven and Pedrycz [118] used the principles of fuzzy logic in AHP. Over the last decade, Fuzzy-AHP was one of the widely accepted techniques for multi-criteria decision making in various applications, such as groundwater quality assessment [119], landslide [120], ecosystem and environment vulnerability [121], groundwater vulnerability analysis [122], and extracting the water quality [123]. In this study also, the competence of the proposed GIS-based Fuzzy-AHP models made the application more realistic and reliable.

Conclusions
One of the most important research areas in arid and semi-arid climatic regions is the study of groundwater resources. In the Itwad-Khamis watershed (situated in South-West Saudi Arabia), there was a need to evaluate the groundwater potential study due to their flash-flood events and agricultural development [73]. In this study, the weighted overlay analysis approach using the GIS-based Fuzzy-AHP technique was utilized to explore groundwater potential zones in the Itwad-Khamis watershed of Saudi Arabia. The advantage of such a mapping approach lied in its ability to use fewer validation data (e.g., such as locations of well and yield data) and thematic layers of hydro-geomorphology, still maintaining optimal mapping accuracy. A total of twelve (12) thematic layers (viz. geology, slope, elevation, rainfall, NDVI, drainage density, LULC, lineament density, topographic wetness index (TWI), proximity to surface water bodies, hydraulic conductivity (K), and soil texture) were identified and prepared as thematic layers to be integrated and analyzed in the GIS setting. Subsequently, potential groundwater areas were identified and drawn into five classes: very good potential, good potential, moderate potential, poor potential, and very poor potential. The estimated GWPZ was verified by examining the distribution of existing open wells and the yield data of selected wells. With this quality-based zoning, the 82% (195 wells out of 240 wells) existing open wells sites were located over the very good and good potential area covered by 43% of the watershed. The statistical analysis showed that 14.59% and 28.78% of the total area were under very good and good, while 27.28% and 20.18% were accounted for the moderate and poor potential zone, respectively. The efforts of geo-hydrological investigation identified some interesting points about the studied watershed. Since the significant part of drainage contributes to Lake Tindaha and Upper Bisha Dam, the mapping groundwater potential would improve the sustainable management of the Aseer region. The sensitivity analysis was performed to understand the role of each parameter in detail. The results showed that geology and slope had a significant influence on the variations in the GWPI assessment (31.33% and 20.56%, respectively). Because of the almost flat terrain and the relatively high infiltration rate, the areas with a 0-6° slope fell into the 'very good' category. In the present study area, 66% of the land was covered by a 0-6° slope. The gentle slope percolated more rainwater than the steep slopes. Alternatively, the least important considerations were TWI and soil texture.
The approach presented in this paper should be further verified by conducting groundwater well discharge, groundwater depth level data, and also step-drawdown pumping well tests at various locations of the watershed to evaluate the specific yield of groundwater wells in various GWPZ. In this way, sustainable groundwater from the site can be extracted. Nonetheless, in the absence of recharge and groundwater head data and rigorous validation, the proposed methodology might be used as a primary estimate of the groundwater prospect, and the identified potential zones for the drilling of water wells should be preferred.
The present study made a significant contribution to understanding the potential of groundwater areas in the watershed studied that can be used as the initial framework by engineers, hydrologists, decision-makers, and regional planners to the replenishment of this precious life-sustaining resource. Since this study's approach was based on logical conditions and generic characteristics, the same would apply to other semi-arid and arid regions of the world with/without necessary modifications. With respect to developing effective policies and systems for the sound use and exploration of groundwater resources, a cohesive strategy is also required, in particular, by ministries, government organizations, NGOs, and people. Future studies would aim to explore groundwater potential by other advanced artificial intelligence methods, such as fuzzy-topics, artificial neural networks, etc. Furthermore, groundwater quality should be integrated as a criterion for decision-making to better understand potential drinking water.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Table S1.  Acknowledgments: Authors thankfully acknowledge the Deanship of Scientific Research for providing administrative and financial supports.

Conflicts of Interest:
The authors declare no conflict of interest.