Delineation of Groundwater Potential Zones (GWPZs) in a Semi-Arid Basin through Remote Sensing, GIS, and AHP Approaches

Abstract

Groundwater occurrence in semi-arid regions is variable in space and time due to climate patterns, terrain features, and aquifer properties. Thus, accurate delineation of Groundwater Potential Zones (GWPZs) is essential for sustainable water resources management in these environments. The present research aims to delineate and assess GWPZs in a semi-arid basin of San Luis Potosi (SLP), Mexico, through the integration of Remote Sensing (RS), Geographic Information System (GIS), and Analytic Hierarchy Process (AHP). Seven thematic layers (geology, lineament density, land use and land cover, topographic wetness index (TWI), rainfall, drainage density, and slope) were generated in raster format. After the AHP procedure and rank assignment, the thematic layers were integrated using the raster calculator to obtain the GWPZs map. The results indicated that 68.21% of the area is classified as low groundwater potential, whereas 26.30% is classified as moderate. Validation was done by assessing the water residence time data from 15 wells distributed in the study area. Furthermore, the Receiver Operating Characteristics (ROC) curve was obtained, indicating a satisfactory accuracy prediction (AUC = 0.677). This study provides valuable information for decision-makers regarding the conservation and sustainable management of groundwater resources.

1. Introduction

Water is the most dynamic natural resource, playing a meaningful role in human life, economic and social development, and maintaining ecological systems [1,2]. Currently, global water resources are under stress principally sourced by climatic and anthropogenic factors. The water demand has increased due to population growth, rapid urbanization, industrialization, and agricultural activities [3,4]. The groundwater-associated problems are critical in many regions with high population density and economic development [5]. In arid and semi-arid regions, water scarcity has increased significantly due to insufficient surface water [6]. Studies have reported that groundwater resources provide more than 70% of the water supply [7,8] and are gradually depleting at about 545 km³/year due to overexploitation [1,9]. According to [10], most of the groundwater used in dry regions will be fossil water and not sustainable for the future. Excessive pumping has decreased groundwater levels [11].

Groundwater is a fundamental component of the hydrological system which occurs in the subsurface geological formations (aquifers). According to [2,4,12], the occurrence and availability of groundwater depend on the recharge process controlled by several factors such as physiography, lithological composition, drainage pattern, use and land cover, and climatic factors such as precipitation, temperature, evapotranspiration, etc., and geological setting as fractures and lineament features. Therefore, groundwater potential varies meaningfully in space and time, sometimes by a few meters, even within the same aquifer. This confirms the contrast in groundwater potential found from place to place [13,14].

The delineation of groundwater potential zones is complicated due to a lack of common understanding of the several environmental, climatic and topographical factors [1,15]. Moreover, delimitation of potential areas includes appraising many geospatial factors based on scientific methods [16]. Most previous investigations have applied conventional techniques for prospecting groundwater based on geophysics, geology, and hydrogeology. For example, in the San Luis Potosí Valley (SLPV) [17], a hydrogeological study proposed a new model conceptual of SLP aquifer based on vertical electrical sounding and magnetic surveys. Additionally, it analyzed and integrated the land-use changes. They reported that natural and anthropogenic sources had mined the aquifer system. Scarce groundwater recharge was identified and attributed to impermeable rock strata in the SLPV. This was explained by water extraction with water residence time between 3300 and 3600 years, indicating low water renewal. On the other hand [18], a numeric flow simulation model was applied to the San Luis Potosí aquifer. The variations in potentiometric levels over the past 30 years were studied, and future water level decline related to increasing water extraction was modelled. They found a significant drawdown in groundwater levels at the east of the urban zone and proposed a decrease of 30% in the water extraction of wells. Meanwhile, to compensate for this reduction, four new wells nearby of the Sierra de San Miguelito (SSM) were suggested. This is because SSM was considered to be a groundwater recharge zone. Meanwhile, another study realized by [19] analyzed the groundwater flow system in the San Luis Potosí Valley (SLPV) to examine the effect of increasing water and land-use changes on the deep aquifer. They implemented a transient groundwater flow model. They reported that the change from agricultural to urban land use affected the natural recharge areas decreasing the recharge shallow aquifer. Moreover, the groundwater extraction of new zones was recommended. Additionally, ref. [20] presented a methodology of exploration and groundwater prospecting. They applied geophysical methods such as: aeromagnetic surveys, ground magnetic surveys and vertical electrical sounding. Subsurface zones with high probability of being fractured were located and related to permeable areas with groundwater potential. However, more studies were suggested in order to confirm the methodology. Nevertheless, these methods are costly, laborious and time-consuming [21,22,23]. Additionally, development in remote sensing (RS) and GIS technology has supported the delineation of groundwater potential zones in big regions quickly [2,4]. These methods have become very useful and cost-effective, mostly in inaccessible locations, through high-resolution satellite images. This has been accomplished with technological advance, which provides temporal trends based on datasets.

Research has combined the use of RS, GIS, and multi-criteria decision analysis (MCDA) to demarcate groundwater potential zones, using several hydrogeological, geological, and environmental parameters [3,5,21,24]. Analytic Hierarchy Process (AHP) was proposed by [25] to solve complex decision-making through pairwise comparisons. AHP is a useful method for prospecting groundwater potential zones (GWPZs) [8,26]. For example, ref. [27] identified groundwater potential zones, integrating groundwater potential index and water quality index based on AHP. A few studies have used remote sensing, GIS, and AHP for monitoring and assessing GWPZs in Mexico. Studies realized by [28,29,30,31,32] have been reported.

The San Luis Potosí Basin (SLPB), located in a semi-arid region of the State of San Luis Potosí, Mexico, is a highly populated region with significant urban and industrial development [33], and it principally relies on groundwater for survival [34]. Currently, in the SLPB, spatial identification of suitable groundwater potential zones is necessary. Previous studies have identified a large depression cone in the urban area due to overexploitation of groundwater resources [17,19]. Meanwhile, in some parts of the SLPB, water with a residence time greater than 1000 years has been extracted [35]. Therefore, accurately delineating groundwater potential zones is a key element for sustainable groundwater resources management in the basin.

This study aims to identify and delineate freshwater zones for the drinking, industrial, and agricultural water supply for the population in the San Luis Potosí Basin to collaborate on sustainable water resources management in the basin. The methodology combines remote sensing (RS) and GIS with the Analytic Hierarchy Process (AHP). Multiple thematic layers such as geology, lineament density, land use and land cover, topographic wetness index (TWI), rainfall, drainage density, and slope were generated. After the AHP procedure and rank assignment, the thematic layers were integrated using the raster calculator to obtain the Groundwater Potential Zones (GWPZs) map. Additionally, the resulted GWPZ´s map was cross-validated with the water residence time of 15 wells and the receiver operating characteristic curve (ROC).

2. Materials and Methods

2.1. Study Area

The SLPB is situated in Mexico’s Central Alluvial Basins Groundwater Region, in the southeast part of the State of San Luis Potosi [36,37]. The basin covers ≈ 1836 km² with an elevation between 1743 and 2821 m above sea level (masl). It is characterized by the presence of two mountain ranges, the Sierra de San Miguelito (SSM) to the west and the Sierra de Álvarez (SA) to the east, delineating the urban zone and settlements of the San Luis Potosi Valley (SLPV) [33,38] (Figure 1). The climate is semi-arid, with the mean annual evaporation (2038 mm) exceeding the mean annual precipitation (402 mm) and an average annual temperature of 17.5 °C [39,40].

2.1.1. Hydrographic Framework

The SLPB belongs to the hydrologic region of El Salado and is classified as an endorheic basin, lacking any perennial runoff due to the semi-arid conditions [39,40,41]. The runoffs originated from the SSM and SA feed several intermittent rivers such as Españita, Paisanos, El Potosino, La Parada, Mexquitic and the main collector of the basin: the Santiago River [18,39,40]. To capture runoffs, the San José and El Peaje dams were built, which supply water to the urban areas in the SLPV; meanwhile, the El Potosino and Cañada del Lobo dams were built for flood control [42].

On the other hand, the runoffs also feed the principal streams: Grande, La Virgen, Calabacitas, San Antonio, Paraíso and Portezuelo. These last two form the Santa Rita and Laguna Seca lagoons in the rainy season [42,43].

2.1.2. Geological and Hydrogeological Framework

The SLPB is a graben structure filled over the years with volcanic rocks, lava flows, pyroclastic material, and alluvial and lacustrine sediments, similar to other northeastern basins of Mexico [37,44,45]. The oldest rocks belong to the Cretaceous from the Sierra Madre Oriental. They are marine-derived sedimentary rocks that outcrop in the Sierra de Álvarez, mainly composed of limestone with dolomites and shales [44,46,47]. Cenozoic volcanic rocks (including rhyolites and ignimbrites) overlaid these Cretaceous rocks. They presented the Cenozoic tectonic trademark, formed by a series of NW-SE normal faults with a dominant and regional inclination between 15–20° NE [20,45].

A Quaternary granular sequence (alluvial, sands, silts, gravels, and clay) was deposited on top of the volcanic units as basin-fill material [45,48]. This sequence includes a compact, fine-grained clay-sand layer of low hydraulic conductivity (10⁻⁹ m/s) under most flat areas (except at the edges). It allows the presence of two nor-interconnected aquifers: (1) a shallow aquifer and (2) a deep aquifer [37,39,48,49]. Figure 2 shows the conceptual cross-section model of the aquifer system in the SLPB.

The shallow aquifer has a maximum thickness of 250 m and presents textural variations tending toward SSM conglomerates in a clayey matrix, whereas toward the SA, it is silts and sands [45,50]. This aquifer is sensitive to seasonal effects with the presence of contaminants and with a very dynamic behavior [17,38]. In contrast, the deep aquifer is constituted by strongly fractured volcanic rock, and it has an irregular distribution due to a system of pillars and trenches in the valley. This confined type aquifer is bordered by the SSM and the SA and has an approximate maximum thickness of 300 to 350 m [38,43,50].

2.1.3. Groundwater Extraction and Water Supply

According to [45], by 1960, surface resources supplied 59% of water, and the remainder was collected from the aquifer system. However, nowadays, 84% of the water used to supply public-urban, agricultural, industrial, and mining activities comes from groundwater and only 16% from surface water [43,45]. The groundwater extraction is carried out through wells and deep wells, with an active number of 282 and 370, respectively [40,46]. The shallow aquifer has been the most exploited and is almost depleted, so pumping wells with depths up to 1000 m have been used to obtain groundwater from the deep aquifer [19,38]. At present, 96% of the total groundwater volume is contributed by the deep aquifer, and only 4% comes from the shallow one [40]. On the other hand, there is groundwater withdrawal of approximately 153.42 Mm³, with an annual recharge volume of 78 Mm³, defining the aquifer system as overexploited [40,45]. In addition, the highest density of extraction wells is concentrated in the urban area, which has produced a large depression cone [17,19].

2.2. Methodology

The procedure for delineating GWPZs is presented in Figure 3. Multiple parameters have been selected to delineate groundwater potential zones (See Supplementary Materials, Table S1). Based on the Table S1, seven thematic layers maps (Geology, Lineament Density, Land Use and Land Cover, Rainfall, Drainage Density, Slope, and Topographic Wetness Index) were generated in raster format using remote sensing data and GIS (ArcMap v. 10.4.1 and QGIS v. 3.4.0 Madeira).

The thematic layers were projected with WGS84/UTM Zone 14 N datum coordinate system in a resampled resolution of 15 × 15 m. The description for each parameter, the weighting procedure using the Analytic Hierarchy Process (AHP), and the GWPZs map generation are given below.

2.2.1. Geology

The geological characteristics of a zone affect the availability and occurrence of groundwater, which involves the recharge process [51,52,53,54]. A geological map of the study area was elaborated in QGIS, referencing the Geological Mining Chart from the Mexican Geological Service and the map elaborated by [47]. The thematic map was converted from vector format to raster for homogenizing the layer before the weight and rank assignment.

2.2.2. Slope

The slope refers to the inclination of a surface to the horizontal. It is a crucial variable that directly influences infiltration rate and surface runoff in a specific area, determining the groundwater recharge process [55,56,57,58]. A steep slope retains little time the water on the surface, decreasing infiltration and groundwater recharge due to the runoff flowing rapidly. Meanwhile, flat areas have high infiltration rates and low runoff, which favors the recharge [59,60].

The slope map (degrees) was obtained from the Digital Elevation Model (DEM) through the “Slope Tool” of QGIS. The classification was made according to the Soil Terrain (SOTER) model established by [61] (See Supplementary Materials, Table S2).

2.2.3. Lineament Density

Lineaments are straight or nearly straight landforms that are enhanced by the terrain’s permeability and extensively spread over the Earth’s surface [62,63]. They reflect a general manifestation of underground fractures, faults, or joints with intrinsic permeability characteristics and porosity [6,64,65]. The lineaments facilitate the infiltration of water into the subsurface and play a great relevance in the movement and storage of groundwater [52,66].

A Principal Component Analysis (PCA) was carried out in ENVI v. 5.3 to improve the Landsat image [67], followed by the automatic line extraction in PC1 using Geomatica 2016. After the extraction, verifying discontinuities of lineaments [22] was made by the “Split Line Tool” in ArcMap. Finally, the lineament density map (km/km²) was generated employing the “Line Density Tool”.

2.2.4. Drainage Density (D)

The drainage density is defined as the ratio between the total length of the watercourses in a basin and the surface area of the drained basin [51,68,69]. Typically, groundwater occurrence is inversely proportional to drainage density; a high drainage density leads to lower infiltration and minor GWPZs, whereas a low density leads to major GWPZs [58,70,71,72,73,74]. However, the drainage system depends on several variables such as slope gradient, the absorption capacity of soils, rainfall, vegetation cover, climate, topography, and subsurface characteristics [65,75,76]. The knowledge of this variable provides a suitable numerical measure and allows understanding and assessing characteristics of runoff potential, relief, groundwater infiltration, and permeability information [62,71,73].

Drainage density (km/km²) was calculated with the Line Density Tool and the Stream Network using Equation (1) [2].

$$D={\displaystyle \sum }_{i-1}^n\frac{D_i}{A}$$

(1)

where D_i is the total length of the entire streams in stream order i (km), and A is the basin area (km²). Stream Network was obtained by the “Hydrology Tool” in ArcMap, following the procedures for Fill DEM, Flow Direction, Flow Accumulation, Stream Order, and Stream to Feature.

2.2.5. Rainfall

Rainfall is the primary water source in arid and semi-arid regions and has a relevant influence on the hydrological cycle [77,78,79]. In the region, this variable was identified as a dominant factor in the groundwater recharge process [53]. A region with high rainfall levels is considered an area with higher water potential and vice-versa [56]. Data on average annual rainfall from the last 40 years (1981–2021) was obtained from the Climate Hazards Group Infra-Red Precipitation with Station data (CHIRPS) datasets. This source provides a spatial resolution of 0.05°, with a daily, pentadal, and monthly dataset from 1981 to the present and has been widely used in hydrological and meteorological areas with reliable results [80,81]. At last, the corresponding map was elaborated in QGIS.

2.2.6. Land Use and Land Cover (LULC)

The LULC of an area provides information about the environment and determines the dependence on groundwater through water infiltration [22,65]. Initially, a visual analysis of the study area was carried out. The band combinations of Landsat 8 are (1) true-color (bands 4/3/2), (2) color infrared (5-4-3), (3) agriculture (6/5/2), and (4) false-color for urban zones (7/6/4) allowed to generate a general idea about the nature of the site [82].

In satellite images, the pixels can be grouped into one of the LULC classes using classification processes [83,84]. A supervised classification was performed in QGIS using the Semi-Automatic Classification Plugin developed by [85]. In this type of classification, the analyst selects the training pixels, which are used to obtain different LULC characteristics and allow the classification [83]. The Maximum Likelihood algorithm was applied because it is considered the most conventional approach for remotely sensed images [83,86,87].

2.2.7. Topographic Wetness Index (TWI)

The primary data source to perform the TWI is a DEM [88]. This index developed by [89] determines the topographic influence on the hydrologic response of a basin and is the most common index used in hydrological-based studies [88,90,91].

TWI describes the tendency of a cell to accumulate water and allows to locate of favorable areas with concentration and slow runoff and is calculated by the following Equation (2) [22,90]:

$$\mathrm{TWI}=ln\left(\frac{\alpha }{tan \left(b\right)}\right)$$

(2)

where TWI is Topographic Wetness Index, α is the specific catchment area, and b is the local slope. The TWI function from “Terrain Analysis-Hydrology” of the SAGA-GIS Tool was employed for the corresponding thematic map generation in QGIS.

2.2.8. Weight Calculation Using AHP

The AHP is a Multi-Criteria Decision-Making (MCDM) method developed by [25], which is widely accepted for analyzing spatial decision problems in the field of natural resources management, including groundwater issues [56,92]. AHP method is used to assess the weight of different layers [12]. The first step is the construction of a Pairwise Comparison Matrix (PCM) using Saaty’s scale (1–9) of relative importance [93] (

Table 1. Saaty’s scale of relative importance.

Scale	Definition	Explanation
1	Equal importance	Two activities contribute equally to the objective
3	Moderate importance of one above other	Experience and judgment strongly favor one activity over another
5	Essential or strong importance	Experience and judgment strongly favor one activity over another
7	Very strong importance	An activity is strongly favored, and its dominance demonstrated in practice
9	Extreme importance	The evidence favoring one activity over another is of the highest possible order of affirmation
2,4,6,8	Intermediate values between the two adjacent judgments	When compromise is needed

). The weights were allocated based on previous studies shown in Table S1 and the study made by [53] in the region.

After comparing each layer based on their relative importance, a PCM for seven variables was constructed (

Table 2. Pairwise Comparison Matrix of seven variables for AHP method.

Variable	Geology	Slope	Lineament Density	Drainage Density	Rainfall	LULC *	TWI *
Geology	1	2	3	4	7	8	9
Slope	0.500	1	3	2	4	5	8
Lineament Density	0.333	0.333	1	2	5	3	4
Drainage Density	0.250	0.500	0.500	1	3	3	4
Rainfall	0.143	0.250	0.200	0.333	1	4	5
LULC *	0.125	0.200	0.333	0.333	0.250	1	3
TWI *	0.111	0.125	0.25	0.250	0.200	0.333	1
SUM	2.462	4.408	8.283	9.917	20.450	24.333	34

Note(s): * LULC = Land Use and Land Cover; TWI = Topographic Wetness Index.

). The second step of this method is to calculate the normalized weights to reduce the associated subjectivity. A sum of the values in each column was carried out by Equation (3) [11] and is also shown in Table 2.

$$L_{ij}={\displaystyle \sum }_{n=1}^n{C}_{ij}$$

(3)

where L_ij is the total column value of the PCM and C_ij is the variable employed for the analysis.

Subsequently, all the column values were divided by the sum of the column to create the Normalized Pairwise Comparison Matrix (NPCM) [4,11]. The normalized weight (N_Wt) of each variable was obtained by averaging all the values of the corresponding row in the NPCM (

Table 3. Normalized Pairwise Comparison Matrix.

Variable	Geology	Slope	Lineament Density	Drainage Density	Rainfall	LULC *	TWI *	Total	Nwt
Geology	0.406	0.454	0.362	0.403	0.342	0.329	0.265	2.561	0.37
Slope	0.203	0.227	0.362	0.202	0.196	0.205	0.235	1.630	0.23
Lineament Density	0.135	0.076	0.121	0.202	0.244	0.123	0.118	1.019	0.15
Drainage Density	0.102	0.113	0.060	0.101	0.147	0.123	0.118	0.764	0.11
Rainfall	0.058	0.057	0.024	0.034	0.049	0.164	0.147	0.533	0.08
LULC *	0.051	0.045	0.040	0.034	0.012	0.041	0.088	0.312	0.04
TWI *	0.045	0.028	0.030	0.025	0.010	0.014	0.029	0.182	0.03

Note(s): * LUCL = Land Use and Land Cover; TWI = Topographic Wetness Index; N_Wt = normalized weight.

). The sum of all normalized weights is always 1 [4,60].

During the AHP method application, a certain level of inconsistency could arise because it is conducted through subjective or personal judgments [60,94]. To evaluate the accuracy, the Consistency Ratio (CR) was calculated.

At first, each column of PCM was multiplied by the respective variable weight. Then, the sum of the rows allowed us to obtain the weighted sum value. After that, a division between the weighted sum value and the variable weight was made, obtaining a λ value [4]. The maximum eigenvalue (λ max) was calculated using Equation (4) [25]:

$$\lambda \max =\frac{C1+C2+C3\dots Cn}{n}$$

(4)

where C1 to Cn is the λ value and n is the number of criteria. In this research, a λ max value of 7.607 was obtained (See Supplementary Materials, Table S3)

Afterward, the Consistency Index (CI) value was obtained by the following Equation (5) [25]:

$$CI=\frac{\lambda \max -n}{n-1}$$

(5)

where λ max denotes the maximum eigenvalue of the judgment matrix, and n is the number of criteria. A Consistency Index value of 0.101 was obtained in this study.

Finally, the Consistency Ratio (CR) was obtained by Equation (6):

$$CR=\frac{CI}{RI}$$

(6)

where CI is the Consistency Index, and RI is the Random Consistency Index established by [93] (

Table 4. Random Consistency Index (RI) values for n variables.

n	1	2	3	4	5	6	7	8	9	10
RI	0.00	0.00	0.58	0.89	1.12	1.25	1.32	1.40	1.45	1.49

). In this study, the RI value is 1.32 since seven variables were employed.

If the value of the Consistency Ratio (CR) is ≤0.10, the inconsistency is acceptable. On the other hand, judgments should be revised if the CR is ≥0.10. In this study, a valid CR value of 0.076 was obtained.

2.2.9. Mapping Groundwater Potential Zones (GWPZs)

After evaluating weights using the AHP method, and to standardize each raster, rating values of 1–5 (very low, low, medium, high, and very high, respectively) were assigned for each feature class and normalized with the geometric mean criteria [95,96]. The rating values reflected the suitability of groundwater potential. For the geology variable, the rank assigned was carried out according to other authors, as well as the hydraulic properties of the geological units [4,37,52,58,60,62,97]. The normalized weight for each layer and the normalized rank for each class are presented in

Table 5. Assigned normalized weights and ranks of thematic layers.

Variable	Units	Normalized Layer Weight	%	Classes	Class Rank	Normalized Class Rank
Geology	-	0.37	37	Alluvium	5	0.18
				Basalt	2	0.07
				Breccia	2	0.07
				Limestone	4	0.14
				Conglomerate	3	0.11
				Ignimbrite	1	0.03
				Lacustrine	3	0.11
				Latite	1	0.03
				Shale	3	0.11
				Rhyolite	1	0.03
				Tuff	2	0.07
Slope	degree	0.23	23	0–2	5	0.33
				2–8	4	0.26
				8–15	3	0.20
				15–30	2	0.13
				>30	1	0.06
Lineament Density	km/km2	0.15	15	0–0.39	1	0.06
				0.39–0.78	2	0.13
				0.78–1.17	3	0.20
				1.17–1.56	4	0.26
				1.56–1.95	5	0.33
Drainage Density	km/km2	0.11	11	0.05–3.64	5	0.33
				3.64–7.23	4	0.26
				7.23–10.82	3	0.20
				10.82–14.41	2	0.13
				14.41–18	1	0.06
Rainfall	mm/year	0.08	8	302–342	1	0.06
				342–381	2	0.13
				381–421	3	0.20
				421–461	4	0.26
				461–501	5	0.33
Land Use and Land Cover	-	0.04	4	Urban area	1	0.06
				Bare ground	2	0.13
				Water	3	0.20
				Vegetation	4	0.26
				Agricultural land	5	0.33
TWI	%	0.03	3	2.22–6.85	1	0.06
				6.85–11.48	2	0.13
				11.48–16.10	3	0.20
				16.10–20.73	4	0.26
				20.73–25.35	5	0.33

.

The GWPZ’s map was obtained after the integration of all parameters according to their importance through the Groundwater Potential Index (GWPI) in Equation (7) [2,4,12]:

$$GWPI={\displaystyle \sum }_{w=1}^m{\displaystyle \sum }_{j=1}^n\left(Wj\times Xi\right)$$

(7)

where Wj is the normalized weight of the jth variable and Xi is the normalized weight of ith class of the variable. The corresponding integration was carried out using the Raster Calculator Tool in QGIS.

2.2.10. Validation of Results

A validation approach for the GWPZ’s map was based on the available water residence time from 15 groundwater wells located in the study area. A projection of the points was made in which the water residence time (age) was classified into five classes.

In addition, the Receiver Operating Characteristic (ROC) curve was also applied. As stated by [98], this curve is obtained by considering the cumulative percentage of the area (x axis) and the cumulative percentage of the number of wells (y axis). The Area Under the Curve (AUC) was calculated to determine the prediction accuracy [56,99]. This curve is a widely used metric for predicting model accuracy in classification issues [100]. In several studies, this method has been used to assess the accuracy of the groundwater potential mapping using a GIS-based AHP method [6,98,101,102,103,104,105].

3. Results and Discussion

3.1. Description of Thematic Layers

The results obtained for seven thematic layers are described and related to their influence on GWPZs. The area coverage of each class is shown in

Table 6. Area coverage of each class for seven thematic layers.

Classes	Area (km2)	Area (%)
Geology
Alluvium	887.95	48.4
Basalt	5.04	0.3
Breccia	2.14	0.1
Limestone	174.59	9.5
Conglomerate	61.11	3.3
Ignimbrite	374.7	20.4
Lacustrine	7.47	0.4
Latite	78.82	4.3
Shale	1.20	0.1
Rhyolite	233.37	12.7
Tuff	9.44	0.5
Slope (°)
0–2	781.45	42.55
2–8	502.66	27.37
8–15	235.75	12.84
15–30	268.25	14.61
>30	48.32	2.63
Lineament Density (km/km2)
0–0.39	309.4	16.85
0.39–0.78	831.2	45.28
0.78–1.17	1.22	0.07
1.17–1.56	621.95	33.88
1.56–1.95	72.1	3.93
Drainage Density (km/km2)
0.05–3.64	1513.41	82.41
3.64–7.23	233.74	12.73
7.23–10.82	56.41	3.07
10.82–14.41	20.82	1.13
14.41–18	12.05	0.66
Rainfall (mm/year)
302–342	152.28	8.29
342–381	803.01	43.73
381–421	578.81	31.52
421–461	219.23	11.94
461–501	82.8	4.51
Land Use and Land Cover
Urban area	309.4	16.85
Bare ground	831.2	45.28
Water	1.22	0.07
Vegetation	621.95	33.88
Agricultural land	72.1	3.93
Topographic Wetness Index (%)
2.22–6.85	817.87	44.54
6.85–11.48	596.63	32.49
11.48–16.10	391.31	21.31
16.10–20.73	26.98	1.47
20.73–25.35	3.64	0.20

and addressed during the description of the layers.

3.1.1. Geology

Geology plays a key role in groundwater potential because water infiltration and percolation are controlled by the hydraulic properties of the rock [106,107]. According to the importance, this variable is ranked 1st with a normalized weight of 0.37 (Table 3). The geological classification of the study area is presented in Figure 4. Sedimentary rocks (e.g., limestone, lacustrine, shale, conglomerate, and breccia) are widespread exposed along the SA and belong to the formations: Indidura, La Peña, El Abra, Tamabra, and Cuesta del Cura [17,20,50]. On the other hand, alluvium encompasses almost half of the total basin area (48.4%) (Table 6). This material is common in basins within mountain ranges, filling the central part (Fitts, 2013), and is considered suitable for groundwater occurrence (Table 5) [94,108,109]. Towards the SSM, volcanic rocks such as ignimbrite, rhyolite, and latite outcrop, occupy 20.4, 12.7, and 4.3% of the basin area, respectively (Table 6). These types of rocks most present low permeability and porosity (<1%) [52,66] and are considered with a very low and low groundwater potential (Table 5) [30,109].

3.1.2. Slope

The terrain slope is a crucial parameter in groundwater potential and has an inverse relation with infiltration [110,111]. The slope is the 2nd most important variable in this research, with a normalized weight of 0.23 (Table 3). The corresponding map is presented in Figure 5. Flat areas (0–2°) are in the center, representing 42.55% of the basin area (Table 6). These zones are associated with very high groundwater potential due to the high infiltration rates and low runoff (Table 5) [112,113]. In contrast, steeper slopes (15–30° and >30°) are found surrounding the flat area, most of them located towards the SSM, which is characterized as a complex terrain [114]. In this area, the slope positively correlated with runoff, indicating high runoff and fewer infiltration rates [53]. Due to the aforementioned, steeper slopes are unsuitable for groundwater occurrence and are classified with very low-low groundwater potential (Table 5) [2,115].

3.1.3. Lineament Density

Lineaments manifest underlying structural features such as fractures, faults, or joints over an area [6,116]. A suitable groundwater potential zone is characterized by a high density of lineaments, while low potential areas have a lower density [12,94]. In this study, lineament density is ranked 3rd with a normalized weight of 0.15 (Table 3). The corresponding map is shown in Figure 6. A density range of 0.39–0.78 km/km² covers almost half of the basin area (45.28%) and is classified as low groundwater potential (Table 5 and Table 6). Meanwhile, higher densities cover 33.88% (1.17–1.56 km/km²) and 3.93% (1.56–1.95 km/km²) of the area and are classified as high and very high potential, respectively (Table 5 and Table 6). Most of the high densities are situated approaching the SSM, which has a volcanic composition, and lineaments allow us to understand groundwater occurrence in these areas with poor primary porosity [117]. This frontier is characterized by well-exposed normal faults, classified as a “domino system” because they show uniform fault dip direction and similar bed dip angles [118,119]. On the other hand, the volcanic rocks that compose this Sierra often present fractures, giving a secondary porosity [66].

3.1.4. Drainage Density (D)

The amount of surface water recharge into groundwater is controlled by the features of the drainage system [120]. The drainage density measures how well a basin is drained by stream channels and is controlled by several factors such as slope angle, geology, vegetation, and soil type [75,121]. This variable is the 4th most important in this study, with a normalized weight of 0.11 (Table 3). Figure 7 shows that the lowest density (0.05–3.64 km/km²) encompasses almost all of the basin, covering 82.41% (Table 6). However, most of this area corresponds to the SSM and SA, which present geological and slope features that allow high runoff rates and less infiltration. A negative correlation of drainage density with slope and hydraulic properties of geological units has been reported [76,122]. The runoff from these mountain ranges feeds several intermittent rivers in the SLPV, corresponding to the greatest drainage densities and thus low groundwater potential (Table 5). After extreme rains, floodings originate and affect the urban zone where the recharge is limited [43].

3.1.5. Rainfall

Rainfall is the primary hydrological source of groundwater recharge in semi-arid areas and tends to present patterns of high frequency-amount and short duration [59,123,124]. In this research, rainfall is the 5th essential variable with a normalized weight of 0.08 (Table 3). In the SLPB, the rainy season encompasses from June to September [125], and according to [46], the complex orography determines the rainfall patterns over the area. As shown in Figure 8, the rainfall distribution clearly indicates that towards the Sierras, high levels of rain (from 421 to 501 mm) occur, classifying them as high-very high groundwater potential (Table 5). By contrast, the general rainfall in the valley is lower (from 302 to 381 mm) due to the low humidity rates and high temperatures and evaporation (2038 mm) [40,46], classifying those levels as low-very low (Table 5).

3.1.6. Land Use and Land Cover (LULC)

The spatial distribution and characteristics of the land are assessed by a LULC analysis [92]. Long-term and seasonal changes in LULC significantly affect the hydrological dynamics such as groundwater storage and recharge [10,126,127]. In this study, LULC is ranked 6th with a normalized weight of 0.04 (Table 3). Over time, the SLPB has presented several land-use changes (e.g., mining, urban, agricultural and industrial) [128]. Figure 9 shows the LULC map for the study area derived from the Supervised Classification. Bare ground encloses 45.28% of the basin and generally refers to areas with rocky outcrops and a lack of vegetation; thus, its groundwater potential is low (Table 5 and Table 6) [2,58,59]. Vegetation represents 33.88% of the area (Table 6) and comprises chaparral, natural grassland, scrub, and in the highest altitudes of the SSM, a pine-oak forest ecosystem [18,129]. This class was ranked as having high groundwater potential (Table 5). The urban area encompasses 16.85% of the basin, and between 2007 and 2020, it has grown by 7.82% [33]. Agricultural areas are located principally in the northeast, where the main crops are corn, tomato, oats, nopal, and bean [18]. Finally, water represents only 0.07% and corresponds mainly to the dams distributed in the southwest zone of the basin.

3.1.7. Topographic Wetness Index (TWI)

TWI represents possible groundwater infiltration into the subsurface, depending on topographical characteristics and their impact on the surrounding terrain [88]. A higher TWI value represents a high groundwater potential and vice-versa; therefore, there is a positive correlation [106,109]. TWI is ranked 7th with a normalized weight of 0.03 (Table 3). The corresponding thematic map is presented in Figure 10. The lowest range values of 2.22–6.85 and 6.85–11.48 covers almost all of the basin area with 44.54% and 32.49%, respectively (Table 6). Those results are in accordance with the irregular topography of the flanks, classified with very low and low groundwater potential (Table 5); by contrast, higher values are present towards the center. According to [130,131], higher TWI can be found in plain regions with recurrent flooding in a basin.

3.2. Mapping Groundwater Potential Zones (GWPZs)

Groundwater occurrence in semi-arid environments is variable in space due to climate factors (low rainfall and high temperatures and evaporation rates) [15,30,132], geological characteristics of the basin [54], and LULC changes [133], among others. Thus, groundwater potential zone delineation in these environments turns essential because it enables more precise hydric resource research and a better understanding of their long-term use.

The groundwater potential map of the SLPB was generated by integrating all the thematic layers employing Equation (7), GWPI values ranged from 0.0785 to 0.2594. The final map was classified into four classes according to the groundwater potential: Very Low, Low, Moderate, and High (Figure 11), representing 3.54, 68.21, 26.30, and 1.95% of the area, respectively (

Table 7. Groundwater Potential Zones classification and their area coverage.

Classification	Area (km2)	Area (%)
Very Low	64.94	3.54
Low	1252.17	68.21
Moderate	482.88	26.30
High	35.72	1.95

).

Very low potential zones are randomly distributed in the northwest and southeast parts of the basin (Figure 11), corresponding mainly to the volcanic formations of San Miguelito Rhyolite, Cantera Ignimbrite, and Portezuelo Latite, which are made up of interbedded layers of volcanic rocks with low permeability and porosity [52,66]. Meanwhile, low potential zones cover almost the entire basin area (Figure 11 and Table 7). In the central part, very low and low zones correspond with the location of the depression cone reported by [17,19] in the shallow aquifer due to the population and industrial growth that have been exacerbating the water scarcity in the study area [17,38]. On the other hand, increasing impervious surface with urbanization tends to reduce infiltration and groundwater occurrence [134,135,136]. As mentioned above, the urban area in the SLPV has grown by 7.82% annually [33].

Moderate zones are predominantly located in the flank areas, as shown in Figure 11. According to [17], the SA presents different hydrogeological features. Towards the Gulf of Mexico, reef and platform origin rocks predominate, considered suitable for water storage and transfer. This agreed with the results obtained because moderate to high potential zones are within the SA similarly. In contrast, groundwater potential decreases as one approaches the SLPV due to facies with significant clay content, performing as a hydrological barrier and decreasing groundwater occurrence [17,46], which also agreed with the very low to low potential areas present in the center.

Towards the SSM, although the GWPZ’s map shows zones with moderate potential, a recent study demonstrates that these areas are not working optimally as a groundwater recharge zone due to high evapotranspiration rate, soil type, terrain characteristics, and geological properties [53]. Besides, in the case of the fractures, at the local scale, many of them are sealed by minerals of hydrothermal origin [17,38,137], acquiring low permeability characteristics and limiting groundwater recharge potential.

In the central part of the basin, some moderate potential zones could be attributed to agricultural activities. Research has reported that agricultural practices positively correlate with groundwater recharge because irrigation results in an increase in groundwater occurrence in many semi-arid regions; however, it depends on its patterns of evapotranspiration and infiltration rates [138,139].

Groundwater potential is variable in space and time, even by a few meters in the same basin, which confirms the heterogeneity found in this study. The results generated from this model reflect that the main variables such as geology and slope exert greater control over the groundwater occurrence in the study area.

Validation of Results

Results have been validated with the help of water residence data of 15 wells distributed in the study area, the age in years was classified into five classes (Figure 11).

As a first validation, the map of located wells was overlain on the main map of GWPZs to verify the accuracy of the AHP method.

Table 8. Water residence age of wells and their groundwater potential.

Well ID	Coordinates (UTM 14 N)		Water Residence Time (Years)	Expected	Location	Result
	X	Y
1	291,884	2,466,140	1000	Very low to low	Low	Agree
2	299,813	2,470,005	1400	Very low to low	Low	Agree
3	309,847	2,467,525	3300	Very low	Low	Disagree
4	307,922	2,463,513	3300	Very low	Low	Disagree
5	293,651	2,453,135	6500	Very low	Low	Disagree
6	300,031	2,452,851	1300	Very low to low	Low	Agree
7	297,702	2,451,805	1300	Very low to low	Low	Agree
8	301,995	2,459,686	<55	Moderate	Moderate	Agree
9	297,849	2,458,327	<55	Moderate	Low	Disagree
10	290,738	2,450,551	<55	Moderate	Moderate	Agree
11	291,091	2,448,932	<55	Moderate	Moderate	Agree
12	298,733	2,444,393	<55	Moderate	Moderate	Agree
13	303,055	2,444,473	<55	Moderate	Moderate	Agree
14	311,302	2,448,205	<55	Moderate	Moderate	Agree
15	308,689	2,459,737	<55	Moderate	Low	Disagree

presents the description of the wells together with the groundwater expectation before the output map, and the actual result obtained.

Water residence time ranges from <55 to 6500 years. The spatial distribution of the wells showed that 60% of the points are located in areas with low groundwater potential and encompass water residence time ≥ 1000 years. Furthermore, the spatial location of wells in the urban zone corresponds with the reported by [17,19], where the oldest water coincides with the depression cone in the San Luis Potosi aquifer. In contrast, 40% of the wells are located in moderate potential areas and include residence time of <55 years. An accuracy percentage was calculated between the expected groundwater potential and the actual one, obtaining a 67% of agreement.

As a second validation, the Receiver Operating Characteristic (ROC) curve and the Area Under the Curve (AUC) were used to validate the prediction accuracy of the groundwater potential map. The results showed a value of 0.677, corresponding to a prediction accuracy of 67% (Figure 12).

Similar results have been reported in other studies which applied the AHP method for groundwater potential zonation [4,140,141]. According to the AUC value obtained, the GIS-based AHP model’s prediction accuracy could be accepted as satisfactory (0.6–0.7) [142]. However, this result can be improved by adding more wells, flow rates, and depth data. In addition, water residence time was a valuable tool to evaluate the GWPZ’s map because it gives an idea of how limited the hydric resource is in the study area.

3.3. Limitations of the Study

Although this research followed a carefully systematic analysis, some limitations are present. (1) The data availability of wells is scarce in the study area, resulting in a lack of information for a more robust validation. (2) AHP is considered an empirical method; the weight and rank assignations are based on an expertise point of view. Although the model’s accuracy was accepted as satisfactory, a more precise map could be developed by rearranging the ranks or adding more variables influencing the groundwater occurrence (e.g., geomorphology, curvature, roughness, soil type). (3) Other Multi-Criteria Decision Making techniques could be applied for better results. Among the most employed for identifying GWPZs are Multi Influencing Factor Analysis (MIF) [14,72,98,143], Analytical Networking Process (ANP) [103,144,145], and Fuzzy-Analytic Hierarchy Process [113,131,146,147]. (4) Thematic maps often generalize the nature of the site, so a variation on a local scale could be found. This must be considered in decision making.

4. Conclusions

In this study, freshwater potential zones in the SLPB were identified and analyzed with remote sensing, GIS, and Analytic Hierarchy Process (AHP) using seven thematic layers such as geology, lineament density, LULC, TWI, rainfall, drainage density, and slope. The results accomplished from this methodology revealed that more than 50% of the basin’s surface area showed low groundwater potential zones situated in the central part, where SLP city is situated. This is explained by the changes from agricultural to urban land use. This area was considered a recharge zone toward a shallow aquifer but was altered by land use change. This has led to the over-extraction of water from the deep aquifer. Currently, a large drawdown cone is located. Very low potential zones (3.5%) are located in the northwest and southeast parts of the basin. This is attributed to impermeable volcanic rock formations that confine the recharge. Meanwhile, high (1.95%) and moderate (26.3%) potential zones are located in the flank areas and within the SA. This can be explained by reef and platform origin rock formations in these zones being suitable for water storage. On the other hand, some moderate potential zones situated in the central part of the SLPB could be explained as agricultural activity increasing groundwater occurrence. The results were validated by the water residence time of 15 wells ubicated within the SLPB and ROC curve. It observed a correlation between low potential zones and the oldest water. This could indicate a scarce recharge or slow water renewal in zones. 60% of the wells were located in low groundwater potential zones and were related to water residence time ≥ 1000 years. Meanwhile, 40% of the wells were located in moderate potential areas and were related to water residence time < 55 years. This was confirmed with the ROC curve, which showed a value of 0.677, indicating a prediction accuracy of 67% (satisfactory). The identification and delineation of freshwater potential zones based on remote sensing (RS), GIS, and the Analytic Hierarchy Process (AHP) have not been determined in the SLP basin. This study could contribute to suitable groundwater management. The spatial location of very low, low, moderate, and high potential zones could help develop strategies to mitigate and prevent groundwater contamination and suitable groundwater resources management.