The Feija Plain is located in southeastern Morocco within the Middle Drâa Valley (MDV). This region, historically characterized by subsistence agriculture and seasonal pastoralism, is ecologically fragile and climatically arid (Schulz and Manfred, 2013). It lies along the western flank of the Drâa River near the city of Zagora and forms part of an interconnected oasis system extending over 200 km (Karmaoui and Adil, 2016).

Over the past 2 decades, Feija has undergone a rapid transformation toward intensive, export-oriented agriculture (Lamqadem and Pradhan, 2019). In particular, watermelon farming has emerged as the dominant land use, incentivized by government-subsidized drip irrigation systems and the widespread availability of high-yield hybrid seeds from international agribusinesses (Fico, 2022; Silva-Novoa Sánchez, 2024). The region’s warm winter climate allows for two harvests per season prior to the main production peaks in competing regions, granting it a comparative advantage in the national market (Bossenbroek et al., 2023).

Hydrogeologically, the Feija Plain is primarily underlain by shallow, unconfined phreatic aquifers hosted in Quaternary alluvial deposits, which have traditionally been accessed through hand-dug wells (Aoubouazza and Elmeknassi, 1996; Klose, 2013; Bassin Hydraulique de Souss Massa et de DraaABHSMD Agence du, 2014). However, the intensification of irrigated agriculture, particularly for high-demand crops such as watermelon, has triggered a dramatic increase in groundwater abstraction. As a result, deep boreholes, often exceeding 100 m in depth, are now being drilled into older, confined aquifer systems embedded within Ordovician and Cambrian geological formations. These deep aquifers exhibit low recharge rates due to their lithological composition and structural confinement, rendering them highly vulnerable to overexploitation. The cumulative impact has been a significant decline in piezometric levels, raising critical concerns over long-term aquifer sustainability, ecological degradation, and the socioeconomic viability of agricultural livelihoods in the region (Moumane et al., 2021; 2024).

The area is also traversed by ephemeral watercourses, including Oued Feija and Oued Boutious, which become active during high-intensity rainfall events. During the 2024–2025 hydrological season, both streams experienced significant flash floods Figures 2, 3, with at least five episodic flow events recorded. A particularly destructive storm occurred in early May 2025, when torrential rains and hail struck the Feija Plain. This event devastated nearly 4,000 to 5,000 tons of watermelon crops, only 1 week before the main harvest. According to field reports, more than 700 farmers lost 90%–95% of their production, leading to severe financial distress and widespread concerns about the repayment of agricultural loans (Freshplaza, 2025). Many cultivation areas were rendered inaccessible due to waterlogged tracks, and no water harvesting or artificial recharge systems were in place to capture the millions of cubic meters of floodwater lost to the desert.

The dataset used in this study integrates both remote sensing-derived features and in situ observations to support accurate classification of groundwater potential zones. The in situ component consists of field-based groundwater potential assessments or borehole productivity data, which serve as ground-truth labels for supervised machine learning. Predictor variables were extracted from multiple geospatial and satellite data sources, including stream distance, drainage density, soil permeability, rainfall, topographic wetness index (TWI), geomorphology, lineament density, curvature, slope, elevation (DEM), and the Normalized Difference Vegetation Index (NDVI). These features were obtained from processed satellite imagery (e.g., Sentinel-2), digital elevation models (e.g., SRTM), and thematic layers such as geological or hydrological maps. All spatial layers were resampled and reprojected to a common spatial resolution and coordinate system to ensure consistency. After preprocessing and integration, the dataset was formatted as a structured table where each entry corresponds to a geospatial unit associated with the input variables and a labeled groundwater potential class. Summary statistics for all input features, including their mean, standard deviation, minimum, and maximum values, are presented in Table 1.

To improve model robustness and reduce redundancy in the input space, feature selection was conducted using both linear and nonlinear dependency analyses. First, a Pearson correlation matrix (Figure 4) was computed to identify pairs of features with high collinearity. Features exhibiting strong correlations (i.e., r > 0.9 ) were examined, and one feature from each correlated pair was removed to avoid multicollinearity, which can distort model learning and increase variance.

In parallel, mutual information (MI) was computed between each input feature and the target groundwater potential class to capture nonlinear dependencies Figure 5. Features with low MI scores, indicating weak predictive relevance, were excluded from the final feature set. This dual approach allowed for the retention of variables that were both independent and informative, thereby optimizing the input space used for model training and reducing computational complexity.

In this study, a combination of traditional ML algorithms and modern DL architectures was employed to classify groundwater potential zones based on a set of geospatial and environmental features. The chosen models were selected to provide a comprehensive comparison between interpretable, efficient classifiers and more complex, representation-learning-based models. Specifically, the ML models used include CatBoost and AdaBoost. These algorithms are known for their robustness, generalization capabilities, and suitability for structured data.

To complement these, three DL models, MLP, TabNet, and TabTransformer, were implemented to evaluate their ability to capture high-order feature interactions and complex patterns. TabNet and TabTransformer are recent architectures designed specifically for tabular data, offering attention-based mechanisms that improve interpretability and learning efficiency. The diversity of models allows for a comprehensive evaluation of predictive performance, computational cost, and interpretability across different learning paradigms.

AdaBoost (adaptive boosting) is an ensemble method that combines multiple weak classifiers, typically decision stumps, to form a strong classifier. It improves performance by focusing on the training instances that were previously misclassified. At each iteration, AdaBoost assigns a weight to the weak learner based on its accuracy and updates the weights of the training samples to emphasize difficult examples.

Given a dataset of n samples x 1 , y 1 , … , x n , y n , where y i ∈ − 1 , + 1 , the final strong classifier H x is defined as follows: H x = sign ∑ t = 1 T α t h t x where:

• T is the total number of boosting rounds,

After each iteration, the weights of the training samples are updated according to their classification outcome, thereby guiding the model to focus on harder-to-classify examples. This adaptive mechanism helps AdaBoost improve the overall prediction accuracy while maintaining good generalization.

AdaBoost has demonstrated robust performance across various classification tasks and has solid theoretical foundations, particularly in reducing both bias and variance (Schapire, 1999).

CatBoost is a gradient boosting decision tree (GBDT) algorithm developed to handle categorical features efficiently without extensive preprocessing (Dorogush, Ershov, and Gulin, 2018). The model constructs an ensemble of decision trees, where each tree T m corrects the errors of its predecessors by minimizing a differentiable loss function L : F m x = F m − 1 x + γ m T m x , where γ m is the learning rate. CatBoost incorporates ordered boosting and novel techniques to mitigate prediction shift, making it highly effective for heterogeneous tabular datasets with mixed data types.

MLP is a fundamental DL architecture composed of multiple layers of interconnected neurons, where each neuron applies an affine transformation followed by a nonlinear activation function (Goodfellow, Bengio, and Courville, 2016). Given an input feature vector x ∈ R d , the hidden layer output h is computed as follows: h = σ W 1 x + b 1 , where W 1 and b 1 are the weight matrix and bias vector, respectively, and σ · is a nonlinear activation function such as a rectified linear unit (ReLU). This hidden representation is then passed through additional layers or directly to the output layer for final prediction: y ^ = ϕ W 2 h + b 2 , where ϕ · denotes the activation function used in the output layer (e.g., softmax for classification or identity for regression). MLPs are powerful function approximators capable of modeling complex nonlinear relationships between input features and target variables.

The TabTransformer is a DL architecture tailored for tabular data, effectively modeling both categorical and numerical features through contextual embeddings and self-attention mechanisms (Huang et al., 2021). In this framework, each categorical feature x i is mapped to a dense embedding vector e i ∈ R d , forming an embedding matrix E ∈ R C × d for C categorical features: E = e 1 ; e 2 ; … ; e C

These embeddings are processed through Transformer encoder layers utilizing multi-head self-attention to capture inter-feature dependencies. The output embeddings are then concatenated with normalized numerical features x num to form a combined representation h : h = Flatten Transformer E ; x num

This representation h is subsequently passed through an MLP for downstream tasks such as classification or regression. The architecture supports self-supervised pretraining strategies, including masked feature modeling, to enhance performance on tasks with limited labeled data (Vyas and Bertsimas, 2024).

Hyperparameter optimization was critical for maximizing the predictive performance of the models. Two different strategies were applied depending on the model type: grid search for DL models and particle swarm optimization (PSO) for classical ML models.

We employed several standard evaluation metrics to comprehensively assess the performance of the developed classification models, including accuracy, precision, recall (sensitivity), F1-score, specificity, and Cohen’s kappa score. These metrics are essential for evaluating both overall and class-wise performance, especially in the presence of class imbalance (Sokolova and Lapalme, 2009; Chicco and Jurman, 2020; Kumar and Singh, 2022).

The following equations define the metrics used:

Table 4 evaluation metrics (mathematical definitions). Accuracy = T P + T N T P + T N + F P + F N  Precision = T P T P + F P  Recall = T P T P + F N F 1 ‐ Score = 2 · Precision · Recall Precision + Recall  Specificity = T N T N + F P  Kappa = p o − p e 1 − p e where p o is the observed accuracy and p e is the expected accuracy under random chance. Here, T P , T N , F P , and F N denote true positives, true negatives, false positives, and false negatives, respectively. The kappa score evaluates the degree of agreement beyond chance, where a value of κ > 0.6 generally indicates substantial agreement (McHugh, 2012).

All metrics were computed using weighted averages to accommodate class imbalance, and models were evaluated on both training and testing data to ensure generalization. In addition, receiver operating characteristic (ROC) and precision-recall (PR) curves were plotted to assess classification thresholds and model discrimination capacity (Davis and Goadrich, 2006).

The spatial distribution of the selected groundwater conditioning factors is illustrated in Figures 6–8. These factors, derived from remote sensing and geospatial datasets, were carefully selected based on their relevance to groundwater recharge dynamics in arid and semi-arid environments. Topographic features such as elevation, slope, curvature, and the TWI influence surface runoff, water accumulation, and infiltration capacity. Soil permeability and geomorphological classes reflect lithological and structural controls that govern the percolation of water through subsurface formations. The Normalized Difference Vegetation Index (NDVI) serves as an ecological proxy for vegetation cover, which is often indicative of groundwater availability in shallow aquifers. Hydrological parameters, including rainfall, lineament density, and stream distance, capture climatic inputs and structural pathways that facilitate recharge. Together, these twelve factors provide a comprehensive representation of the hydro-environmental variability across the Feija Basin and form a robust input set for data-driven groundwater potential modeling.

All performance results presented in this section are based on the optimal hyperparameter configurations obtained through the tuning procedures outlined in the Methodology section. To ensure a comprehensive evaluation, the classification performance of the five models (TabNet, TabTransformer, MLP, CatBoost, and AdaBoost) was assessed separately on the training and testing datasets. A range of standard evaluation metrics was computed, including accuracy, recall, specificity, precision, F1-score, and Cohen’s kappa. Figure 9 presents radar plots that visualize and compare the performance of each model across these metrics for both phases, providing clear insight into learning effectiveness and generalization capability.

Overall, TabNet and TabTransformer exhibited superior and consistent performance, attaining high scores across all evaluation metrics during the testing phase. AdaBoost and CatBoost also demonstrated competitive results, whereas MLP displayed moderate performance, reflecting limited generalization capacity.

Figure 10 presents the confusion matrices for each model on the test dataset. TabNet and TabTransformer achieved the most balanced predictions across the three groundwater potential (GWP) classes, misclassifying only a small number of samples. CatBoost also performed reasonably well but showed minor misclassification between Class 2 and Class 3. MLP and especially AdaBoost displayed noticeable confusion, particularly with Class 2 instances being predicted as Class 3.

To assess the discriminative power of the classifiers beyond accuracy-based metrics, we conducted a ROC analysis using the micro-average AUC approach, suitable for multi-class classification scenarios. Figure 11 presents the ROC curves for all models on both training and testing datasets.

The ensemble-based and DL models (TabNet, TabTransformer, and CatBoost) achieved near-perfect AUC values of 0.99 on both training and testing phases, reflecting excellent separability among the GWP classes. MLP and AdaBoost exhibited slightly lower AUC scores, with values of 0.96 and 0.97 on the test set, respectively. These results corroborate earlier observations from confusion matrices and classification metrics, indicating that TabNet and TabTransformer not only classify accurately but also maintain strong discriminatory performance across all classes.

The correlation matrix and mutual information analysis (not shown here) provided insights into the relationships and relevance of input features. These analyses highlighted that factors such as elevation, slope, and soil permeability have a strong influence on GWP classification. Feature importance scores obtained from each model (Figure 12) further confirmed these findings, with elevation consistently ranked as the most influential variable.

The distribution of predicted groundwater classes across models is visualized in Figure 13. Most models exhibited relatively balanced classifications among the three GWP categories. However, AdaBoost disproportionately predicted Class 3 (42.2%), indicating bias toward higher groundwater potential zones.

The spatial predictions generated by each model are presented in Figures 13, 14, 15. Models such as TabNet, TabTransformer, and CatBoost produced coherent and geographically consistent patterns that aligned well with known hydrogeological and topographic features of the study area. These models accurately delineated high-potential groundwater zones typically located in low-lying alluvial plains, valley corridors, and areas characterized by high soil permeability and vegetation density features often associated with groundwater accumulation. The spatial transitions between different GWP classes were smooth and well-localized, reflecting each model’s capacity to learn complex spatial and geophysical relationships from the input data.

From a technical standpoint, these models exhibited stronger generalization capabilities, benefiting from deep representation learning (TabNet), ensemble-based robustness (TabTransformer), and gradient-boosted refinement (CatBoost). In contrast, MLP and AdaBoost produced noisier and more fragmented spatial outputs. These inconsistencies were particularly evident in transition zones, such as the interfaces between moderate and high GWP areas, where these models often failed to capture subtle environmental gradients. AdaBoost tended to overpredict Class 3 (high potential), generating spatial overestimation and reducing practical utility for hydrogeological planning. These results emphasize the necessity of selecting models that are not only accurate in classification metrics but also capable of preserving spatial coherence and geographic relevance in the context of groundwater potential mapping.

Two robust evaluation methods were applied to provide a statistically grounded comparison of model performance: the Friedman test and the TOPSIS ranking method. The Friedman test (Figure 16, left) showed statistically significant differences among models ( χ 2 = 30.0 ; p < 0.001 ), with TabNet achieving the best average rank (1.00), followed by CatBoost (3.00), TabTransformer (4.00), and MLPr (5.00). AdaBoost ranked the lowest (6.00), confirming its underperformance.

To complement this analysis, the TOPSIS method was employed to evaluate the closeness of each model to the ideal solution across all metrics. The results (Figure 16, right) aligned with the Friedman rankings. TabNet scored highest (1.0000), indicating optimal performance, followed closely by CatBoost (0.8935) and TabTransformer (0.80). MLP and AdaBoost were assigned significantly lower scores, reinforcing their inferior metric consistency.

The aggregated evaluation from both tests is summarized in Table 5, where models were qualitatively categorized based on their combined scores.

Table 5 summarizes the integrated ranking of models based on statistical (Friedman) and multi-criteria (TOPSIS) evaluations. TabNet achieved the best overall performance, while AdaBoost was ranked lowest in both methods.