Integration of Geospatial Technologies and Fuzzy-AHP Analysis to Assess Groundwater Potential in the Sirwa Massif, Anti-Atlas Region, Morocco

Tazi, M. J.,^1* El Azzab, D.,¹ El Moutaouakkil, N.² and Charroud, M.¹

¹Sidi Mohamed Ben Abdellah University, Faculty of Science and Technology, Intelligent Systems Georessources and Renewable Energy Laboratory, Fez, Morocco

²Mohammed V University Rabat, Faculty of Sciences Rabat, Geology Department, Geosciences, Water and Environment Laboratory, Morocco

*Corresponding Author

Abstract

Climate change has significantly impacted the availability of freshwater in several countries, particularly the African continent. In parallel, the increasing demand for water supply for domestic uses, irrigation and industry is putting more pressure on water reserves, leading to even more water scarcity. Groundwater on the other hand, still holds great potential for water supply reinforcement. aquifers remain largely underexploited in many regions, especially in inaccessible hard rock terrains. In this sense, remote sensing offers a cost-effective way to map groundwater occurrence in these geologically complex terrains. This research aims to delineate groundwater potential in the Sirwa massif, Anti-Atlas, Morocco. initially, a spatial analysis incorporating remote sensing and geophysical data was performed to create nine factors that impact groundwater occurrence including the lithology, fault density, lineament density, node density, slope, distance from rivers, drainage density, topographic wetness index, land use/land cover. Followed by a fuzzy analytic hierarchy process to evaluate the degree of influence of each factor on groundwater occurrence. The groundwater potential map shows that less than 8% of the total study area has high to very high groundwater potential. The remaining surface area (92%) is split between moderate (35.59%) and poor to very poor groundwater potential (55.51%). The predicted map was validated through receiver operating characteristic (ROC) using active borewell data (a success rate of 79.1%). Ultimately, the results shed light on suitable sites for groundwater resources, helping stakeholders to reduce the scale of prospection and make informed decisions for sustainable water management.

Keywords: Arid Regions, GIS, Groundwater Exploration, Multicriteria Analysis, Remote Sensing

1. Introduction

The global climate change crisis is leading to extreme weather events, such as large temporal and spatial variability of precipitation and rising temperatures. These climatic factors are causing more frequent drought periods that can lead to human and environmental disasters, particularly for countries most vulnerable to water scarcity [1]. Morocco is no exception; the country is experiencing the most severe drought in forty years. The frequency of precipitation has dropped significantly, leading to a drastic decline in available water resources, which now stand at less than 600 m³/year/inhabitant, which is considered below the water scarcity threshold [2]. In this sense, groundwater is an essential element in the water supply of these countries, including for domestic, irrigation and other industrial uses. However, the exploration of this resource is complicated and depends on various parameters, including geological, hydrogeological, topographical, climatic and biological factors. What's more, the regions most vulnerable to water are difficult to reach, making ground-based geophysical or hydrogeological studies impossible. It is therefore crucial to develop new cost-effective techniques that combine both conventional data and remote sensing methods to delineate potential groundwater zones in these regions. Over the past few years, new groundwater technologies and approaches have become available, which enables more data to be collected, over large areas with much greater accuracy at relatively low cost and higher speed.

Hence, this can compensate for the lack of information and weak ground observation infrastructure. For instance, numerous researchers have implemented Gis-based multicriteria techniques (MCDM) including, multi-influence factor [3], weighted linear combination [4], and the analytic network process [5] to identify groundwater potential zones (GWPZ) in different regions of the world. However, techniques of this kind require a large amount of groundwater data to build more robust models. While other peer-leading literature [6] and [7] strongly recommended the use of the Hierarchical Analytical Process (AHP) technique for groundwater forecasting. In short, this method is simple to apply, well-structured and flexible. moreover, it enables the user to integrate multiple criteria, conduct pairwise comparisons and consistency verification then synthesize the results. Nevertheless, recent researchers have pointed out the subjectivity of the judgment matrix and the difficulty of quantifying the fuzziness of human judgment using the basic weighting scale (ranging from 1 to 9). To address this situation, the Fuzzy analytic hierarchy (FAHP) was proposed in the literature by [8] as an advanced extension of the classical AHP [9]. It incorporates uncertainty into decision-making by employing fuzzy sets [10] rather than simple numerical weights. Thus, it can eliminate the imprecise judgments of decision-makers. Some of the recent research in the groundwater field, including Moroccan studies [11] and [12] have demonstrated the suitability of the FAHP method for delineating groundwater potential zones especially when the number of factors is high and the relative importance of these factors is close to each other. In contrast, the Sirwa massif has not been the focus of adequate scientific research on groundwater. The existing literature [13] has only provided a limited outline of this massif’s hydrogeological framework. Besides, the complex topographical and geological settings of this region. have made it impractical to accurately map groundwater using ground-based investigations.

Consequently, this research presents a robust methodology to delineate (GWPZ) in the Sirwa Massif. The identification of critical groundwater factors was based on an extensive literature review and expert judgment. Related thematic maps were subsequently extracted from remote sensing and aeromagnetic data and incorporated using GIS-based FAHP analysis. The results of this research will help water managers target the most favorable zones for future groundwater drilling.

2. Study Area

The Sirwa massif is located between the High Atlas and the Anti-Atlas Mountain belts (in the center of Morocco) (30°30'N to 30°45'N, 7°30'W to 7°45'W), at 200 km south of Marrakech (Figure 1). This extended domain is part of a crustal uplift where the highest altitude reaches more than 3140 m, resulting from the contribution of both tectonic movements and the presence of a hot upper mantle anomaly [14]. The area presents a diversity of landscapes marked by large rocky terrains (metamorphic and volcanic rocks) and less extensively, natural vegetation-covered zones, especially in areas irrigated by the main rivers.

Figure 1: Location map of Sirwa massif. (a) Moroccan view, (b) the Anti Atlas belt, and (c) the elevation model of the studied area

Figure 2: The average annual rainfall data recorded during a 36-year monitoring period, in numerous climatic stations located within or around the study area

The Sirwa massif has an arid to sub-arid climate which is influenced by both humid influences from the west and hot dry Saharan influences from the southeast. the temperature varies from winter to summer from -4°C to +46°C. while the average annual Rainfall varies from 260 mm/year in the highest peaks to less than 177 mm/year in the plain areas (Figure 1). From a hydrogeological point of view, the study area is characterized by a complex system of aquifers and springs that are mainly controlled by the geological and tectonic features of the region. The underlying formations are made up of crystalline Proterozoic rocks of plutonic and metamorphic origin strongly altered and fractured [13]. These hard rocks present a mechanical behavior and relatively homogeneous overall hydraulic properties and are characterized mainly by discontinuous permeability. The latter has an important implication from a water exploration point of view since they act as conductors or barriers to underground hydrodynamics [15]. The principal groundwater systems in the Sirwa massif are as follows:

3. Materials and Methods

The current study consisted of four major steps as illustrated in Figure 2 as follows: 1) Geospatial data collection (conventional maps, remote sensing and geophysical data). 2) a quantitative evaluation of groundwater influencing factors using the FAHP algorithm. 3) data transformation and integration through a GIS environment. 4) The final step is ground-truth validation to ensure the reliability of the results. The methodology employed for the current study is presented in Figure 3.

3.1 The Aeromagnetic Data

The Aeromagnetic data covering the Sirwa massif were obtained from the Ministry of Energy Transition and Sustainable Development, the survey was conducted at a flight altitude of 60 m above the earth's surface and 500 m line spacing. The data was processed to transform anomalies and remove asymmetry caused by tilt using pole reduction [16] as shown in Figure 4(a). Subsequently, the horizontal derivative filter in Figure 4(b) was applied to the RTP map to extract the subsurface lineaments as illustrated in Figure 4(c). A detailed description of the techniques can be found in [17]. Field missions were also carried out to study the surface indices of the extracted lineaments.

3.2 The ASTER GDEM (Version 3) and Landsat 8-9 OLI/TIR Level 2

For this study, an ASTER GDEM version 3 and a Landsat 8 OLI scenes (from the period 01/07/2023) covering the study area with no cloud cover were freely acquired from the USGS Earth Explorer ( http://earthexplorer.usgs.gov) and NASA Earth data ( https://www.earthdata.nasa.gov). The technical description of these images is presented in (Table 1).

Figure 3: Flowchart of the methodology adopted for the identification of groundwater potentiality

Figure 4: The Aeromagnetic processing methods: (a) reduction to the pole (b) Total Horizontal derivative (c) the extracted lineaments

Table 1: Characteristics of the ASTER GDEM ([22])

Table 2: >Characteristics of Landsat 8 OLI Bands

The GDEM data provide valuable spectral and geomorphological information very useful for extracting the surface lineaments[18]. On the other hand, Landsat 8 OLI data provides nine spectral bands including a panchromatic band in the visible, near-infrared, and shortwave infrared ranges. These bands have been proven very efficient for detecting and mapping the Land surface variation in several previous studies (i.e., [19]. Before image analysis and interpretations, it is required to apply an atmospheric and radiometric calibration of the Landsat 8 Oli raw scene to correct the image deformation caused by the sensor-earth geometry errors and to calibrate the original image data values to radiance, reflectance, or brightness temperatures.

3.3 Generation of Geospatial Layers

After an extensive review of similar studies with similar hydrogeological and climatic conditions [20][21] and [22] among others, in line with our field observation, nine crucial factors including lithology (Lith), fault density (FD), lineament density (LD), node density (ND), slope (SL), topographical wetness index (TWI), drainage density (DD), distance from the main rivers (DFR) and Land use/Land cover (LULC) were extracted from different sources.

3.3.1 Lithology (Lith)

Groundwater occurrence is dependent largely on the hydrogeological properties of the rock including porosity, permeability and hydraulic conductivity and the aquifer’s storage capacity [23]. Generally, sedimentary rocks such as sandstone and limestone are highly permeable and porous and can transmit significant amounts of water. Thus, they are considered an excellent aquifer material. On the other hand, igneous and metamorphic rocks have medium to low permeability/porosity. Generally, fissures and crack zones provide secondary porosity. however, the circulating water is usually sparse and discontinuous [24]. Figure 5 presents a simplified lithological map of the study area derived from Sirwa's geological map (1/50000). various sedimentary, metamorphic and igneous rocks have been mapped and classified from the Paleoproterozoic age to more recent Meso-Cenozoic to Quaternary periods. These rocks exhibit distinct hydrogeological properties. the Paleoproterozoic to Cryogenian crystalline basement formations are composed of medium to high-grade metamorphic rocks (i.e., Gneiss and Schists) sandstone-quartzites and basaltic flows forming the Tachdamt-Bleida group. these rocks are naturally impermeable. The effective permeability of these rocks is dependent exclusively on the weathering and fracturing process [25]. Conglomerates, tuffs, and pelites from the early Ediacaran period, along with the late Ediacaran volcano-clastic sequence of “The Ouarzazate group”, consisting of coarse-grained conglomerates, rhyolites, ignimbrites, andesites and several granitic intrusions are characterized by discontinuous permeability. Aquifers with continuous permeability and recharge, are restricted to highly altered and fractured zones [13]. The Infra-Cambrian to Cambrian formations related to the Adoudounian series are made up of a thick sedimentary sequence of fractured limestones and dolomites.

Figure 5: Simplified lithological map of the study area

From a hydrogeological perspective, these rocks are permeable, although aquifers underlying these rocks are often deep and difficult to reach. The Sirwa volcano-sedimentary complex is distinguished by voluminous trachy-phonolitic rocks and various pyroclastic flows, which are more or less porous and often permeable, resulting in a multi-layered aquifer system. Finally, Quaternary-age detrital and sedimentary deposits, such as silts, gravels and limestones, are characterized by very high permeability and porosity. These formations are generally found close to main watercourses and develop aquifers with continuous permeability and maintained recharge [13].

3.3.2 Fault Density (FD), Lineament Density (LD) and Node Density (ND)

Surface and deep discontinuities such as faults are the key geological features that control water flow patterns in subsurface aquifers, especially in hard rock terrains. they can act as barriers by blocking or slowing the groundwater flow when lithification and mineral precipitation close the open fractures over time. In contrast, the open faults may act as a favorable pathway for groundwater movement, depending on their length and density[26].

Firstly, the surface structural map in Figure 5 was constructed using a semi-automated approach. This approach is sufficiently rapid and has been proven efficient by many authors [22]. In this technique, a shaded relief filter is applied to an ASTER GDEM image, to enhance the subtle expression of the lineaments. Subsequently, the "LINE algorithm" in PCI Geomatica was used to automatically extract lineaments from the shaded image. The algorithm functions as an edge detection filter, applying the Canny edge detection operator to create a gradient image. Then, the linear and curve-linear features are extracted from the binary edge image and saved as a vector layer [26]. This process is controlled by six parameters as shown in (Table 3). A manual examination and validation of the raw results were performed using existing geological and topographical maps, Google Earth imagery, and ground-truth inspection to distinguish geological structures (such as faults, fractures, etc.) from non-geological objects such as rivers, roads, and railroads as shown in Figure 6 (6(d), 6(e), 6(f) and 6(g)).

Table 3: Parameters used for the PCI LINE algorithm

Figure 6: (a)The extracted lineaments superimposed to faults lines (b) lineaments rose diagram (c) rose diagram of geological faults and (d, e, f and g) examples of local fractures identified in the field

Figure 7: (a) fault density (b) lineament density

Figure 8: Node density variation in the study area

Subsequently,the magnetic lineament and surface faults density maps in (Figures 7(a) and 7(b)) were derived from Figure 4(d) and Figure 6(a) using the density analysis tool within the ArcGIS environment, as defined by Equation 1.

Equation 1

FD: the fault density (or, Lineament density), L_i: the total length of all lineaments (km), i: the number of each lineament in the study area, A: the effective area of lineaments cell grids (km²). Tectonic nodes created by fault intersections and connectivity also affect groundwater flow in hard rock terrains. These nodes increase the effective permeability of the hosting rock and provide localized conductive zones where groundwater can flow and be stored [27]. The node density map presented in Figure 8 represents the number of intersection points/ unit area.

3.3.3 Slope (Sl) and TWI

Topographical indices such as the slope have a major impact on the surficial water infiltration. The gentle slope decreases water flow energy, which provides residence time for rainwater to percolate to the ground and favors groundwater recharge [28]. On the other hand, terrains with a higher degree of slope results present a high runoff coefficient and very low infiltration rates. The slope map (Figure 9(a)) was produced using ArcGIS's Spatial Analyst tool and classified into five classes ranging from flat to gently sloping (0 - 7.33°), moderately sloping (7.33° - 14.14°), strongly sloping (14.14° - 21.21°), steeply sloping (21.21° - 29.85°) and very steeply sloping (29.85°- 66.78°). highest ranks were assigned to low slope values and lower ranks to higher slopes. TWI (topographical wetness index) indicates the spatial distribution of the soil moisture along a hillslope. In general, a higher TWI value indicates a high surface water accumulation and a high probability of the presence of groundwater and vice versa [30].

Figure 9: (a) Slope map (b) TWI map of the study area

Figure 10: (a) Drainage density map (b) distance from the river map of the study area

TWI map (Figure 9(b)) was generated by integrating the slope, flow direction, flow accumulation, and tangent of the slope as input parameters in a GIS environment using Equation 2:

Equation 2

Where: α is the cumulative upslope drainage area per unit contour length and β = Topographic gradient (Slope).

3.3.4 Drainage Density (DD), Distance from Rivers (DFR)

Another important factor that controls groundwater distribution is (DD). It influences the percolation and infiltration rate of water to the underground [31]. Regions with low drainage density favor the infiltration rate, which indicates a favorable condition for groundwater potential. However high to very high drainage density implies high surface run-off and less infiltration. The drainage density map (Figure 10(a)) was achieved through the line density tool in ArcGIS, using the following Equation 3:

Equation 3

∑CI = Sum of all channel lengths in each surface unit (km). At = total area (km²).

Furthermore, (DFR) was taken into account in this study, because unconsolidated alluvial layers are highly permeable and porous and are located essentially near the river courses [32]. These layers allow the groundwater to flow through the pore spaces and form good aquifers. To prepare the distance from the river map (Figure 10(b)), we first produced the hydrographic network using a topographical map of the Sirwa region (at 1/50000 scale) and recent Google Earth images. Thereafter, the Euclidean distance tool in ArcGIS was calculated.

3.3.5 Land use Land Cover (LULC)

In arid regions, vegetation can be a reliable indicator of groundwater occurrence, because it slows the surface runoff allowing more time for the infiltration of rainwater and snowmelt to the underground [33]. Usually, natural plants are found along rivers or in places where groundwater is available at shallow depths with optimal salinity percentages.

Table 4: Theoretical confusion matrix of LULC classification

Figure 11: (a) LULC Classification (b)Training samples using Google Earth images

Here, we've focused mainly on the main earth’s surface features including water bodies, cultivated land, sparse vegetation, bare land and urban areas which provide information on the likelihood of groundwater availability. In this respect, the LULC map in Figure 11(a) was generated using supervised classification from the Landsat 8 image covering the study area. This classification involves collecting a sufficient number of spectral signatures (training samples) for each land cover class using different RGB (Red, Green and Blue) band combinations. These band combinations were proven effective by numerous authors (i.e., [34]. For instance, the (7, 6, 4) combination enabled us to distinguish cropland from sparse vegetation. Additionally, the natural color composite combination (4, 3, 2) was used to identify water bodies and urban areas. Finally, the band combination (7, 6, 2) was employed to delimit bare lands. subsequently, a maximum likelihood classification was employed using ArcGIS to automatically generate an initial brut LULC map.

To assess the accuracy of the LULC map, we carried out cross-validation by comparing several test pixels (test points) from each mapped class with the actual ground-truth data (e.g., Google Earth images and field inspection). Based on the obtained results, the confusion matrix of the current classification was created (Table 4). The Latter summarizes the number of accurate and inaccurate predictions with count values for each sample class and identifies the nature of classification errors [35]. Subsequently, the overall accuracy (OA) was deduced using the Equation 4:

Equation 4

Where TAP = Sum of the correct products for each class (of diagonal). Moreover, we computed the Kappa coefficient (also known as Cohen's kappa) [36] using Equation 5:

Equation 5

Where: m = number of rows and columns in the matrix

n = total number of observations (pixels)

P_ii= observation in row iand column i

P_i+ = marginal total of row i

P_+i = marginal total of column i

This coefficient is commonly employed for accuracy assessment in the remote sensing field [37]. it compares two or more distinct land use classification sets to determine the degree of concordance between raters. Generally, kappa’s values ranging from 0 to 0.20 are considered weak, from 0.21 to 0.40 are acceptable, from 0.41 to 0.80 are substantial, and values exceeding 0.8 indicate excellent accuracy [38].

3.4 Data Integration Using Fuzzy-AHP Method (FAHP)

The FAHP model was used to integrate all groundwater influencing factors to create an accurate GWPZ map. This model, like many other existing MCDM techniques, allows the ultimate goal to be organized into a leveled hierarchical design [39]. The criteria and sub-criteria are evaluated and assigned different scores based on their influence and importance to the desired goal at each level. Finally, using GIS software, the results are converted to a graphical map. Historically, the FAHP approach was developed to deal with subjective judgments of the pairwise comparison matrix in the classical AHP method. Thus, FAHP includes the uncertainty and impreciseness of human judgments by using fuzzy linguistic terms which are quantified to fuzzy numbers [10]. The latter are defined by a membership function μD(x) ∈ [0, 1]. Each value of x is between 0 and 1 and refers to the degree of membership of x in fuzzy set D. Each membership function has a different graphical presentation, including the generalized bell-shaped function, sigmoid function, Gaussian function, trapezoidal function, and triangular function Among these functions, the triangular membership function [8] is mathematically simpler and fits perfectly with the context of this study. It is defined by three parameters, l, m, and u, which correspond to the lower, medium and upper vertices of μD(x) along the X-axis, respectively, as illustrated in Figure 12.

Figure 12: Triangular membership function of the fuzzy logic method

The proposed method proceeds in four stages: Stage 1: construction of the pairwise comparison matrix à = (ã)_n×n of criteria based on Equation 6 and using fuzzy triangular scale as illustrated in Table 2:

Equation 6

Where: indicates the fuzzy judgments of i^th criterion over j^th criterion evaluated by the decision maker, and i, j = 1, 2, 3, …, n

Stage 2: the computation Geometric mean [40] for all the criteria as shown in the Equation 7:

Equation 7

Stage 3: the calculation of the fuzzy weight for each criterium using Equation 8:

Equation 8

Step 4: defuzzification of the weights using the center of area Equation 9:

Equation 9

To evaluate the judgment accuracy of both the pairwise comparison matrix, we computed the consistency ratio (CR) through the following steps [41]: First, we calculated the consistency index (CI) using Equation 10:

Equation 10

Where: λ_max is the maximum eigenvalue of the matrix [42] which is calculated using equation 11 and n is the number of criteria considered in the study

Equation 11

Subsequently (CR) is calculated through Equation 12:

Equation 12

Where: RCI is the random index derived from the table developed by [42].

3.5 Potential Groundwater Zones Mapping and Validation

Statistical validation methods such as the Receiver operating characteristic curve’ (ROC) and ‘area under the curve’ (AUC) are widely used to assess the reliability and credibility of the groundwater potential zones map [43] and [44]. These techniques diagnose the statistical agreement between the final map and the ground-truth data (i.e., borewell yield data). The ROC curve helps identify the best threshold for making a decision. It plots the True positive rate (TTR) against the False positive rate (FPR) which represents the location of wells/springs and the cumulative percentage of the final groundwater potential map. these two parameters are calculated using Equations 13 and 14 [46].

Equation 13

TPrepresent the number of true positive values; FN is the number wrong predicted pixels:

Equation 14

FP represents the number of false positives. TN is the number of true negatives.

Furthermore, the AUC values are computed trough Equations 15:

Equation 15

TC is the number of correctly classified pixels, TD is the number of incorrectly classified pixels, A is the total number of groundwater pixels, and B is the total number of non-groundwater pixels. According to the aforementioned authors, an AUC value between 0.5 and 1 indicates good to excellent accuracy of the model. However, values <0.5 indicate low accurateness of the model.

3.6 Weights Assignment

The outcomes from expert judgments and a thorough literature review enabled us to critically evaluate the contribution of each criterion and sub-criteria to groundwater occurrence. Accordingly, (FD) and (Lith) were identified as the most important factors influencing groundwater availability in the Sirwa massif. The (FD) factor was given slightly higher weight, because faults and fractures are key components for secondary porosity and effective permeability, especially in areas dominated by impermeable hard rock terrains. Subsequently, (DFR) (ND) were allocated the third and fourth higher weights, respectively. Followed by, (LULC), which has a significant impact on drainage density, water retention and recharge rates. Furthermore, (LD) and (SL) were ranked as the sixth and seventh most important factors in this study. The latter controls the surface runoff, infiltration rate and soil drainage capacity. Finally, (DD) and (TWI) have the least influence on groundwater occurrence. Therefore, they were assigned the lowest weights. Subsequently, a pairwise comparison matrix was delicately computed as shown in Table 5. We used a fuzzy scale to rank the criteria and the Saaty (1=Equal importance, 3=Weak importance of one over the other, 5=Essential or strong importance, 7=Demonstrated importance, and 9=Absolute importance to class the sub-criteria. Additionally, the scale allows for intermediate values (2, 4, 6, and 8) [43]. The matrix consistency was further checked using consistency ratio (See Equation 11 in section matrix consistency), taking into account the following indices: λmax = 9.95, n = 9 (number of factors), RCI = 1.45 (based on Saaty’s Random Consistency Indices) [42]. Consistency index (CI)= 0.11 (based on equation 10). As a result, (CR) equals 0.08 < 0.1 (indicating a satisfactory level of consistency).

Table 5: Fuzzy pairwise comparison matrix

Figure 13: (a) Potential groundwater zones map of Sirwa massif, superimposed to the water point’s location. (b) ROC curve of the groundwater potential map. (c) the surface of each class (in Km and %)