Mapping Potential Water Resource Areas Using GIS-Based Frequency Ratio and Evidential Belief Function

Abstract

Groundwater is a critical freshwater resource that is necessary for sustaining life. Thus, targeting prospective groundwater zones is crucial for the extraction, use, and management of water resources. In this study, we combined the remote sensing, GIS-based frequency ratio (FR), and evidential belief function (EBF) techniques into a model to delineate and quantify prospective groundwater zones. To accomplish this, we processed Shuttle Radar Topography Mission (SRTM), Landsat-8 Operational Land Imager (OLI), Sentinel-2, and rainfall data to reveal the geomorphic, hydrologic, and structural elements and climatic conditions of the study area, which is downstream of the Yellow River basin, China. We processed, quantified, and combined twelve factors (the elevation, slope, aspect, drainage density, lineament density, distance to rivers, NDVI, TWI, SPI, TRI, land use/cover, and rainfall intensity) that control the groundwater infiltration and occurrence using the GIS-based FR and EBF models to produce groundwater potential zones (GWPZs). We used the natural breaks classifier to categorize the groundwater likelihood at each location as very low, low, moderate, high, or very high. The FR model exhibited a better performance than the EBF model, as evidenced by the area under the curve (AUC) assessment of the groundwater potential predictions (FR AUCs of 0.707 and 0.734, and EBF AUCs of 0.665 and 0.690). Combining the FR and EBF models into the FR–EBF model increased the accuracy (AUC = 0.716 and 0.747), and it increased the areas of very high and moderate potentiality to 1.97% of the entire area, instead of the 0.39 and 0.78% of the FR and EBF models, respectively. The integration of remote sensing and GIS-data-driven techniques is crucial for the mapping of groundwater prospective zones.

1. Introduction

Water scarcity is a global issue because of the interdependence of water, food, and energy, as well as the impact that these elements have on human livelihoods and national and international economies [1,2]. Water has always been the most crucial natural resource for maintaining life on Earth [3]. Groundwater is one of the most precious natural resources in every climate zone on the globe and is a prominent water source [4,5]. Moreover, groundwater is essential for socioeconomic development, even though its quality and quantity widely vary, because it meets a variety of human demands, including those for drinking water, irrigation, forestry, industrial uses, and animal support [6]. Because groundwater is less susceptible to environmental degradation than surface water, we can use it in a more sanitary and reliable manner [7,8], and particularly in areas with polluted water surfaces. The infiltration of river water improves the riparian groundwater capacity for denitrification. Therefore, the riverbank infiltration process is widely accepted. For the more than 2.5 billion people who live on the planet, these fresh water supplies are a daily requirement [9], and the demand for groundwater will increase due to global population expansion and its effects on urbanization, industrialization, and agricultural operations, which may be the main contributing factors to the world’s groundwater scarcity [10]. Moreover, the increased need for groundwater and ongoing climatic change have inspired the research community to develop quantitative methods to predict groundwater availability [11,12].

Groundwater is a type of water resource that fills joints, voids, and pore spaces in the soil found inside geologic formations and layers. The presence and infiltration of groundwater in rocks depends on the hydraulic conductivity of the lithologic materials, which is a result of the materials’ porosity, permeability, and fluid flow through geologic structures [2,13,14]. Researchers describe the use of groundwater mapping as a tool for water resource development and planning in [15,16,17]. The topography, geological structure, landform, drainage pattern, land use/cover, and climate are some of the geoenvironmental factors that affect groundwater availability [18,19,20]. The zones with considerable groundwater capacities that can be predicted and exploited are known as “groundwater potential zones” (GWPZs) [21,22,23].

Because groundwater resources are buried beneath strata, we need to use predictive models to investigate and uncover them [17,24,25]. The groundwater discovery techniques that are currently used are time-consuming and costly. For groundwater exploration and prediction, as well as for regional estimations, we can use remote sensing (RS) and geographic information systems (GISs) [26,27]. We can use the GIS-based knowledge and/or data-driven techniques to combine and analyze large volumes of geographical data to predict and discover new water sources [28,29]. Researchers have demonstrated the value of using RS and GISs to identify possible groundwater resource locations in several studies [15,25,30,31]. For instance, multi-criteria decision making is a fast and cost-effective method [17,25,29,32].

Researchers have effectively accomplished the exploration of GWPZ maps in various environmental contexts using the knowledge-driven and data-driven methodologies. They always employ GIS technology, which can process vast amounts of spatial data and combine different data types to predict and locate more water supplies. The knowledge-driven techniques include overlay analyses [16], the analytical hierarchy process (AHP) [27,29,33], Boolean logic [34], index overlays, and fuzzy methods [35,36,37]. The data-driven techniques include, among other, the frequency ratio [38,39,40,41,42,43], machine learning models [28], logistic regression [38,44], random forests [41,45], the weight-of-evidence method [45,46,47], the evidential belief function [48,49,50,51,52,53,54,55,56], support vector machines [57,58], artificial neural networks [46,58], and maximum entropy. We can use these techniques to successfully model groundwater availability and flood hazards [55,59,60,61].

Researchers frequently use data-driven models and techniques (e.g., frequency ratio, evidential belief function, WOE, machine learning models, logistic regression, weights of evidence, linear regression decision tree analyses, and neural networks) [55,62]. We used the frequency ratio (FR) technique, which is a data-driven model and bivariate statistical method, to calculate the spatial relationships between the independent variables and classes of thematic layers and the dependent variable, which was the existing wells (spring wells, observation wells), to assign rating (r) values to each class [6,63,64]. In many diverse situations, researchers commonly use the FR to map the groundwater potential [6,50,62,63,65,66,67,68]. For predicting groundwater occurrences, they use several criteria in this technique, including geological, topographic, climatic, and hydrologic data. Several authors have applied the EBF to GWPZs [28,50,54,69].

In the present study, our aim was to apply the FR and EBF data-driven techniques using GIS-based approaches for the delineation and identification of prospective groundwater resource areas through multi-criteria analyses derived from remote sensing data. This technique was useful for groundwater extraction, management, and prediction, as well as for the reliability assessments of the combined FR–EBF model compared with the individual models.

2. Study Area

The study area is situated in Shandong Province’s lower sections of the Yellow River, and it borders the Bohai Sea to the east (Figure 1), covering about 7438 km². The muddy waters of the Yellow River give the area its name. In 1958, 2.1 billion tons of silt were dumped into the ocean, and every year since, around 30 km² of new land forms. This region experiences a warm temperate semi-humid monsoon climate. The annual average rainfall is between 530 and 630 mm, which is considerably less than the yearly average evaporation of between 1470 and 2246 mm [70,71]. In western Shandong, the Yellow River flows along a levee that is higher than the surrounding landscape. After traversing various regions of the province, it drains into the sea around Shandong’s northern shore. The topography of the area is characterized by steep hills in the southwest and level terrain in the northeast. The alluvial marine plains in the coastal locations have a ground elevation of only 2–5 m, compared with the northern Shandong hills, with elevations around 10–30 m [48]. The average temperature ranges from 5 to 1 °C (from 23 to 34 °F) in January and from 24 to 28 °C (from 75 to 82 °F) in July. The two main methods for shallow groundwater recharging are precipitation and river infiltration. The Quaternary sediment thickness in the current study area is about 350 m, and the Holocene marine and deltaic sediment thickness is about 26 m [71].

3. Data Used and Methods

In this study, we integrated remote sensing data from a variety of sensors with topographic, hydrologic, and meteorological data to depict the possible water resource locations. We used GIS techniques to merge several thematic maps (elevation, slope, aspect, terrain roughness index, SPI, lineaments, drainage density, distance to river, TWI, NDVI, LU/LC, and rainfall intensity) produced from these data. These twelve factors are important contributors to groundwater occurrence that represent the climatic, hydrologic, land-cover, and topographic factors. We describe these factors in Section 4. We present the identification approaches of the FR and EBF models used to map the groundwater potential zones in Figure 2.

We initiated the digital elevation models (DEMs) using SRTM (Shuttle Radar Topography Mission) data (90 m cell size). We used the Shuttle Radar Topographic Mission SRTM 30 m resolution NASADEM 1 arc second WGS84 data to examine the topographic changes (NASADEM 1 arc second WGS84). The DEMs allowed for the extraction of the slope, aspect, and TRI. We produced the stream network using the 8D approach [72], and we used a GIS application to produce the stream density map, TWI, SPI, and distance to rivers. We calculated the distances to rivers using the Euclidean distance spatial tool in ArcGIS.

Landsat 8, which was deployed on 11 February 2013, carries the operational land imager (OLI) and thermal infrared sensor (TIRS). In the mapping of the land cover and use, the visible-near-infrared and shortwave-infrared wavelength ranges are used. A Landsat-8 OLI scene (path 121/row 034) that was acquired on 30 December 2021 underwent image alterations and enhancement procedures. We stacked bands 2, 3, 4, 5, and 7 to a 30 m spatial resolution and projected them to WGS-84, UTM Zone 50 N. We used the data from Landsat 8 to reveal the vegetation and water resource signatures. We computed the normalized difference vegetation index (NDVI) to reveal the vegetation areas [25], and we estimated them by applying the visible infrared bands (NDVI = NIR (band 5) − R (band 4) and NIR (band 5) + R (band 4)). In addition to Landsat OLI, we used a Landsat-5 scene (path 121/row 034) that was acquired on 24 August 1992 (bands 7, 4, and 2 in R, G, and B, respectively) in the change detection in the present study (

Table 1. Data used in present study.

No	Type of Data	Source	Date	Resolution
1	Landsat-8 OLI	USGS	30 December 2021	bands 2, 3, 4, 5, 6, and 7 (30 m)
2	Sentinel-1	ESA/Copernicus	18 October 2022	bands 2, 3, 4, 8, (“10” m, 11, and 12 (“20” m)
3	SRTM DEM	USGS	11–22 February 2000	C-band (30 m)
4	Climatic Research Unit data	crudata.uea.ac.uk	1 January 2011–1 January 2020	0.5 degrees in latitude and longitudes

).

The Sentinel-2A satellite was launched on 23 June 2015, and the first data were collected a few days later. The Sentinel-2 sensors collect the VNIR, SWIR, and TIR data. The spatial resolution of these bands is 10–60 m. Sentinel-2 records a total of 13 bands of the VNIR and SWIR spectra. The coastal B1 (443 nm) band has a 60 m pixel geometry; however, the VNIR blue B2 (490 nm), green B3 (560 nm), red B4 (665 nm), and infrared B8 (842 nm) bands have 10 m pixels. The SWIR bands (B11: 1610 nm; B12: 2190 nm) have 20 m wide pixels. Sentinel-2 satellites have temporal resolutions of 10 and 5 days, which makes them extremely valuable for future investigations. Sentinel-2 scenes are delivered as zip-compressed files in Sentinel’s own SAFE format. The spectral bands are stored as jpg files in this SAFE file at three different geometric resolutions (10, 20, and 60 m). The jpg files of bands B2, B3, B4, and B8, with spatial resolutions of 10 m, and of bands B11 and B12, with 20 m resolutions, are stacked into a single GeoTIFF file of a uniform pixel size (10 m). We obtained a subset of these data during preprocessing using SNAP software to minimize the computational time and data. We extracted the LC/LU differences using images from two dates (3 November 2018 and 18 October 2022).

The Climatic Research Unit obtained the data on the average rainfall with a 0.5° resolution. The period covered by the collected average rainfall data spanned from January 2011 to 2020, and we interpolated them using the kriging tool of Arc GIS (https://crudata.uea.ac.uk/cru/data/hrg/ accessed on 01 November 2022).

In this study, we applied the data-driven application of the EBF and FR methodologies to obtain more accurate estimates of the potential groundwater areas. We compared the EBF results with the results of well-known statistical methods, including logistic regression (LR) and the frequency ratio (FR). The EBF model produces better results in terms of groundwater [50] and flood hazards [73], and it produces better results over fuzzy logic models in mapping landslides [74]. Moreover, it produces better results in revealing land subsidence [75]. We can use the model to analyze the effects of all classes of each factor and the correlation between each component, which gives it an advantage over other statistical methods [73]. Additionally, the frequency ratio (FR) is one of the better-known statistical methods [76], and it is more efficient and dependable than the knowledge-driven AHP model [77].

We produced the groundwater potential map using a data-driven FR–EBF model built on a geographic information system (GIS), and we set or resampled all the input data to a 30 m resolution following the spatial resolution of digital elevation data. In numerous prediction techniques, researchers frequently employ this kind of multicriterion decision-making process. The user chooses the relative importance of each observation in the model, which is based on remote sensing, hydrologic, and geologic data. In the GIS approach, a raster combination is used in which each layer’s pixel is associated with a particular geographic location. As a result, the combination process is more suited to the integration of traits from several datasets into an output layer. Moreover, for each pixel in the study area, we can calculate the GWPZ indices using the EBF model. We calculated the likelihood of a well occurrence in each unit cell to infer the spatial relationship between the wells and each thematic evidentiary layer. To calculate the data-driven EBF application, we used the number of well cells (with at least one cell) in each evidentiary layer.

The frequency ratio (FR) model, which is used to represent the incidence probability for a specific attribute, is a straightforward geospatial assessment tool that we can use to ascertain the final relationship between the GRW potential and its effective factors [62,68]. To determine the frequency ratio and to apply it to the overall recharge, we determined the recharge occurrence ratio for each subclass of conditioning factors. We calculated the surface ratio for each class and measured it against the total watershed area in Step 9. We calculated the FR values for each subclass of GWR potential effective factors based on how well they correlated with the GWR potential inventory.

4. Results

4.1. Elevation

The elevation is a topographical feature that is used as a surface indication to investigate the groundwater potential. Elevation change can affect the climate, which can change the rainfall, soil quality, vegetation, land uses, and vegetation types [62]. The elevation is an important layer that governs the direction of the water flow over the land, and it controls the groundwater occurrence and recharge potential [78]. We categorized the elevation map of the study area (Figure 3a) into five zones: (1) 0–3; (2) 3–7; (3) 7–12; (4) 12–19; (5) 19–75, which covered 29.35, 48.04, 19.55, 2.86, and 0.20% of the entire basin, respectively.

4.2. Slope

The slope is one of the indications of prospective groundwater occurrences, and the topography and/or slope gradient have a direct impact on how much rainwater is infiltrated [79]. Additionally, the slope can provide a general indication of the groundwater movement direction [80]. While a moderate-to-steep slope increases the runoff water, a low or almost level slope has a high infiltration rate and little runoff, which leads to good groundwater recharging [81]. We can determine the slope via either digital elevation models from physical measurements, such as STRM DEMs [82,83,84,85], or from survey base points and topographical contours, both of which allow for the creation of a DEM [63]. We classified the slope angles of the study area, which ranged from 0 to 44 degrees, into four classes: (1) 0–1; (2) 1–2; (3) 2–3.6; (4) 3.6–44 (Figure 3b), which covered 40.87, 36.36, 19.26, and 3.51% of the entire area, respectively (Table 1). Areas with low slope angles are groundwater potential areas with high potentiality; however, steep slopes promote runoff [86].

4.3. Aspect

The infiltration rate is substantially influenced by the slope aspect, which also affects the GWR, which is a measure of the solar radiation. The slope aspect, which ranges from 0 to 360 degrees in a clockwise rotation, is the slope’s orientation. We created the slope factor using the DEM in ArcGIS 10.6, and we used the DEM in ArcGIS 10.8 to calculate the aspect. The aspect measurements included ten categories: (1) flat (or no aspect direction); (2) north; (3) northeast; (4) east; (5) southeast; (6) south; (7) southwest; (8) west; (9) northwest; (10) north. The northern hemisphere has more water resources on the slopes that face north and east than on the slopes that face south and west. The slopes of the mountains to the east and north receive less sunshine than the slopes to the south and west [87]. Despite the high soil moisture levels on the slopes that face north and east, the transpiration is modest, and vegetation has increased on the northern and eastern faces as a result. An increase in vegetation can improve the surface infiltration and groundwater recharge in some areas [26,87]. We divided the aspect map (Figure 3c) into ten classes: (1) flat (−1); (2) N (0–22.5); (3) NE (22.5–67.5); (4) E (67.5–112.5); (5) SE (112.5–157.5); (6) S (157.5–202.5); (7) SW (202.5–247.5); (8) W (247.5–292.5); (9) NW (292.5–337.5); (10) north (337.5–360), which covered 8.77, 9.68, 8.74, 8.89, 10.52, 8.24, 9.91, 9.49, 9.91, and 15.84% of the entire area, respectively.

4.4. TRI

The topography roughness index (TRI), which is a geomorphometric statistic that we used to describe and quantify the spatial distributions of the hills and valleys in the research area, also has an impact on the groundwater potential. The TRI was created to assess landscape heterogeneities, and we can use it to identify groundwater [4,88]:

$$TRI=\sqrt{(ma{x}^2}-mi{n}^2)$$

(1)

We classified the TRI result map of the study area, which ranged from 0 to 0.89, into five classes: (1) 0.00–0.11; (2) 0.11–0.31; (3) 0.32–0.47; (4) 0.48–0.62; (5) 0.63–0.89 (Figure 3d), which covered 6.12, 8.84, 31.30, 38.06, and 15.68% of the entire area, respectively (Table 1).

4.5. Drainage Density (Dd)

Drainage is a crucial hydrogeological control mechanism. Surface and subsurface structures are reflections of drainage patterns [89,90], while the relationship between the permeability and drainage density is inverse. The type of vegetation, the soil’s ability to absorb rainwater, infiltration, and the slope gradient are all factors that influence an area’s drainage system. The bedrock type and structure also have impacts [91]. When there is a high drainage density in a low-permeable-surface region, there is also high precipitation runoff from that area [89,92]. A high-drainage-density zone contributes to a high amount of surface runoff with a low groundwater recharge volume [93,94], and the overall length of the drainage densities correlates with the groundwater recharge volume [95]. Areas with less drainage density have greater infiltration and less surface runoff, which means that groundwater development is appropriate in places with low drainage densities [96,97]. Furthermore, because it is used to measure the surface runoff, the drainage density indirectly reveals the groundwater recharge [98]. Low infiltration and higher surface runoff due to a high Dd result in lower groundwater potential [63,99]. We classified the Dd result map of the study area, which ranged from 5.82 to 244.6, into five classes: (1) 5.82–53.57; (2) 53.58–101.3; (3) 101.3–149.1; (4) 149.1–196.8; (5) 196.9–244.6 (Figure 4a,b), which covered 12.26, 33.02, 32.53, 18.48, and 3.71% of the entire area, respectively (Table 1).

4.6. Topographic Wetness Index (TWI)

We used the TWI map as an indicator of the slope, elevation, and landform effects on the groundwater resources [100]. In the TWI, the upslope area is an indicator of flowing water (local slope) or a gauge of the subsurface lateral transmissivity [45,50]. The regional diversity of the hydrological conditions, including the soil moisture, substantially depends on the TWI [101]. The extent and zoning of saturated areas are influenced by the TWI, which, in turn, affects the incidence of springs [44,45]; thus, the higher the TWI, the greater the groundwater potential [7]. The TWI calculations include a summary of the topographic roughness, hillslope, and foothill effects on the lateral groundwater flow [102]. Areas with high TWIs allow us to locate areas with infiltration potential and soil moisture accumulation, which are unusual in foothills [103,104].

$$TWI=\ln (\frac{Ac}{tans})$$

(2)

We classified the TWI result map of the study area into three classes: (1) from −8.16 to −2.97; (2) from −2.97 to 0.25; (3) from 0.25 to 13.51 (Figure 4c), which covered areas of 56.19, 31.89, and 11.92% of the entire area, respectively (Table 1).

4.7. SPI

The SPI is an important index for mapping groundwater infiltration zones. The SPI is an indicator of the slope, elevation, and landform effects on the groundwater resources [62]. Using slope and flow accumulation parameters, we computed the SPI using Arc GIS software version 10.5. The runoff influence increased as the SPI value increased [105,106]. We classified the SPI result map of the study area, which ranged from 0 to 34.002, into three classes: (1) 0–0.001; (2) 0.001–0.1; (3) 0.1–34.002 (Figure 4d), which covered areas of 94.67, 4.79, and 0.54% of the entire area, respectively (Table 1).

4.8. Distance to River

In hydrogeological research, the distance from hydrographic networks is crucial because the local alluvial layers are typically found close to river courses, and especially in semiarid environments [107]. Riverfront locations are ideal for efficient infiltration, and thus, for groundwater recharge [44,55]. The groundwater potential is stronger in areas near rivers, and particularly those with continuous or extended flows [43,59,108]. However, in places farther than 600 m, it is challenging to locate the alluvial layers [4]. Rivers contribute to the GWR, and hence, affect its potential in watersheds. To initiate the distance categories, we utilized the Euclidean distance tool from the spatial analyst tools in Arc GIS 10.6 [62]. We classified the map of the distance to river results of the study area, which ranged from 0 to 0.019, into five classes: (1) 0.00–0.0039; (2) 0.004–0.0078; (3) 0.0079–0.012; (4) 0.013–0.016; (5) 0.017–0.019 (Figure 5a), which covered 28.96, 34.66, 21.77, 31.12, and 1.48% of the entire area, respectively (Table 1).

4.9. Lineaments

Linear earth features that result from geologic structures, which are known as lineaments, serve as overall surface depictions of the subsurface cracks [109]. Geological features with linear or curved shapes play a major role in the occurrence and transport of groundwater over crystalline terrain. Lineaments, such as cracks, fissures, and joints, cause the infiltration of surface runoff and the replenishment of hard-rock aquifers, and they often form because of tectonic stress/strain. Numerous academics have highlighted the connection between lineaments and the presence of groundwater to emphasize how the lineament density affects well yields [110,111,112]. Lineaments fall under the category of secondary porosity, and we can distinguish them from other terrain characteristics on satellite imagery by tonal differences. A lineament can be a fault, fracture, or master joint, or a lengthy linear geological structure, topographic linearity, or a stream’s straight course [113]. Lineaments have a substantial impact on the groundwater storage and flow and surface runoff infiltration into the subsurface [114]. We classified the lineament density map of the study area (Figure 5b,c) into five zones: (1) 0–12; (2) 12–24; (3) 24–35; (4) 35–47; (5) 48–59, which covered 34.09, 22.17, 23.50, 14.40, and 5.84% of the entire basin, respectively. We gave dense lineament zones that were deemed potential groundwater recharge zones and that promoted well production to a higher rank (“numeric value”).

4.10. Land Use/Cover

Land use/cover includes a region’s vegetation cover, aquatic bodies, residential neighborhood distribution, urbanization, and soil deposit types [115]. Land use/cover has a substantial impact on the groundwater recharge, occurrence, and availability [116,117]. For instance, agricultural and plantation lands are good locations for groundwater investigations because trees and plants can collect water on their foliage and allow it to seep into the earth through their roots and biological processes, which replenishes the water table [118]. However, because of the usage of concrete floors, which only allow for surface drainage, settlements are poor groundwater areas [85]. Typically, researchers use topographic maps and a combination of field verification and remotely sensed data to map land-use and land-cover patterns. False color composites (FCCs) are typically created from satellite images utilizing a variety of band combinations [84,119].

The recharge in vegetated areas is substantially lower than that in non-vegetated areas, according to Gee et al. [120]. Furthermore, compared with permanent lands, such as shrub and forest areas, recharging is higher in agricultural and grassland areas [121]. To identify the different types of LUs, we used Landsat-8 data. Additionally, using a supervised classification strategy and the maximum likelihood approach, we classified the LU/LC result map of the study area into three classes: (1) water, (2) vegetation, and (3) urban (Figure 5d), which covered 13.19, 20.24, and 46.57% of the entire area, respectively (Table 1).

4.11. NDVI

Researchers frequently use the NDVI for groundwater potential zones [122,123]. Using a map of the normalized difference vegetation index (NDVI), we depicted the vegetation density and coverage. We used Envi 5.3 to prepare the NDVI layer. The NDVI has a value range from −1 to 1 (often known as from 0 to 255). A greater NDVI value indicates thick vegetation. We classified the NDVI result map of the study area, which ranged from 0 to 255, into five classes based on the natural breaks method: (1) 0.00–66; (2) 66–131; (3) 131–174; (4) 174–210; (5) 210–255 (Figure 5e), which covered 8.52, 8.61, 28.13, 39.45, and 15.28% of the entire area, respectively (Table 1).

4.12. Rainfall

The rainfall amount is one of the key elements in the identification of groundwater potential zones. Rainfall is a hydrologic technique for refilling aquifers [97]. The groundwater likelihood in a given land area increases with the rainfall [124]. Researchers can track, observe, and measure the analyzed watershed’s downpour rates using rainfall data from the TRMM satellite. Due to the substantial precipitation and subsequent flooding from storms, the infrastructure in the area has been severely damaged, and there have been losses in the provision of water resources. The identification of areas that are prone to water accumulation may be possible using rainfall information because precipitation is crucial for recharging groundwater [125]. We classified the mean annual rainfall intensity (2011–2020) map of the study area (Figure 5f) into five zones: (1) 567.4–575.8; (2) 575.9–584.2; (3) 584.3–592.5; (4) 592.6–600.9; (5) 601–609.3, which covered 32.05, 38.89, 11.33, 9.31, and 8.42% of the entire basin, respectively. We gave the dense rainfall zones that were deemed potential groundwater recharge zones and that promoted well production to a higher rank (“numeric value”).

5. Results

5.1. Application of FR Model

Researchers use the FR to describe the likelihood that a phenomenon possesses a particular characteristic [28]. We entered the calibration well locations (the dependent variable) and conditioning factors (such as the lithology, slope, distance to streams, and distance to faults) using the FR approach [126]. We calculated the well occurrence chance in each class for all the parameters using the FR [28]. We used the following equation to determine each class’s impact on the independent variable:

$$\mathrm{Fr}=\frac{\left(P_t/T{p}_t\right)}{\left(N_c/T_{Nc}\right)}$$

(3)

where FR is the rate at which each class of each parameter impacts the parameter; P_t denotes the proportion of points in each class; T_pt is the total number of points across all classes; N_C is the number of classes; T_NC is the total number of pixels. The FR value produced for each class of conditioning factors determines the weight of each class of conditioning factors in the subject layers for mapping the groundwater potential. After applying the FR to each class, the new thematic layer is the input value for the hybrid model.

We identified the spatial correlations between the well locations (a proxy for the groundwater potential) and conditioning factors using the FR model (Table 1). According to Lee and Pradhan [127], FR values below 1 indicate low correlations, and those above 1 indicate strong correlations. The elevation classes from 19 to 75 had the highest FRs (9.06), followed by the elevation classes from 12 to 19 (3.81), from 7 to 12 (1.48), from 3 to 7 (0.87), and from 0 to 3 (0.55). Because this range indicates a flat plateau, the areas with high elevations (between 19 and 75) are more likely to retain groundwater. The zones with slopes from 0 to 1 had high FRs (1.157), as did the E, NE, Flat, and SE zones (2.65, 1.45, 1.24, and 1.03, respectively), as the area’s general slope is due east. Additionally, the areas with higher TRIs (from 0.63 to 0.89) and low terrains (from 0.11 to 0.31) had FRs of 1.27 and 1.23, respectively.

In terms of the hydrologic characteristics, the zones with low Dd values (from 53.58 to 101.3) had higher FRs (1.26); however, the areas of high Dd values had the lowest FRs (Table 1). The groundwater potential probability was the highest for the TWI class from −2.97 to 0.25 (FR = 1.026), followed by the TWI classes from −8.16 to 2.97 (FR = 1.003) and from 0.25 to 13.51 (FR = 0.915). The SPI classes from 0.0 to 0.001 (FR = 1.018) and from 0.001 to 0.1 (FR = 0.760) had high probabilities; however, the areas of high SPIs (from 0.1 to 34.002) had no potentiality (FR = 0). Furthermore, the areas nearest to the rivers from 0 to 0.0039 had the highest FRs (1.130), followed by those from 0.004 to 0.0078 (0.997); however, the areas farther from the rivers had no potentiality (FR = 0).

The lineament density values from 36 to 47 (FR = 1.136) were the highest, followed by those from 0 to 12 (FR = 1.120); however, the moderate LD areas (from 25 to 35) had the lowest FR values (FR = 0.851). Using the LU/LC, the groundwater potential probability was the highest for the vegetation class (FR = 1.08), as well as for areas with high NDVIs (from 210 to 255) (FR = 1.42); however, the areas from 0 to 66 had low GWPs (FR = 0.427). In terms of the rainfall classes, the areas that received from 114.42 to 148.13 mm had the highest precipitation (FR = 1.56), followed by those that received from 148.13 to 244.56 mm (FR = 1.512) and from 570 to 580 mm (FR = 1.36); however, the areas that received from 580 to 590 mm had the lowest FRs (0.608).

5.2. Application of EBF Model

For the spatial prediction of possible groundwater areas, we used the Dempster–Shafer-based evidential belief function (EBF) [128,129,130]. This strategy is a potentially effective way to address some of the challenges that are associated with the combination of evidence and data fusion [131]. We can calculate the EBF as follows [132]:

$$\mathrm{BeI}{C}_{ij}=\frac{1W_{C_{ij}D}}{{{\displaystyle \sum }}_j^m=1W_{C_{ij}D}}$$

(4)

$$W_{C_{ij}D}=\frac{N\left(C_{ij}{{\displaystyle \cap }}^{\hspace{1em}}D\right)/N\left(C_{ij}\right)}{N\left(D\right)-N\left(C_{ij}{{\displaystyle \cap }}^{\hspace{1em}}D\right)/N\left(T\right)-N\left(C_{ij}\right)]}$$

(5)

The numerator used to calculate the $W_{C_{ij}D}$ is the likelihood that the D (groundwater) occurs in the presence of C_ij (a particular class of one conditioning factor), which demonstrates that the D more frequently occurs (or is more present) in the C_ij than would be predicted by chance. The denominator used to calculate the $W_{C_{ij}D}$ is the conditional probability that D exists given the lack of C_ij, which indicates that D more frequently occurs outside of C_ij than would be predicted by chance. Therefore, the parameter $W_{C_{ij}D}$ is the weight of the C_ij in terms of D being more present than absent as anticipated by chance. Therefore, the relative strength of the $W_{C_{ij}D}$ for each j^th C_ij class of evidence in map Xi is BeI C_ij, as indicated in the equations above. N (Cij∩D) represents the number of wells that occur in a subclass; N (Cij) is the number of pixels in a subclass; N (D) is the total number of wells; N (T) is the sum of the pixel domain for a class. Accordingly, for the $W_{C_{ij}D}$, we can calculate the BeI Cij, as shown in Table 1. The values of the Bel range are between 0 and 1 for each subclass.

The groundwater potential probability was the highest for the elevation classes from 19 to 75 m with the highest EBFs (0.565); however, it was the lowest for the elevation classes from 0 to 3 m (EBF = 0.029) (

Table 2. Groundwater factors of the adopted models in the study area.

Topography	No. Pixels in Domain	Domain %	No. Wells	No. Wells %	FR	N (Cij∩D)/N (Cij)	N (D) −N (Cij∩D)	N (T) −N (Cij)	N (D) −N (Cij∩D)/N (T) −N (Cij)	WcijD	Bel
0 to 3	2,883,463	0.29345745	9	0.164	0.558	0.0000031	46	6,942,367	6.62598 × 10−6	0.471	0.029
3 to 7	4,720,324	0.48039952	23	0.418	0.870	0.0000049	32	5,105,506	6.26774 × 10−6	0.777	0.048
7 to 12	1,921,284	0.19553402	16	0.291	1.488	0.0000083	39	7,904,546	4.93387 × 10−6	1.688	0.104
12 to 19	281,047	0.02860288	6	0.109	3.814	0.0000213	49	9,544,783	5.13369 × 10−6	4.159	0.255
19 to 75	19,712	0.00200614	1	0.018	9.063	0.0000507	54	9,806,118	5.50677 × 10−6	9.212	0.565
Slope
0 to 1	4,010,495	0.40868454	26	0.473	1.157	0.0000065	29	5,802,685	4.99769 × 10−6	1.297	0.260
1 to 2	3,568,474	0.36364094	16	0.291	0.800	0.0000045	39	6,244,706	6.24529 × 10−6	0.718	0.144
2 to 3.6	1,889,910	0.19258895	9	0.164	0.850	0.0000048	46	7,923,270	5.80568 × 10−6	0.820	0.164
3.6 to 44	344,301	0.03508557	4	0.073	2.073	0.0000116	51	9,468,879	5.38607 × 10−6	2.157	0.432
Slope ASPECT
Flat	860,817	0.08772049	6	0.109	1.244	0.0000070	49	8,952,363	5.47342 × 10−6	1.273	0.117
North	950,349	0.09684414	4	0.073	0.751	0.0000042	51	8,862,831	5.75437 × 10−6	0.731	0.067
Northeast	857,199	0.08735181	7	0.127	1.457	0.0000082	48	8,955,981	5.35955 × 10−6	1.524	0.140
East	872,576	0.08891878	13	0.236	2.658	0.0000149	42	8,940,604	4.69767 × 10−6	3.171	0.291
Southeast	1,032,724	0.10523847	6	0.109	1.037	0.0000058	49	8,780,456	5.58058 × 10−6	1.041	0.095
South	808,454	0.08238451	3	0.055	0.662	0.0000037	52	9,004,726	5.77475 × 10−6	0.643	0.059
Southwest	972,047	0.09905525	2	0.036	0.367	0.0000021	53	8,841,133	5.99471 × 10−6	0.343	0.031
West	931,399	0.09491307	4	0.073	0.766	0.0000043	51	8,881,781	5.74209 × 10−6	0.748	0.069
Northwest	972,718	0.09912363	5	0.091	0.917	0.0000051	50	8,840,462	5.65581 × 10−6	0.909	0.083
North	1,554,897	0.15844986	5	0.091	0.574	0.0000032	50	8,258,283	6.05453 × 10−6	0.531	0.049
Terrain Roughness Index (TRI)
0.0 to 0.11	602,302	0.06124132	3	0.055	0.891	0.0000050	52	9,232,594	5.63222 × 10−6	0.884	0.169
0.11 to 0.31	869,074	0.08836636	6	0.109	1.235	0.0000069	49	8,965,822	5.4652 × 10−6	1.263	0.241
0.32 to 0.47	3,078,693	0.31303768	15	0.273	0.871	0.0000049	40	6,756,203	5.92049 × 10−6	0.823	0.157
0.48 to0.62	3,743,124	0.3805962	20	0.364	0.955	0.0000053	35	6,091,772	5.74545 × 10−6	0.930	0.177
0.63 to 0.89	1,541,703	0.15675844	11	0.200	1.276	0.0000071	44	8,293,193	5.30556 × 10−6	1.345	0.256
Drainage Density
5.821–53.57	10,128	0.12264619	6	0.109	0.889	0.0005924	49	72,451	0.000676319	0.876	0.210
53.58–101.3	27,266	0.3301808	23	0.418	1.267	0.0008435	32	55,313	0.000578526	1.458	0.350
101.4–149.1	26,864	0.32531273	16	0.291	0.894	0.0005956	39	55,715	0.000699991	0.851	0.204
149.2–196.8	15,257	0.18475642	10	0.182	0.984	0.0006554	45	67,322	0.000668429	0.981	0.235
196.9–244.6	3064	0.03710386	0	0.000	0.000	0.0000000	55	79,515	0.000691693	0.000	0.000
TWI
−8.16 to −2.97	5,514,106	0.56190817	31	0.564	1.003	0.0000056	24	4,299,074	5.5826 × 10−6	1.007	0.341
−2.97 to −0.25	3,128,986	0.31885546	18	0.327	1.026	0.0000058	37	6,684,194	5.53545 × 10−6	1.039	0.352
0.25 to 13.51	1,170,088	0.11923637	6	0.109	0.915	0.0000051	49	8,643,092	5.66927 × 10−6	0.904	0.307
SPI
0 to 0.001	9,290,340	0.94672063	53	0.964	1.018	0.0000057	2	522,840	3.82526 × 10−6	1.491	0.665
0.001 to 0.1	469,594	0.0478534	2	0.036	0.760	0.0000043	53	9,343,586	5.67234 × 10−6	0.751	0.335
0.1 to 34.002	53,246	0.00542597	0	0.000	0.000	0.0000000	55	9,759,934	5.63528 × 10−6	0.000	0.000
Dist to Rivers
0–400	23,750	0.28963768	18	0.327	1.130	0.0007579	37	58,249	0.000635204	1.193	0.294
400–800	28,423	0.34662618	19	0.345	0.997	0.0006685	36	53,576	0.000671943	0.995	0.246
800–1200	17,855	0.21774656	11	0.200	0.918	0.0006161	44	64,144	0.000685957	0.898	0.222
1200–1600	10,759	0.13120892	7	0.127	0.970	0.0006506	48	71,240	0.000673779	0.966	0.238
1600–1900	1212	0.01478067	0	0.000	0.000	0.0000000	55	80,787	0.000680803	0.000	0.000
Lineaments
0–12	3,273,187	0.34086999	21	0.382	1.120	0.0000064	34	6,329,263	5.37187 × 10−6	1.194	0.240
13–24	2,128,686	0.22168155	11	0.200	0.902	0.0000052	44	7,473,764	5.88726 × 10−6	0.878	0.176
25–35	2,256,503	0.23499242	11	0.200	0.851	0.0000049	44	7,345,947	5.9897 × 10−6	0.814	0.163
36–47	1,382,833	0.14400835	9	0.164	1.136	0.0000065	46	8,219,617	5.59637 × 10−6	1.163	0.234
48–59	561,241	0.05844769	3	0.055	0.933	0.0000053	52	9,041,209	5.75144 × 10−6	0.929	0.187
Land cover
Water	1,090,317	0.13192203	4	0.073	0.551	0.0000037	51	7,174,542	7.10847 × 10−6	0.516	0.186
Vegetation	3,325,930	0.40241824	24	0.436	1.084	0.0000072	31	4,938,929	6.27666 × 10−6	1.150	0.415
Urbans	3,848,612	0.46565973	27	0.491	1.054	0.0000070	28	4,416,247	6.34023 × 10−6	1.107	0.399
NDVI
0–66	729,152	0.08520368	2	0.036	0.427	0.0000027	53	7,828,600	6.77005 × 10−6	0.405	0.084
66–131	736,815	0.08609913	4	0.073	0.845	0.0000054	51	7,820,937	6.52096 × 10−6	0.833	0.172
131–174	2,407,709	0.2813483	18	0.327	1.163	0.0000075	37	6,150,043	6.01622 × 10−6	1.243	0.257
174–210	3,376,062	0.39450337	19	0.345	0.876	0.0000056	36	5,181,690	6.94754 × 10−6	0.810	0.167
210–255	1,308,014	0.15284551	12	0.218	1.427	0.0000092	43	7,249,738	5.93125 × 10−6	1.547	0.320
Rainfall
567.4–575.8	2,944,227	0.32047756	24	0.436	1.362	0.0000082	31	6,242,772	4.96574 × 10−6	1.642	0.282
575.9–584.2	3,573,148	0.38893528	13	0.236	0.608	0.0000036	42	5,613,851	7.4815 × 10−6	0.486	0.083
584.2–592.5	1,041,162	0.11332994	3	0.055	0.481	0.0000029	52	8,145,837	6.38363 × 10−6	0.451	0.077
592.5–600.9	855,039	0.09307054	8	0.145	1.563	0.0000094	47	8,331,960	5.64093 × 10−6	1.659	0.285
601–609.3	773,423	0.08418669	7	0.127	1.512	0.0000091	48	8,413,576	5.70507 × 10−6	1.586	0.272

). The zones with the highest slopes (from 3.6 to 44) and lowest slope angles (from 0 to 1) had the highest EBFs (0.43 and 0.26, respectively). The east, northeast, and flat zones had the highest EBFs (0.29, 0.14, and 0.11, respectively); however, the southwest zone had the lowest (EBF = 0.095). Furthermore, the higher TRI areas (from 0.63 to 0.89) had high EBFs (0.256). In the case of the Dd, EBF values = 0.35 were higher in the classes of low density (from 53.58 to 101.3). We assigned the highest EBF values to the two TWI classes: from −2.97 to 0.25 (EBF = 0.35) and from −8.16 to −2.97 (EBF = 0.34). Furthermore, we assigned the SPI classes from 0.0 to 0.001 and from 0.001 to 0.1 EBFs of 0.66 and 0.33, respectively, as well as groundwater probabilities. According to the LD, the values from 0 to 12 had the highest EBFs (0.24). Using the LU/LC, the groundwater potential probability was the highest for the vegetation class (EBF = 0.415), as well as for the areas with high NDVIs (from 210 to 255), with the highest EBF 0.32. The class from 592.5 to 600.9 mm had the highest EBF value (0.285).

We evaluated the EBF weights (Bel values) for the conditioning factors and their classes, and according to the results, the trends for the EBF values for the classes were correlated as for the FR values (Figure 6), with the exception of a few classes in the Dd, NDVI, lineaments, and rainfall.

5.3. Groundwater Potential Mapping and Validation of Built Models

We classified the combined groundwater potential zone (GWPZ) map using the FRs (Figure 7a) combined from twelve thematic maps using the natural breaks approach into five classes: (1) very low (23.67%); (2) low (40.60%); (3) moderate (26.44%); (4) high (8.90%); (5) very high (0.39%). Furthermore, the GWPZ map using the EBF (Figure 7b) that ranged from very low to very high was occupied by 23.60, 41.62, 24.44, 9.57, and 0.78%, respectively. We combined the FR and EBF models to produce the combined FR–EBF model (Figure 7c), which we divided into five classes ranging from very low to very high, and covering areas of 22.89, 39.94, 26.22, 8.98, and 1.97% of the entire area (

Table 3. GWPZ areas.

GWPZ	Area FR %	Area EBF %	Area EBF + FR %
Very low	23.67	23.60	22.89
Low	40.60	41.62	39.94
Moderate	26.44	24.44	26.22
High	8.90	9.57	8.98
Very high	0.39	0.78	1.97

). In the combined FR–EBF model, the areas of the very high and moderate potentiality increased (Figure 8). Based on the Landsat-8 band composites 7, 5, and 3, most of the very high to extreme GWPZs are occupied by vegetation and water resources (Figure 7e,f).

We present the ROC curves for the predicted groundwater map in Figure 8. We can use the AUC to describe the system’s capacity to precisely anticipate both the “groundwater” and “no-existence of groundwater”, which highlights the usefulness of its prediction. The AUC ranges from 0 to 1, with values smaller than 0.5 indicating model integrity, and models with larger values indicating higher accuracy. The AUC results of the models, including the FR, EBF, and combined FR–EBF models, were 0.707, 0.665, and 0.716, respectively. We achieved better accuracy by fusing the two models (the RF–EBF model). Due to the uncertainty around the presence of groundwater in the results, the combined model had better accuracy (AUC = 0.716).

6. Discussion

Water scarcity has now become a major global issue, with substantial implications for the long-term sustainability of water supplies. In the research on spatial groundwater recharge potential zones, researchers emphasize the importance of the efficient planning and management of water resources. According to remote sensing data analyses, GWPZs are controlled by topographic, geomorphic, hydrological, and climatic conditions [28,110,133]. According to the applied FR and EBF models, the rainfall, NDVI, lineament density, and distance to river were positively correlated, and there were more water well locations in the high-altitude areas, which consist of plateaus that are enriched in vegetation and annual precipitation, as well as surface water, near to the Yellow and Zihe Rivers (in the south). The highest well locations were on the east, northeast, and flat slope aspects, which is because the northern hemisphere and north- and east-facing slopes contain more water resources than the south- and west-facing slopes. Less sunlight hits the northern and eastern mountain slopes than it does the southern and western ones. On the north- and east-facing inclines, the transpiration is low despite the high soil moisture levels. The effect is vegetation growth on the northern and eastern faces. In certain places, the increased vegetation enhances the groundwater recharge and surface infiltration [26,87]. The LC/LU types have a substantial impact on the permeability, evapotranspiration, and runoff. The recharge in vegetated areas is substantially lower than that in non-vegetated areas, according to Gee et al. [120]. Furthermore, the recharging is higher in agricultural and grassland areas than on permanent lands, such as shrub and forest areas [121].

The trends of the EBF values for the classes were related to the FR values, according to the evaluation of the EBF weights (Bel values) for the conditioning factors and their classes (Figure 6). The only exception to the similarity were a few classes in the Dd, NDVI, lineaments, and rainfall. The applied FR was better validated than the EBF based on the AUC values of 0.665 and 0.707, respectively. In various studies in the literature, including [75], researchers have demonstrated that the weights generated using the FR and EBF models are consistent [134]. Furthermore, the FR surpasses the EBF techniques, according to numerous studies [28,50]. Combining the two models (FR–EBF model) produced better results than those of each individual model, as validated by the AUC value of 0.716, which was higher than those of the FR and EBF models (Figure 9). Moreover, the combined FR–EBF model increased the areas of very high perspectivity to 1.97, compared with the 0.39 and 0.78 of the FR and EBF models, respectively.

The Yellow River is extensively contaminated by factory waste and sewage from the expanding cities, and hence, it is polluted. The main sources of pollution are industrial, agricultural, and domestic sewage, along with industrial waste gas [135]. Additionally, the increasing wastewater discharge into the river [136] threatens the sustainability of the surface water use. As a result, groundwater has become an appropriate and clean water resource for the long-term sustainability of water supplies in areas with low water quality. Our increasing understanding of the spatial groundwater recharge potential zones highlights the importance of the effective planning and management of water resources [62].

7. Conclusions

In this article, we applied the FR and EBF models, and we integrated them into a hybrid ensemble FR–EBF model that we can use to reveal groundwater potential areas. We prepared, analyzed, and integrated twelve evidential factors. We investigated the spatial recharge based on the topographic, hydrologic, climatic, and land-use elements, and we assessed their potential impacts on GWPZs. According to the findings, the AUC function for the hybrid model (0.716) was reasonably better than those of the individual FR (0.707) and EBF (0.665) models. According to the results, we classified the GWPZs using the natural breaks classifier into five classes that ranged from very low to very high, and that covered areas of 22.89, 39.94, 26.22, 8.98, and 1.97% of the entire area, respectively. The applied modeling techniques and output maps can considerably improve decision makers’ abilities to manage groundwater resources sustainably, to advance efforts to accurately predict the promising areas for extraction and recharge capability. This study presented valuable and better information about the utilize of ensemble FR–EBF model for revealing groundwater potential zones compared with the findings of individual models of FR or/and EBF. Overall, the findings of this study may be applied to the management, planning, and regulation of water resources under various conditions.