Modelling Groundwater Potential of Tropical Hard-Rock Aquifers using Spatial Cross-Validation augmented with Uncertainty and Sensitivity Analysis

Modelling Groundwater Potential of Tropical Hard-Rock Aquifers using Spatial Cross-Validation augmented with Uncertainty and Sensitivity Analysis | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Modelling Groundwater Potential of Tropical Hard-Rock Aquifers using Spatial Cross-Validation augmented with Uncertainty and Sensitivity Analysis Suraj Dule, Arabinda Sharma This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-9249309/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Sustainable groundwater management in complex hard-rock terrains requires the precise delineation of potential zones. Traditional learning approaches used for groundwater potential mapping often fail to account for the spatial dependence of hydrogeological data and the inherent uncertainty in predictions. This study presents a rigorous comparative framework for Groundwater Potential Zonation (GWPZ) evaluating six architectures: Random Forest, CatBoost, Support Vector Machines (SVM), standard XGBoost, spatially augmented XGBoost (XGB-SF), and Deep Learning Multi-Layer Perceptron (MLP-DL). To mitigate the spatial autocorrelation bias common in conventional ML/DL applications, the study proposes an innovative spatial cross-validation strategy integrated within an Optuna hyperparameter tuning workflow. Results indicated that while the XGB-SF model achieved the highest internal accuracy ( R 2 = 0.853), the standard XGBoost model proved optimal for operational mapping. It delivered the most suitable balance of predictive accuracy ( R 2 = 0.846), algorithmic stability (Bootstrap STD = 0.129), and ordinal ranking fidelity for groundwater potential (Spearman’s 𝜌 = 0.827). Furthermore, feature importance analysis revealed that tree-based ensembles successfully navigated severe multicollinearity without manual feature elimination, consistently identifying geomorphology and surface runoff as the dominant physical controls on hard-rock aquifer storage. Conversely, SVM demonstrated exceptional stability with zero overfitting, serving as a robust baseline for data-scarce regions. Crucially, this research advances decision support by generating spatially explicit uncertainty maps via quantile regression. This enables planners to distinguish high-confidence zones suitable for immediate exploitation from high-uncertainty areas requiring geophysical verification, providing a scientifically defensible framework for exploration in heterogeneous aquifer systems. Groundwater Potential Mapping Machine Learning Deep Learning Spatial Cross-Validation Spatially Explicit Uncertainty Hard-Rock Aquifers Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Figure 11 Figure 12 Figure 13 Figure 14 1. Introduction Groundwater is a critical freshwater resource for domestic, agricultural, and industrial use, serving as a primary buffer against climate variability in semi-arid regions (Famiglietti 2014). This dependence is particularly acute in India, the world's largest consumer of groundwater, where the resource supports over 60% of irrigated agriculture and 85% of rural drinking water supplies (Mukherjee et al. 2015). The resource management challenge is most severe in hard-rock terrains, which cover nearly two-thirds of the Indian peninsular shield. Unlike extensive alluvial aquifers where water is stored in continuous pore spaces, groundwater in crystalline basements (e.g., granites, gneisses) occurs discontinuously. Its storage and movement are governed strictly by secondary porosity elements, including fractures, joints, and the variable thickness of the weathered regolith (Machiwal et al. 2011 ). In these heterogeneous environments, aquifer boundaries are ill-defined and well yields can vary drastically over short distances, making the accurate delineation of Groundwater Potential Zones (GWPZ) a complex but indispensable prerequisite for artificial recharge planning. Traditional groundwater exploration relies on precise but capital-intensive hydrogeological field surveys, which are logistically prohibitive for regional assessments (Jha et al. 2010). While early geospatial efforts utilizing Remote Sensing (RS) and GIS integrated surface indicators via heuristic, knowledge-driven approaches (e.g., AHP, TOPSIS) (Dule and Sharma 2025 ; Kar et al. 2025 ; Ray 2025 ), these methods remain inherently subjective and struggle with complex, non-linear interactions (Jenifer and Jha 2017). Consequently, the field has shifted decisively toward data-driven Machine Learning (ML) algorithms, which demonstrate superior capability in modelling spatial heterogeneity (Rana et al. 2025 ; Naghibi et al. 2016 ; Saha and Pal 2025 ). A primary deficiency in current literature is treating geospatial data as independent. Most studies utilize standard random k-fold cross-validation to tune and evaluate models. This approach violates the fundamental assumption of independent and identically distributed (i.i.d.) data, as geospatial variables are invariably governed by Tobler’s First Law of Geography (Tobler 1970 ), leading to spatial autocorrelation. Random splitting leads to "spatial information leakage," where training and testing points are geographically adjacent and statistically dependent. This results in inflated, over-optimistic performance metrics that fail to generalize to unseen locations (Roberts et al. 2017). Secondly, GWPZ mapping is often mis-framed as a binary classification task or simple regression, neglecting the ordinal nature of potential zones (e.g., distinguishing 'Excellent' from 'Good') required for prioritization. Few studies validate the preservation of this ordinal rank (e.g., using Spearman’s 𝜌) in model outputs. Although recent studies have explored Deep Learning techniques such as Convolutional Neural Networks (CNNs) for spatial mapping (Panahi et al. 2020 ), the potential of modern architectures to capture geological nuances via entity embeddings remains underutilized (Guo and Berkhahn 2016 ). Most critically, standard ML outputs are deterministic, providing single predictions without quantifying epistemic uncertainty. In decision-making contexts, a prediction of 'High Potential' carries significantly different risks depending on whether the model is 90% confident or only 50% confident. The lack of spatially explicit uncertainty quantification renders many current GWPZ maps unsuitable for high-stakes resource allocation, where the cost of a dry borehole is substantial (Kosmidis et al. 2020). While recent studies have applied artificial intelligence (AI) and advanced machine learning techniques to hydrology, and others have successfully utilized optimization techniques such as Genetic Algorithms for broad regional mapping, these assessments often prioritize temporal trends over the rigorous quantification of spatial prediction uncertainty (Mussa et al. 2024 ; Dikshit et al. 2021 ; Saha et al. 2024 ). Consequently, the simultaneous integration of spatial autocorrelation correction and uncertainty analysis remains largely unexplored in hard-rock terrains (Liu et al. 2025 ; Vishwakarma et al. 2021 ). To address these gaps, this study presents a comparative machine learning framework for GWPZ mapping in the Khatra sub-division, West Bengal, driven by four methodological innovations. First, Spatial Cross-Validation was integrated with Optuna hyperparameter tuning to eliminate spatial autocorrelation bias and ensure robust model generalization. Second, model selection prioritized ordinal ranking fidelity (Spearman’s 𝜌) alongside traditional binary metrics (AUC) to ensure proper hierarchical classification of potential zones. Third, deterministic predictions were transformed into risk-aware decision-support tools by generating spatially explicit uncertainty maps via Quantile LightGBM regression. Finally, the study evaluated advanced predictive architectures, including an entity-embedded Deep Learning network (MLP-DL) and a spatially augmented XGBoost (XGBoost-SF) designed to capture multi-scale topographic context. By benchmarking these advanced architectures alongside established ensemble methods (RF, CatBoost, SVM, and standard XGBoost), this work establishes a new standard for reliable, transparent, and uncertainty-aware groundwater modelling in complex hard-rock terrains. 2. Data Sources & Methodology 2.1 Study Region The present study centres on the Khatra Sub-division of Bankura District, West Bengal, India (Fig. 1 ), encompassing a total area of approximately 2,433 km² between longitudes 86°36'52" E to 87°19'04" E and latitudes 22°37'61" N to 23°16'56" N. Administratively, the sub-division comprises eight Community Development (CD) blocks: Khatra, Indpur, Hirband, Ranibandh, Simlapal, Taldangra, Raipur, and Sarenga. The region marks the transition from the Chotanagpur Plateau to Gangetic alluvial plains. The area lies geologically on the erosional eastern fringe of the Chotanagpur Granite Gneiss Complex (CGGC), dominated by Archean crystalline rocks (granite gneiss, mica schist) with negligible primary porosity. Groundwater occurrence is strictly confined to secondary porosity elements like fractures and weather regolith (Central Ground Water Board [CGWB], 2022). Topographically, the region exhibits a westward-to-eastward gradient (400 m to 60 m Above Mean Sea Level or AMSL), featuring undulating pediments in the west that merge into gently sloping alluvial plains in the east (CGWB, 2022). The study area experiences a tropical sub-humid regime with an average annual rainfall of approximately 1,386 mm, the vast majority of which occurs during the southwest monsoon (June to September). The drainage architecture is controlled by two prominent river systems, the Kangsabati and the Shilabati. While the Kangsabati River is regulated by the Mukutmanipur Dam in the Khatra block to support canal irrigation, surface water availability remains critically insufficient for agricultural demands during the dry rabi and pre- kharif seasons (CGWB, 2022). This scarcity, combined with high surface runoff facilitated by the undulating terrain, has driven an increased dependence on groundwater extraction via borewells, placing unsustainable stress on the limited storage capacity of the fractured aquifers. 2.2 Geospatial Database and Ground Truth Inventory The spatial database for this study, anchored by a Survey of India (SOI) topographic base map (1:50,000), integrates multi-source satellite imagery, ancillary maps, and hydrological field data. Thirteen thematic predictor layers were selected based on established literature (Naghibi et al., 2016 ) and processed within a QGIS environment. Topographic and structural variables—Elevation, Slope, Drainage Density, and Lineament Density—were extracted from the ALOS PALSAR Digital Elevation Model (12.5 m resolution). Land Use/Land Cover (LULC) and the Normalized Difference Vegetation Index (NDVI) were derived from Landsat 9 OLI imagery (30 m, 2023) using the specified bands, with the LULC classification rigorously validated against ground control points via Google Earth Pro (Overall Accuracy: 90.8%, Kappa coefficient: 0.846). Geology and soil texture parameters were digitized from the Geological Survey of India (GSI) District Resource Map (1:250,000) and the National Bureau of Soil Survey and Land Use Planning (NBSS&LUP) soil texture map (1:50,000), while the geomorphological vector layer was sourced from the Bhukush GSI geoportal. Hydrometeorological forcing incorporated 2023 CHIRPS rainfall data (0.05°) and surface runoff, which was explicitly modelled using the NRCS-Curve Number (CN) method by integrating NASA HYSOG data, rainfall, and LULC. Finally, subsurface conditions, specifically Mean Water Level (2019–2023) in the form of point data was interpolated via Inverse Distance Weighting (IDW) and the Aquifer Depth layer was digitised from aquifer deposition 2013 data from the India Water Resources Information System (WRIS). To establish the target variable, groundwater yield (Litres Per Minute, LPM) was derived from NRSC ISRO prospect maps based on extensive field observations (April–May 2009 and January–February 2010). Although dynamic surface predictors utilized recent 2023 imagery, this historical baseline remains hydrologically valid because groundwater accumulation in hard-rock terrains is primarily driven by static structural controls (e.g., geomorphology, geology) rather than short-term surface changes. From this dataset, 303 points were randomly sampled and categorized into five ordinal classes (1–5 scale): Very Low (30–50), Low (50–100), Moderate (100–200), High (200–400), and Very High (400–800 LPM), then partitioned into 80% training and 20% testing subsets (random seed: 42). To prevent spatial data leakage, a completely independent external validation dataset of 180 points was curated from the same source with strictly zero spatial overlap with the primary data 2.3 Hydrogeological Rationale of Thematic Predictors The delineation of groundwater potential in hard-rock terrains depends on the non-linear interaction between surface boundary conditions and subsurface storage capacity. The thirteen selected thematic layers given in Fig. 2 and Table 1 were categorized into three functional groups: Topographic and Geomorphic Controls: Topography heavily dictates surface water retention. Lower elevations (ranging down to 1 m AMSL) facilitate water accumulation and longer residence times, favouring recharge, whereas higher elevations (up to 248 m AMSL) promote rapid runoff (Fig. 2 a) (Todd & Mays, 2004 ). Similarly, slope exhibits an inverse relationship with infiltration; gentle slopes (< 5°) enhance residence time, while steep slopes (up to 32.53°) accelerate runoff (Fig. 2 b). Geomorphologically, alluvial and flood plains offer high primary porosity, pediment complexes serve as shallow aquifers, and dissected hills act as primary runoff zones (Fig. 2 k) (Machiwal et al., 2011 ). Hydrological Drivers: Rainfall (1168–1360 mm) acts as the primary source of recharge, directly correlating with increased potential in permeable zones (Fig. 2 c) (Adiat et al., 2012 ). Conversely, surface runoff (1047–1351 mm) and high drainage density (up to 3.14 km/km²) serve as inverse indicators; high values signify impermeable lithology and low infiltration capacity (Fig. 2 d, e) (Chowdary et al., 2009 ; Yeh et al., 2016 ). Furthermore, Land Use/Land Cover (LULC) and NDVI (-0.75 to 0.88) dictate percolation efficiency; vegetated and agricultural lands reduce runoff velocity and increase soil porosity via root systems, whereas impervious built-up areas impede recharge (Fig. 2 f, j) (Sander, 2007 ; Shaban et al., 2006 ). Subsurface and Structural Constraints: In crystalline Geological units like Granite Gneiss and Schist, groundwater movement relies entirely on secondary porosity, making high Lineament Density (up to 0.56 km/km²) critical conduits for flow (Fig. 2 g, m) (Mogaji et al., 2015 ; Singhal & Gupta, 2010 ). Pedologically, coarse Soil Textures like sandy loam facilitate rapid recharge compared to clay-rich soils (Fig. 2 l) (Hillel, 2013 ). Subsurface capacity is indicated by Mean Water Level and Aquifer Depth; shallow water levels (< 5 mbgl) signify saturated valley fills, while thicker aquifers (90–150 m) indicate larger storage reservoirs capable of sustaining extraction during dry periods (Fig. 2 h, i) (Prasad et al., 2008 ; Todd & Mays, 2004 ). Figure 2 Spatial distribution of input thematic layers: (a) Elevation, (b) Slope, (c) Rainfall, (d) Runoff, (e) Drainage Density, (f) NDVI, (g) Lineament Density, (h) Mean Water Level, and (a) Aquifer Depth, (j) LULC, (k) Geomorphology, (l) Soil Texture, and (m) Geology. Table 1 Data sources and characteristics of the thematic layers. Parameters Source Data Type Native Resolution/Scale Elevation, Slope, Drainage Density, Lineament Density ALOS PALSAR from Alaska Satellite Facility (ASF) and Japan Aerospace Exploration Agency (JAXA). ( https://datapool.asf.alaska.edu ) Raster 12.5 m Rainfall Climate Hazards Group Infrared Precipitation with Stations (CHIRPS). ( https://data.chc.ucsb.edu ) Raster 0.05° x 0.05° Runoff Obtained using NRCS-CN Method Raster 12.5 m NDVI, LULC Landsat 9 OLI (30 m). ( https://earthexplorer.usgs.gov ) Raster 30 m Geomorphology, Geology Geological Survey of India (GSI) (1:250k). ( https://gsi.gov.in/home ) Vector and Raster 1:250k Soil Texture National Bureau of Soil Survey & Land Use Planning (NBSS&LUP). ( https://icar-nbsslup.org.in ) Vector 1:50k Mean Water Level, Aquifer Depth India Water Resource Information System (WRIS). ( https://indiawris.gov.in ) Vector and Raster Points, 1:250k Groundwater Prospect The National Remote Sensing Centre (NRSC) and ISRO. ( https://maps.wbphed.gov.in ) Raster 1:50k 2.4 Data Preprocessing and Algorithm Architectures To ensure spatial consistency across the multi-source database (Fig. 3 ), all thematic layers were resampled to a unified 12.5 m resolution, using the ALOS PALSAR DEM as the reference master layer. Nearest-neighbour interpolation was employed for categorical layers to preserve discrete class fidelity, while bilinear interpolation was applied to continuous variables to maintain surface gradients in QGIS software. Following the data resampling and interpolation step, six distinct architectures were evaluated. Preprocessing pipelines and feature engineering were tailored specifically to the mathematical requirements of each algorithm. 2.4.1 Random Forest (RF) Random Forest, developed by Leo Breiman ( 2001 ), is an ensemble learning method that constructs a multitude of decision trees at training time. It operates on the principle of bagging (bootstrap aggregating), where each tree is trained on a random subset of the data sampled with replacement. Node splitting within the individual trees is governed by impurity metrics such as the Gini index (Eq. 1 ): $$\:Gini=1-{\sum\:}_{i=1}^{C}{\left({p}_{i}\right)}^{2}$$ 1 Where C is the number of classes and p i is the proportion of class i . RF introduces a second layer of randomness by selecting a random subset of features at each split point (controlled by max_features). This decorrelates the trees, significantly reducing the variance of the model without increasing bias, making it highly robust to noise and overfitting. For this, in RF and the other tree-based ensembles, Min-Max scaling [0, 1] is used for continuous features, while categorical variables are integer-encoded. Because the model is dealing with a continuous regression task, the algorithm must aggregate the independent outputs of the forest into a single continuous value. It achieves this by calculating the arithmetic mean of the predictions from all N individual trees ( f i ), yielding the final groundwater potential prediction ( F ( x )) for any given point ( x ), using Eq. 2 : $$\:F\left(x\right)=\frac{1}{N}{\sum\:}_{i=1}^{N}{f}_{i}\left(x\right)$$ 2 2.4.2 Extreme Gradient Boosting (Standard XGBoost) Unlike RF, which builds independent trees in parallel, Gradient Boosting builds an ensemble sequentially. Each new tree is trained to predict the pseudo-residuals (errors) of the preceding ensemble, effectively converting weak learners into a strong learner (Friedman, 2001 ). eXtreme Gradient Boosting (XGBoost) is an advanced implementation optimized for speed and performance (Chen, 2016 ), minimizing a regularized objective function ( Obj ) that combines a convex loss function ( L ) with a penalty term (𝛺) for model complexity (Eq. 3 ): $$\:Obj={\sum\:}_{i=1}^{n}L\left({y}_{i},\widehat{{y}_{i}}\right)+{\sum\:}_{m=1}^{M}{\Omega\:}\left({f}_{m}\right)$$ 3 where n represents the total number of training instances, \(\:{\text{y}}_{\text{i}}\) is the actual observed groundwater yield, and \(\:\widehat{{\text{y}}_{\text{i}}}\) is the model's prediction. The first term, \(\:\text{L}\left({\text{y}}_{\text{i}},\widehat{{\text{y}}_{\text{i}}}\right)\) , calculates the training loss by measuring the difference between the actual and predicted values. The second term introduces the regularization penalty \(\:{\Omega\:}\) applied across all M sequential trees \(\:{f}_{m}\) in the ensemble. Specifically, the penalty is defined as \(\:{\Omega\:}\left({f}_{m}\right)={\gamma\:}\text{T}+\frac{1}{2}{\lambda\:}{\left|\left|\text{w}\right|\right|}^{2}\) , where T represents the total number of leaves in the tree, w represents the leaf weights (scores), and \(\:{\gamma\:}\) and \(\:{\lambda\:}\) are parameters that control the severity of the regularization. This explicit penalization is critical for preventing the model from overfitting on limited, high-dimensional training data. This study configured the model with tree_method='hist' and enable_categorical=True to natively handle the integer-coded categorical features with efficiency. 2.4.3 Spatially-Augmented XGBoost Standard machine learning algorithms remain inherently "point-based," often failing to account for spatial autocorrelation as described by Tobler’s First Law of Geography (Tobler, 1970 ). To address this limitation, a spatially-augmented variant (XGBoost-SF) was developed. It shares the identical mathematical foundation as Standard XGBoost but incorporates engineered spatial textures. Three features were derived from the Elevation raster: a Focal Mean (5 x 5-pixel window) to smooth noise and highlight regional trends; a Focal Standard Deviation to quantify local topographic roughness; and a Sobel Gradient (3 x 3 filter) to assess multi-directional slope steepness. By expanding the input matrix with these textures, the gradient boosting objective function explicitly learns regional topographic transitions rather than treating each pixel as an isolated entity. 2.4.4 CatBoost Categorical Boosting (CatBoost) is a GBDT algorithm specifically engineered to handle categorical data (Prokhorenkova et al., 2018 ). Its primary innovation is "ordered target statistics"- a permutation-based strategy for encoding categorical features that prevents target leakage and allows for the full, unbiased utilization of high-cardinality discrete features. It also grows "symmetric" or "oblivious" trees, which act as a form of regularization and significantly speed up inference. This architectural symmetry was deemed highly relevant given the established importance of categorical layers like Geomorphology and Geology in hydrogeological modelling. 2.4.5 Support Vector Machine (SVM) Support Vector Regression (SVR) is a kernel-based method rooted in statistical learning theory (Vapnik, 1995). Unlike traditional least-squares regression, SVR seeks a flat function f ( x ) deviating from the target y i by at most a margin 𝝐 (Drucker et al., 1997). This 𝝐-insensitive loss function only penalizes errors falling outside this specific margin. To model non-linear spatial relationships, data is mapped to a higher-dimensional feature space using the Radial Basis Function (RBF) kernel (Eq. 4 ): $$\:\left({x}_{i},{x}_{j}\right)=\text{exp}\left(-\gamma\:{\left|\left|{x}_{i}-{x}_{j}\right|\right|}^{2}\right)$$ 4 where \(\:{x}_{i}\:\) and \(\:{x}_{j}\) are distinct geographical feature vectors, and \(\:{\left|\left|{x}_{i}-{x}_{j}\right|\right|}^{2}\) calculates the squared Euclidean distance. The kernel coefficient \(\:{\gamma\:}\:\) (gamma) defines the similarity radius of individual support vectors. Together with the regularization penalty \(\:C\) (for margin violations) and 𝝐, it dictates the model's complexity and ability to generalize. Unlike tree ensembles, SVM strictly requires One-Hot Encoding for categorical inputs and Standard Scaling (zero mean, unit variance) to prevent magnitude bias in distance calculations (Hastie et al., 2009 ; Géron, 2022 ). Due to its \(\:\mathcal{O}\left({\text{n}}^{3}\right)\) computational complexity, the SVR was trained on a random sub-sample of 8,000 points. 2.4.6 Multi-Layer Perceptron Deep Learning (MLP-DL) A custom tabular Deep Learning (DL) model was implemented in PyTorch to process mixed data types. To prevent the sparsity of one-hot encoding, high-cardinality categorical variables (e.g., Geology) were passed through nn.Embedding layers (Guo & Berkhahn, 2016 ). This generated dense, low-dimensional vector representations capturing semantic class similarities. These embeddings were concatenated with continuous inputs and routed through a simplified Self-Attention mechanism (Vaswani et al., 2017 ) to dynamically weigh feature importance. The core backbone consisted of a Residual MLP with four hidden layers ([256, 256, 128, 64]) utilizing LeakyReLU activation (Eq. 5 ): $$\:f\left(x\right)={max}\left(0.01x,x\right)$$ 5 where x represents the input signal to the neuron. Unlike standard ReLU, LeakyReLU applies a small positive slope (0.01) to negative values, ensuring continuous gradient flow and mitigating the vanishing gradient problem. Finally, the activated data passed through Dropout (0.1) layers for regularization. 2.5 Hyperparameter Optimization and Spatial Cross-Validation Except for the MLP-DL, all architectures were optimized via Optuna using a 50-trial Tree-Structured Parzen Estimator (TPE) (Akiba et al., 2019 ; Online Resource 1). To prevent the optimistic bias and spatial data leakage typical of random sampling, Spatial Cross-Validation (SpatialKFold) was integrated directly into the objective function. By partitioning the study area into geometrically disjoint blocks (2x2 geographic grid) via coordinate k-means clustering (Fig. 4 ), models were forced to learn generalized hydrogeological rules rather than memorizing local coordinates. This realistically simulates exploration in uncharted terrains. Additionally, the RF and SVM objective functions incorporated a penalty for train-validation divergence to explicitly enforce algorithmic generalization. Figure 4 Spatial Cross-Validation blocks (coloured clusters). 2.6 Feature Diagnostics To ensure the model relies on physically valid hydro-climatic relationships rather than spurious correlations, we implemented a dual-stage diagnostic framework. The first stage focused on multicollinearity screening, specifically restricted to continuous climatic and topographic predictors (e.g., Rainfall, Slope, Elevation). Variance Inflation Factor (VIF) analysis, a linear diagnostic derived from Ordinary Least Squares (OLS) regression, was utilized for this purpose. Applying VIF to high-cardinality categorical features like Geology would require one-hot encoding, which introduces artificial multicollinearity (the 'dummy variable trap') without providing meaningful insight. By restricting VIF solely to continuous raster, physical redundancies—such as those between Elevation and Slope—were successfully identified without introducing numerical artefacts. The second stage evaluated global feature importance to assess the predictive power of the full mixed-type feature set. For this, model-agnostic techniques, specifically SHAP (SHapley Additive exPlanations) and Permutation Importance, were employed. Unlike VIF, these methods quantify the marginal contribution of a feature to the model's decision boundary rather than its linear correlation. This approach allows for a direct, unbiased comparison between continuous drivers (e.g., NDVI) and categorical controls (e.g., Geomorphology) on a unified scale, effectively capturing both linear and non-linear interactions that traditional diagnostics typically miss (Lundberg & Lee, 2017 ; Fisher et al., 2019 ). 2.7 Model Interpretability A multi-stage validation workflow was implemented to assess model robustness, stability, and feature logic beyond simple accuracy metrics. To establish baseline predictive performance, internal validation on the 20% hold-out set utilized the Coefficient of Determination ( R 2 ), Root Mean Square Error (RMSE), and Mean Absolute Error (MAE). Generalization was then rigorously tested against an independent external ground truth dataset. This external phase employed the Area Under the Receiver Operating Characteristic Curve (AUC-ROC) to evaluate binary discrimination capability, alongside Spearman's Rank Correlation (𝜌) to assess the models' ordinal fidelity across the discrete GWPZ classes. Beyond standard error metrics, the workflow explicitly quantified algorithmic uncertainty and input sensitivity. Model stability was evaluated via a 60-iteration bootstrap analysis (Efron & Tibshirani, 1994 ), which calculated the mean per-point standard deviation to establish the deterministic consensus of each architecture. Finally, non-linear relationships were mathematically decomposed using Sobol Global Sensitivity Analysis (Soboĺ, 1993 ; Saltelli et al., 2010 ). By calculating First-Order ( \(\:{S}_{i}\) ) and Total-Order ( \(\:{S}_{T}\) ) indices for the continuous features, this analysis successfully isolated the variance contribution attributable strictly to individual inputs versus complex multi-feature interactions. 2.8 Spatial Autocorrelation Analysis To evaluate spatial independence—a fundamental assumption in geospatial modelling—autocorrelation was assessed using the esda package (Rey & Anselin, 2007 ). A spatial weights matrix was defined via k -Nearest Neighbours \(\:(k=5)\) with significance tested across 999 simulations. The analysis was two-fold: (a) Input Data: Global Moran's I (Moran, 1950 ) was calculated for the raw GWPZ values. A significant clustering \(\:(p<0.05)\) confirmed the necessity of the Spatial Cross-Validation strategy. (b) Residual Analysis: Post-modelling, both Global Moran's I and Geary's C (Geary, 1954 ) were computed on the test set residuals ( e i ). A non-significant result \(\:(p>0.05)\) was required to confirm that the model successfully captured the underlying spatial structure without leaving systematic, clustered errors. 2.9 Spatial Uncertainty Analysis To move beyond a single mean prediction and quantify epistemic (model) uncertainty, a quantile regression workflow was implemented. Three LightGBM models (Ke et al., 2017 ) were trained using the quantile objective to predict the 5th, 50th, and 95th percentiles (Koenker & Bassett Jr, 1978 ). LightGBM was selected for its computational efficiency in estimating conditional quantiles. The final uncertainty map was generated by calculating the 90% Prediction Interval (PI) width \(\:({\mathcal{Q}}_{95}-{\mathcal{Q}}_{05})\) for each pixel, allowing stakeholders to distinguish high-confidence zones from areas requiring further field verification spatially. We extracted the 95% Prediction Intervals (PI) for the independent test set to assess the coverage probability, comparing the predicted mean yield against the width of the uncertainty bound. 3. Results and Discussion This study rigorously evaluated six distinct machine learning architectures for Groundwater Potential Zonation (GWPZ): RF, CatBoost, SVM, MLP-DL, Standard XGBoost, and XGBoost-SF. All models underwent identical Optuna hyperparameter tuning and spatial cross-validation (CV) to mitigate spatial autocorrelation. The assessment covers spatial analysis, internal predictive mechanics, and definitive external validation. 3.1 Internal Predictive Mechanisms Model predictive capability was initially evaluated on the 20% hold-out test set to assess the trade-off between pure accuracy and algorithmic generalization. The spatially-augmented XGBoost-SF demonstrated superior internal accuracy, yielding the highest Test R 2 (0.853) and lowest errors (RMSE = 0.508, MAE = 0.403) with a minimal overfitting gap of 1.1% compared to its Training R 2 (0.864). This outperformed the standard XGBoost (Test R 2 = 0.846, Train R 2 = 0.876, gap = 3.0%, RMSE = 0.520, MAE = 0.425), confirming that the integration of the engineered focal mean elevation feature provided highly valuable spatial context. Visualized via a radar chart (Fig. 4 ), XGBoost-SF achieved the optimal balance of high accuracy and low error. The remaining tree-based ensembles also exhibited strong predictive capabilities and small generalization gaps, indicating that regularization parameters optimized via Optuna successfully prevented memorization. CatBoost achieved a Test R 2 of 0.834 (Train R 2 = 0.863, gap = 2.9%, RMSE = 0.540, MAE = 0.419), followed closely by Random Forest with a Test R 2 of 0.823 (Train R 2 = 0.865, gap = 4.2%, RMSE = 0.558, MAE = 0.432). Conversely, the MLP-DL and SVM architectures displayed highly constrained learning, avoiding overfitting almost entirely. The MLP-DL achieved a Test R 2 of 0.831 (RMSE = 0.544, MAE = 0.423) with a negligible overfitting gap of 0.5% against its Training R 2 (0.836), suggesting the Dropout layers and entity embeddings allowed the network to learn a highly generalized representation. Remarkably, SVM displayed a negative overfitting gap (-2.0%), performing better on the test set (Test R 2 = 0.786, RMSE = 0.613, MAE = 0.521) than the training sub-sample (Train R 2 = 0.766). While this testifies to the powerful global regularization of the SVM’s objective function—establishing it as a robust baseline for data-scarce regions—it resulted in lower overall peak accuracy compared to the boosting models. 3.2 Hydrogeological Drivers and Feature Logic Feature importance was triangulated using Variance Inflation Factor (VIF), SHAP (SHapley Additive exPlanations), Permutation Importance, and Sobol Sensitivity Analysis to ensure physical validity and algorithm transparency. 3.2.1 Multicollinearity Test VIF analysis (Fig. 5 ; Online Resource 1) confirmed the physical independence of the base thematic predictors, with all continuous variables—including Runoff (VIF = 5.25) and Rainfall (VIF = 4.56)—registering well below the severe multicollinearity threshold of 10. Conversely, the introduction of engineered spatial textures in the XGBoost-SF model inherently generated extreme mathematical redundancy, specifically between Slope and GWPZ Gradient (VIF > 71,000) and Elevation and GWPZ Mean (VIF > 8,000). 3.2.2 Feature Importance and Sensitivity Synthesizing global SHAP values (Fig. 6 ) and Permutation Mean R 2 Drop (Fig. 7 ) reveals how algorithms navigated collinearity and classified predictors into hydrogeological tiers. First, tree-based models successfully managed dataset redundancies without manual feature elimination by consistently prioritizing Runoff over Rainfall (mean SHAP < 0.012) (Breiman, 2001 ; Chen, 2016 ). This is physically sound, as Runoff acts as a localized synthesis of soil texture, land use, and slope, whereas rainfall provides a uniform regional input. Beyond resolving collinearity, a dominant triad of Geomorphology, LULC, and NDVI emerged universally. Geomorphology acted as the absolute primary control (SHAP: 0.681–0.804; Mean R 2 Drop > 0.95), dictating static aquifer boundary conditions, while LULC and NDVI served as critical proxies for dynamic surface permeability. Because Geomorphology dominates this discrete space, Sobol Sensitivity Analysis (Fig. 8 ) was utilized to isolate the variance from continuous physical variables. The resulting first-order indices ( S 1 ) identified NDVI as the primary contributor ( S 1 = 0.6667), followed by Drainage Density ( S 1 = 0.0962) and Elevation ( S 1 = 0.0517). Thus, while Geomorphology sets the structural container, vegetation health (NDVI) fine-tunes near-surface soil moisture and shallow groundwater availability during dry seasons. Furthermore, analysis of the spatially augmented XGBoost-SF (Fig. 6 b) validated the spatial feature engineering methodology. The focal mean texture feature (GWPZ_MEAN) emerged as the fourth most important predictor (SHAP = 0.077), demonstrating that providing the gradient boosting algorithm with 5x5-pixel regional topographic awareness significantly influences its decision tree splits. Non-tree algorithms exhibited unique inductive biases. The MLP-DL model (Fig. 6 f) elevated Geology (SHAP = 0.150) and Soil Texture to much greater prominence than the tree-based models, where Geology caused an anomalous 0.094 drop in R 2 . This is directly attributed to the nn.Embedding layers (Guo & Berkhahn, 2016 ), which allowed the neural network to learn dense, multi-dimensional vector representations and exploit complex semantic relationships that are missed by standard integer-encoding. Conversely, the SVM demonstrated a distinct kernel bias towards hydrologic flow. Because SVM required one-hot encoding, the categorical dominance of Geomorphology was fragmented across multiple dummy variables (e.g., Geomorphology_7, Geomorphology_1). This fragmentation allowed Runoff to become a highly influential feature (SHAP = 0.205; Permutation R 2 Drop = 0.112), demonstrating the RBF kernel's sensitivity to continuous, scaled hydrologic gradients. 3.3 Spatial Zonation Analysis To evaluate how the statistical differences identified in the previous section translate into physical zoning, the continuous Groundwater Potential Index (GWPI) generated by all six optimized models was classified using the Natural Breaks (Jenks) method to minimize intra-class variance (Jenks, 1967 ). This aligned with the ordinal-rank validation strategy, classifying the study area into five ordinal zones: "Very Low," "Low," "Moderate," "High," and "Very High." Visual inspection of the classified maps (Fig. 9 ) alongside their comparative areal distributions (Fig. 10 ) reveals critical insights into model stability, geomorphic adherence, and algorithmic inductive bias. XGBoost (Standard & SF): Both the Standard XGBoost (Fig. 9 a) and the XGBoost-SF (Fig. 9 b) produced hydro-geologically plausible maps with remarkably balanced areal distributions. The combined 'Very High' and ‘High’ zones (4.9% in Standard, 5.2% in SF) were clearly delineated along primary river valleys characterized by gentle slopes and saturated sediments, creating smooth, logical transitions to 'Moderate' weathered pediments. The ‘Moderate’ zone constitutes approximately 34.2% (Standard) to 36.3% (SF) of the area. Hydrologically, this aligns perfectly with the reality of hard-rock terrains, where true high-yield aquifers are rare, but a significant portion of the landscape consists of weathered mantles offering moderate, seasonally variable yields. The remaining area falls under the 'Low' (31% in Standard, 29.8% in SF) and 'Very Low' (29.9% in Standard, 28.7% in SF) zones, which are heavily concentrated in the south-western and south-to-southeastern regions. Notably, the XGB-SF map exhibited slightly expanded 'Moderate' zones in upland plateaus due to the engineered focal mean feature. CatBoost: This model generated a highly conservative landscape dominated by ‘Low’ (63.6%) and ‘Very Low’ (26.5%) zones (Fig. 9 c). High potential was strictly confined to narrow river corridors, with ‘High’ and ‘Very High’ zones representing a negligible 5% combined. The ‘Moderate’ potential zone was exceptionally limited (4.9%), restricted to areas in proximity to the southern river sections. MLP-DL: Conversely, MLP-DL demonstrated the most optimistic areal distribution, shifting the spectral curve toward higher potentiality (Fig. 9 f). It identified a combined 27.2% of the region as ‘High’ to ‘Very High’ potential and allocated the largest proportion of any model (40.8%) to the ‘Moderate’ class. The remaining elevated forested land was classified as ‘Low’ (15.9%) and ‘Very Low’ (16%). This tendency to over-predict in marginal areas and produce fragmented outputs reflects the stochastic training nature of neural networks on limited datasets, which often fail to suppress noisy signals. RF: The RF model allocated 5% and 12.3% of the study area to the ‘Very High’ and ‘High’ zones, respectively (Fig. 9 d). While these were generally situated close to waterbodies, the model anomalously predicted high potential across the south-western hilly tract. Almost the entire northern section was identified as ‘Moderate’ (20.9%), whereas a vast area covering the central, eastern, and select western and southern regions was classified as ‘Low’ (38.6%) and ‘Very Low’ (23.2%). SVM: Displaying large, contiguous blocks of potential, the SVM map (Fig. 9 e) reflects the "gap penalty" and RBF kernel properties of the algorithm (Vapnik, 2013 ). While this generalization creates broad zoning, it glosses over the local hydrogeological nuances required for site-specific exploration. It predicted 4.6% and 8.5% of the area as ‘Very High’ and ‘High’ potential, primarily within southern, western, and eastern river valleys. Adjoining areas were classified as ‘Moderate’ (28.7%), while rocky uplands and elevated hills defaulted to the ‘Low’ (34.2%) and ‘Very Low’ (24.0%) potential zones. 3.4 Spatial Robustness and Uncertainty The integration of spatial diagnostics and bootstrap resampling allowed for a rigorous evaluation of model reliability beyond standard error metrics. 3.4.1 Spatial Independence To confirm the mitigation of spatial autocorrelation bias, the distribution of spatial cross-validation blocks was visualised (Fig. 4 ). Assigning ground truth points to geographically distinct 2 x 2 blocks ensures that the training and validation subsets are physically separated by a sufficient distance to minimise information leakage. As indicated by the non-significant Global Moran’s I ( p = 0.218) and Geary’s C ( p = 0.106) with values 0.036 and 0.899 for the test set residuals, this blocking strategy successfully compelled the model to learn regional hydrogeological drivers rather than local spatial dependencies. Consequently, the reported predictive accuracy is a robust measure of the model's ability to generalize to unseen terrains within the crystalline basement of the Khatra region. 3.4.2 Model Stability and Confidence The deterministic consensus of the optimal models was interrogated through a 60-iteration bootstrap analysis (Efron & Tibshirani, 1994 ). A clear trade-off emerged between predictive performance and algorithmic stability. Regarding algorithmic mechanics, SVM proved the most deterministic (Bootstrap Uncertainty STD = 0.119), a result attributed to its global convex optimization which guarantees a unique solution (Vapnik, 2013 ). This was closely followed by the highly stable tree-based ensembles: Standard XGBoost (STD = 0.129), Random Forest (STD = 0.135), and the spatially-augmented XGBoost-SF (STD = 0.141). CatBoost demonstrated moderate stability (STD = 0.197). In contrast, the MLP-DL exhibited the highest instability (STD = 0.340) due to its non-convex loss landscape (Paszke et al., 2019 ), manifesting as high-frequency noise and stochastic variance in the spatial outputs. This suggests that for smaller geospatial datasets (< 10,000 points), the complexity of deep learning may introduce unnecessary stochasticity compared to gradient boosting or kernel methods. Ultimately, this algorithmic stability metric serves as a crucial prerequisite; however, the final model selection must also satisfy rigorous independent testing discussed underneath. 3.5 Model Validation and Selection To evaluate true hydrogeological applicability, the six models were tested against an independent external dataset using AUC-ROC for binary discrimination and Spearman's Rank Correlation (𝜌) for ordinal fidelity. While all models exhibited exceptional binary discrimination (Fig. 11 )—with CatBoost achieving the highest AUC (0.9997), followed by MLP-DL (0.9961), XGBoost-SF (0.9919), SVM (0.9919), Standard XGBoost (0.9892), and RF (0.9856)—relying solely on AUC for five-class zonation is misleading. It often fails to capture subtle, transitional 'Moderate' zones critical for regional planning, excelling primarily at separating absolute extremes ('Very High' vs. 'Very Low'). Conversely, Spearman’s Rho (𝜌) revealed the true predictive hierarchy (Fig. 12 ). Standard XGBoost achieved the highest rank fidelity (𝜌 = 0.8272), followed closely by CatBoost (𝜌 = 0.8244) and XGBoost-SF (𝜌 = 0.8194). MLP-DL (𝜌 = 0.7596) and SVM (𝜌 = 0.7071) showed moderate performance, while Random Forest recorded a significant drop (𝜌 = 0.6014). CatBoost’s near-perfect AUC but lower 𝜌 reflects a binary nature that struggles with intermediate classes, whereas Standard XGBoost accurately captured continuous hydrogeological gradients. Consequently, Standard XGBoost was selected as the optimal architecture. It delivered the highest ordinal fidelity, coherent spatial logic, and high bootstrap stability (STD = 0.129), all without the computational overhead of the XGBoost-SF variant. With Standard XGBoost established for definitive point-estimates, the final critical step is quantifying the model's spatial confidence boundaries, evaluated in the upcoming section. 3.6 Spatial Uncertainty Mapping While Standard XGBoost provided the optimal deterministic GWPZ scores, the resulting 90% Prediction Interval (PI) map (Fig. 13 ) fundamentally alters practical decision support by visually representing the model's confidence levels. Unlike standard models that provide a single baseline estimate, this approach maps the variance in predictions. The output reveals narrow PI zones that indicate high model consensus, identifying safe, low-risk targets suitable for immediate investment. Conversely, the wide PI zones highlight regions of high epistemic uncertainty—areas where the model exhibits lower confidence due to conflicting features or sparse data. By explicitly mapping these confidence boundaries, stakeholders are cautioned against blind drilling; instead, these high-uncertainty areas can be designated for mandatory geophysical verification, effectively bridging the gap between predictive modelling and field-level risk management. 3.7 Limitations and Future Scope While the proposed Uncertainty-Quantified framework demonstrates high algorithmic stability and successfully mitigates spatial autocorrelation, a minor epistemic uncertainty is introduced by the decadal temporal gap between the static ground truth yield points (2009–2010) and the dynamic surface predictors (2023 LULC/NDVI). Consequently, the primary contribution of this research is the demonstration of a spatially robust machine learning pipeline—integrating Spatial CV with Quantile Regression—rather than a contemporary operational audit of the Khatra sub-division. Future studies will focus on deploying this validated ML architecture across real-time, high-density piezometric monitoring networks to further operationalize and refine its predictive capabilities. 4. Conclusion This study established a rigorous comparative framework to evaluate the trade-offs between predictive accuracy, stability, and hydrogeological interpretability for Groundwater Potential Zonation (GWPZ) in hard-rock terrains. Moving beyond single-algorithm applications, the Standard XGBoost model emerged as the superior architecture for definitive spatial mapping. While CatBoost excelled in binary discrimination and XGBoost-SF in internal precision, Standard XGBoost achieved the optimal synthesis of performance and ordinal rank fidelity (𝜌 = 0.827). Methodologically, this study demonstrates that standard random validation is insufficient for heterogeneous geospatial data. Implementing Spatial Cross-Validation revealed that conventional metrics inherently overestimate model reliability. Enforcing spatial independence provides a conservative, physically grounded benchmark reflecting true exploration challenges in uncharted aquifers. Regarding algorithmic dynamics, SVMs proved to be robust baselines for conservative planning, whereas MLP-DL exhibited high stochastic variance, indicating a need for further regularization in data-scarce domains. Physically, all models consistently prioritized Geomorphology and Runoff, confirming that landform geometry and surface permeability remain the primary controls on local aquifer storage. Crucially, this work advances the paradigm from deterministic mapping to risk-aware spatial planning via a novel dual-algorithm framework. By coupling Standard XGBoost for primary GWPZ estimation with Quantile LightGBM regression for spatial uncertainty mapping, this study provides stakeholders with quantifiable confidence bounds. This clearly delineates high-confidence valley fills from high-uncertainty upland transitions, offering actionable guidance for regional planners to prioritize development in 'High Potential/Low Uncertainty' zones while mandating geophysical investigations in 'High Uncertainty' areas. Furthermore, the strong influence of land-use features suggests moderate-potential pediment plains are optimal for Managed Aquifer Recharge (MAR) interventions. Future research should integrate dynamic temporal variables, such as seasonal water level fluctuations, to refine operational yield estimates. Ultimately, this study facilitates a vital shift from static 'map-making' to dynamic 'risk-management,' empowering water managers to reject risky high-yield predictions in favour of stable, verified alternatives. Declarations Ethics approval and consent to participate: Not applicable. This study does not involve human participants, human tissue, or animal subjects. Consent for publication: Not applicable. The manuscript does not contain data from any person. Competing Interests: The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Funding: The first author (Suraj Dule) acknowledges the financial support provided by the University Grants Commission (UGC), Government of India, in the form of a Junior Research Fellowship (Ref. No. 210510010833). Author Contribution Suraj Dule: Methodology, Validation, Resources, Software, Formal analysis, Investigation, Data curation, Writing – original draft, Visualization. Arabinda Sharma: Conceptualization, Methodology, Resources, Supervision, Writing – review & editing. All authors read and approved the final manuscript. Acknowledgement The first author acknowledges the UGC for providing financial assistance through the Junior Research Fellowship. The authors gratefully acknowledge the School of Geography, Gangadhar Meher University, for providing the necessary computational infrastructure. We also thank the Geological Survey of India (GSI) and Central Ground Water Board (CGWB) for making the ancillary geospatial data publicly available. Data Availability The raw spatial, geological, and hydro-meteorological datasets analysed in this study are publicly available from their respective institutional repositories (e.g., Geological Survey of India, USGS EarthExplorer, Climate Hazards Center, and India WRIS) as detailed in Table 1. The processed geospatial matrices, extracted validation point data, and generated machine learning model artifacts supporting the findings of this study are available from the corresponding author upon reasonable request. References Adiat, K. A. N., Nawawi, M. N. M., & Abdullah, K. (2012). Assessing the accuracy of GIS-based elementary multi criteria decision analysis as a spatial prediction tool–a case of predicting potential zones of sustainable groundwater resources. Journal of hydrology , 440 , 75–89. Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. (2019, July). Optuna: A next-generation hyperparameter optimization framework. In Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining (pp. 2623–2631). Breiman, L. (2001). Random forests. Machine learning , 45 (1), 5–32. Brenning, A. (2012, July). Spatial cross-validation and bootstrap for the assessment of prediction rules in remote sensing: The R package sperrorest. In 2012 IEEE international geoscience and remote sensing symposium (pp. 5372–5375). IEEE. Central Ground Water Board. (2022). Report on national aquifer mapping and management plan of Bankura district, West Bengal . Department of Water Resources, River Development and Ganga Rejuvenation, Government of India. https://cgwb.gov.in/cgwbpnm/download/897 Chen, T. (2016). XGBoost: A Scalable Tree Boosting System. Cornell University . Chowdary, V. M., Chakraborthy, D., & Jeyaram, A. (2009). Multi-criteria decision-making approach for watershed management using remote sensing and GIS. Water Resources Management , 23(8), 1553–1573. Dikshit, A., Pradhan, B., & Alamri, A. M. (2021). Pathways and challenges of the application of artificial intelligence to geohazards modelling. Gondwana Research , 100 , 290–301. Drucker, H., Burges, C. J., Kaufman, L., Smola, A., & Vapnik, V. (1996). Support vector regression machines. Advances in neural information processing systems , 9 . Dule, S., & Sharma, A. (2025). Deciphering groundwater potential zones of Sarenga block, West Bengal, India using Geographic Information System and Multi-Criteria Decision-Analysis methods. International Journal of Engineering and Geosciences , 11 (2), 301–320. Efron, B., & Tibshirani, R. J. (1994). An introduction to the bootstrap . Chapman and Hall/CRC. Fisher, A., Rudin, C., & Dominici, F. (2019). All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously. Journal of Machine Learning Research , 20 (177), 1–81. Friedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. Annals of statistics , 1189–1232. Geary, R. C. (1954). The contiguity ratio and statistical mapping. The incorporated statistician , 5 (3), 115–146. Géron, A. (2022). Hands-on machine learning with Scikit-Learn, Keras, and TensorFlow . " O'Reilly Media, Inc.". Guo, C., & Berkhahn, F. (2016). Entity embeddings of categorical variables. arXiv preprint arXiv:1604.06737 . Hastie, T., Tibshirani, R., & Friedman, J. (2009). The elements of statistical learning: Data mining, inference, and prediction (2nd ed.). Springer. Hillel, D. (2013). Introduction to soil physics . Academic press. Jenks, G. F. (1967). The Data Model Concept in Statistical Mapping. International Yearbook of Cartography , 7 , 186–190. Kar, S., Barman, B., & Das, J. (2025). Identification of Groundwater Potential Zones in North-Western Bankura, West Bengal, India, Using GIS Based MCDA Techniques TOPSIS and CoCoSo. Journal of the Indian Society of Remote Sensing , 1–17. Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., … Liu, T. Y. (2017). Lightgbm: A highly efficient gradient boosting decision tree. Advances in neural information processing systems , 30 . Koenker, R., & Bassett Jr, G. (1978). Regression quantiles. Econometrica: journal of the Econometric Society , 33–50. Liu, R., Gou, J., Han, J., Cao, K., Pan, X., Wang, D., & Wei, L. (2025). Optimizing groundwater potential assessment: uncertainty reduction through sample balancing and enhanced hybrid modeling. Stochastic Environmental Research and Risk Assessment , 1–21. Lundberg, S. M., & Lee, S. I. (2017). A unified approach to interpreting model predictions. Advances in neural information processing systems , 30 . Machiwal, D., Jha, M. K., & Mal, B. C. (2011). Assessment of groundwater potential in a semi-arid region of India using remote sensing, GIS and MCDM techniques. Water resources management , 25 (5), 1359–1386. Mogaji, K. A., Lim, H. S., & Abdullah, K. (2015). Regional prediction of groundwater potential using GIS-based Dempster–Shafer theory and frequency ratio models. Environmental Earth Sciences , 73(1), 327–346. Moran, P. A. (1950). Notes on continuous stochastic phenomena. Biometrika , 37 (1/2), 17–23. Mussa, M. M., Lohani, T. K., & Eshete, A. A. (2024). Evaluation of Groundwater Potential Zones Using GIS-Based Machine Learning Ensemble Models in the Gidabo Watershed, Ethiopia. Global Challenges , 8 (12), 2400137. Naghibi, S. A., Pourghasemi, H. R., & Dixon, B. (2016). GIS-based groundwater potential mapping using boosted regression tree, classification and regression tree, and random forest machine learning models in Iran. Environmental monitoring and assessment , 188 (1), 44. Panahi, M., Sadhasivam, N., Pourghasemi, H. R., Rezaie, F., & Lee, S. (2020). Spatial prediction of groundwater potential mapping based on convolutional neural network (CNN) and support vector regression (SVR). Journal of Hydrology , 588 , 125033. Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., … Chintala, S. (2019). Pytorch: An imperative style, high-performance deep learning library. Advances in neural information processing systems , 32 . Prasad, R. K., Mondal, N. C., Banerjee, P., Nandakumar, M. V., & Singh, V. S. (2008). Deciphering potential groundwater zone in hard rock through the application of GIS. Environmental geology , 55 (3), 467–475. Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., & Gulin, A. (2018). CatBoost: unbiased boosting with categorical features. Advances in neural information processing systems , 31 . Rana, M. S. P., Rahman, M. T., & Hassan, M. F. (2025). Mapping groundwater potential zone by robust machine learning algorithms & remote sensing techniques in agriculture dominated area, Bangladesh. Cleaner Water , 3 , 100064. Ray, S. K. (2025). Unveiling groundwater gems: A GIS-powered fusion of AHP and TOPSIS for mapping groundwater potential zones. Groundwater for Sustainable Development , 29 , 101431. Rey, S. J., & Anselin, L. (2007). PySAL: A Python library of spatial analytical methods. The Review of Regional Studies , 37(1), 5–27. Saha, A., & Pal, S. C. (2025). Delineation of groundwater potential zones using data-driven approaches: towards achieving sustainable groundwater management in drought-prone region of Eastern India. Environmental Monitoring and Assessment , 197 (10), 1–22. Saha, R., Chiravuri, S. S., Das, I. C., Kandrika, S., Kumranchat, V. K., Chauhan, P., & Chitikela, V. L. (2024). Urban aquifer health assessment and its management for sustainable water supply: an innovative approach using machine learning techniques. Groundwater for Sustainable Development , 25 , 101130. Saltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., & Tarantola, S. (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. Computer physics communications , 181 (2), 259–270. Sander, P. (2007). Lineaments in groundwater exploration: a review of applications and limitations. Hydrogeology journal , 15 (1), 71–74. Shaban, A., Khawlie, M., & Abdallah, C. (2006). Use of remote sensing and GIS to determine recharge potential zones: the case of Occidental Lebanon. Hydrogeology Journal , 14 (4), 433–443. Singhal, B. B. S., & Gupta, R. P. (2010). Applied hydrogeology of fractured rocks . Springer Science & Business Media. Soboĺ, I. M. (1993). Sensitivity estimates for nonlinear mathematical models. Math. Model. Comput. Exp. , 1 , 407. Tobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. Economic geography , 46 (sup1), 234–240. Todd, D. K., & Mays, L. W. (2004). Groundwater hydrology . John Wiley & Sons. Vapnik, V. (2013). The nature of statistical learning theory . Springer science & business media. Vaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., … Polosukhin, I. (2017). Attention is all you need. Advances in neural information processing systems , 30 . Vishwakarma, B. D., Zhang, J., & Sneeuw, N. (2021). Downscaling GRACE total water storage change using partial least squares regression. Scientific data , 8 (1), 95. Yeh, H. F., Cheng, Y. S., Lin, H. I., & Lee, C. H. (2016). Mapping groundwater recharge potential zone using a GIS approach in Hualian River, Taiwan. Sustainable Environment Research , 26 (1), 33–43. Additional Declarations No competing interests reported. Supplementary Files ESM1.pdf ESM1.pdf Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-9249309","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":623985787,"identity":"e240f70f-1efd-4503-b3f6-e5515cce4883","order_by":0,"name":"Suraj Dule","email":"","orcid":"","institution":"Gangadhar Meher University","correspondingAuthor":false,"prefix":"","firstName":"Suraj","middleName":"","lastName":"Dule","suffix":""},{"id":623985788,"identity":"bccf8248-5365-446c-b954-4127e29155c7","order_by":1,"name":"Arabinda Sharma","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA/UlEQVRIie3PMWvCQBTA8XcELsvZrDkE/QTCEyGT36TLdYlTxclJ9LrURejaYD9EP0LKQaY0rulmdgXBUaG9S52EJI4O95/u4H68ewA2213GIAYKizdX6jN2PADyEt9CSLSKDRlwCUQ2EjDkMxfmbAhALfHW319qMgWHR/tiu5ug31sqPWU2fKwifjYW6j0D6rVHg/4Hoh+kT5ok4bOsMilD1XoFxtchbTM8z4NYEyJVJeleiI8/iXtiZsqmqCd4IYg5pU5J8oYp/fRBKJaB4KvQ4eUuuZ4ianbppC11ZFMQnpuQw+5sPjYqtofZsHr9/36v7qL+uc1ms9ka+gM/OlyWd/RLQAAAAABJRU5ErkJggg==","orcid":"","institution":"Gangadhar Meher University","correspondingAuthor":true,"prefix":"","firstName":"Arabinda","middleName":"","lastName":"Sharma","suffix":""}],"badges":[],"createdAt":"2026-03-28 03:23:11","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-9249309/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-9249309/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":107486888,"identity":"f8d10999-7d1c-4528-94e3-41b1c2ab901b","added_by":"auto","created_at":"2026-04-22 02:39:12","extension":"jpeg","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":368516,"visible":true,"origin":"","legend":"\u003cp\u003eStudy area location map.\u003c/p\u003e","description":"","filename":"1.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/3442cc44331baca0509d86db.jpeg"},{"id":107487257,"identity":"48c1bbdd-7732-4ac9-930d-711e8397a463","added_by":"auto","created_at":"2026-04-22 02:40:14","extension":"jpeg","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":1045089,"visible":true,"origin":"","legend":"\u003cp\u003eSpatial distribution of input thematic layers: (a) Elevation, (b) Slope, (c) Rainfall, (d) Runoff, (e) Drainage Density, (f) NDVI, (g) Lineament Density, (h) Mean Water Level, and (a) Aquifer Depth, (j) LULC, (k) Geomorphology, (l) Soil Texture, and (m) Geology.\u003c/p\u003e","description":"","filename":"2.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/50d754effd8ebffcfea9c2d3.jpeg"},{"id":107486978,"identity":"34a33238-aaa1-46a3-bdea-19f979a12414","added_by":"auto","created_at":"2026-04-22 02:39:26","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":356841,"visible":true,"origin":"","legend":"\u003cp\u003eComprehensive methodological flowchart.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/211a101a35dff25d8c31c20d.png"},{"id":107487144,"identity":"a7f66027-966e-4373-a1d1-137501dd3d73","added_by":"auto","created_at":"2026-04-22 02:39:53","extension":"jpeg","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":197704,"visible":true,"origin":"","legend":"\u003cp\u003eSpatial Cross-Validation blocks (coloured clusters).\u003c/p\u003e","description":"","filename":"4.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/8ce14af00dbcccde4c3b3e60.jpeg"},{"id":107487148,"identity":"1429a749-2166-4afe-8375-660007d17903","added_by":"auto","created_at":"2026-04-22 02:39:53","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":150956,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 4 \u003c/strong\u003eThe radar chart showing performance trade-off.\u003c/p\u003e","description":"","filename":"5.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/1b575cab2051222e5b395bd3.png"},{"id":107486977,"identity":"41eab087-99d3-431b-b75d-45b2271fd07f","added_by":"auto","created_at":"2026-04-22 02:39:26","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":132166,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 5 \u003c/strong\u003eComparative\u003cstrong\u003e \u003c/strong\u003eVIF analysis plot: Base and XGB-SF.\u003c/p\u003e","description":"","filename":"6.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/949faa52d1efaec14ded6934.png"},{"id":107487060,"identity":"98fc11d4-9645-4e25-b913-018906d68629","added_by":"auto","created_at":"2026-04-22 02:39:45","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":591609,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 6 \u003c/strong\u003eSHAP Summary Plots for (a) XGBoost, (b) XGBoost SF, (c) CatBoost, (d) RF, (e) SVM, and (f) MLP DL.\u003c/p\u003e","description":"","filename":"7.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/d07e47ebb38a6330cb223efb.png"},{"id":107487059,"identity":"c51693f1-4578-4bc8-9ece-f3dede8ac462","added_by":"auto","created_at":"2026-04-22 02:39:45","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":163235,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 7\u003c/strong\u003e The Heatmap showing the Permutation Feature Importance across six model architectures.\u003c/p\u003e","description":"","filename":"8.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/15cda7c95af913f7c1d32a2c.png"},{"id":107487123,"identity":"d1439055-88b2-46e7-8574-566bb390849e","added_by":"auto","created_at":"2026-04-22 02:39:47","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":56844,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 8 \u003c/strong\u003eSobol Sensitivity Analysis (First-Order Indices \u003cem\u003eS\u003c/em\u003e\u003csub\u003e1\u003c/sub\u003e).\u003c/p\u003e","description":"","filename":"9.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/8f874b006955ac2a4acad34d.png"},{"id":107486887,"identity":"a247c14a-5226-467d-967f-e4210b9ccd09","added_by":"auto","created_at":"2026-04-22 02:39:12","extension":"jpeg","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":646327,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 9 \u003c/strong\u003eFinal GWPZ maps produced by the six models: (a) XGBoost, (b) XGBoost-SF, (c) CatBoost, (d) RF, (e) SVM, and (f) MLP DL\u003c/p\u003e","description":"","filename":"10.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/37bf79b817463ca495fc0ae1.jpeg"},{"id":107487152,"identity":"1b76f05b-98db-452a-a874-87087e6aee7a","added_by":"auto","created_at":"2026-04-22 02:39:54","extension":"png","order_by":11,"title":"Figure 11","display":"","copyAsset":false,"role":"figure","size":86141,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 10 \u003c/strong\u003ePercentage areal extent of GWPZ classes predicted by the comparative machine learning models.\u003c/p\u003e","description":"","filename":"11.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/c137e30d7cc7b2556a8b7684.png"},{"id":107487261,"identity":"e47ecaa9-a0c5-41b3-8afb-b271d653fc21","added_by":"auto","created_at":"2026-04-22 02:40:16","extension":"png","order_by":12,"title":"Figure 12","display":"","copyAsset":false,"role":"figure","size":126454,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 11 \u003c/strong\u003eAUC ROC curves for the used models.\u003c/p\u003e","description":"","filename":"12.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/7e08fcb05d5b41383d991f28.png"},{"id":107487271,"identity":"21e722bf-f802-439e-b09a-b6e562093221","added_by":"auto","created_at":"2026-04-22 02:40:18","extension":"png","order_by":13,"title":"Figure 13","display":"","copyAsset":false,"role":"figure","size":81079,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 12 \u003c/strong\u003eComparative Bar Chart of Spearman's Rho for all models.\u003c/p\u003e","description":"","filename":"13.png","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/6e68837d1aaf2a6e9ab855a4.png"},{"id":107487251,"identity":"c3529b8c-3c04-40aa-bdfc-d1c4ceb1e1f4","added_by":"auto","created_at":"2026-04-22 02:40:13","extension":"jpeg","order_by":14,"title":"Figure 14","display":"","copyAsset":false,"role":"figure","size":162841,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFig. 13 \u003c/strong\u003eThe GWPZ LightGB-based Quartile Uncertainty map.\u003c/p\u003e","description":"","filename":"14.jpeg","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/6533d9dfcb5f8c6db9584a84.jpeg"},{"id":107705399,"identity":"38a65815-8f98-4e36-b599-bd7a32fe9079","added_by":"auto","created_at":"2026-04-24 09:12:19","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1426065,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/34b517de-4ffe-4f59-bed1-35e287738ed2.pdf"},{"id":107487157,"identity":"7b42052c-1d3c-49c9-919f-8eeeff6b0354","added_by":"auto","created_at":"2026-04-22 02:39:55","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"supplement","size":184435,"visible":true,"origin":"","legend":"","description":"","filename":"ESM1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/6b3685dfa32ac02a95875f30.pdf"},{"id":107487154,"identity":"02266f91-3980-41b3-afdc-19cca3c08f2b","added_by":"auto","created_at":"2026-04-22 02:39:54","extension":"pdf","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":184435,"visible":true,"origin":"","legend":"","description":"","filename":"ESM1.pdf","url":"https://assets-eu.researchsquare.com/files/rs-9249309/v1/de608d7029bce1019dd8458b.pdf"}],"financialInterests":"No competing interests reported.","formattedTitle":"Modelling Groundwater Potential of Tropical Hard-Rock Aquifers using Spatial Cross-Validation augmented with Uncertainty and Sensitivity Analysis","fulltext":[{"header":"1. Introduction","content":"\u003cp\u003eGroundwater is a critical freshwater resource for domestic, agricultural, and industrial use, serving as a primary buffer against climate variability in semi-arid regions (Famiglietti 2014). This dependence is particularly acute in India, the world's largest consumer of groundwater, where the resource supports over 60% of irrigated agriculture and 85% of rural drinking water supplies (Mukherjee et al. 2015). The resource management challenge is most severe in hard-rock terrains, which cover nearly two-thirds of the Indian peninsular shield. Unlike extensive alluvial aquifers where water is stored in continuous pore spaces, groundwater in crystalline basements (e.g., granites, gneisses) occurs discontinuously. Its storage and movement are governed strictly by secondary porosity elements, including fractures, joints, and the variable thickness of the weathered regolith (Machiwal et al. \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2011\u003c/span\u003e). In these heterogeneous environments, aquifer boundaries are ill-defined and well yields can vary drastically over short distances, making the accurate delineation of Groundwater Potential Zones (GWPZ) a complex but indispensable prerequisite for artificial recharge planning.\u003c/p\u003e \u003cp\u003eTraditional groundwater exploration relies on precise but capital-intensive hydrogeological field surveys, which are logistically prohibitive for regional assessments (Jha et al. 2010). While early geospatial efforts utilizing Remote Sensing (RS) and GIS integrated surface indicators via heuristic, knowledge-driven approaches (e.g., AHP, TOPSIS) (Dule and Sharma \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2025\u003c/span\u003e; Kar et al. \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2025\u003c/span\u003e; Ray \u003cspan citationid=\"CR35\" class=\"CitationRef\"\u003e2025\u003c/span\u003e), these methods remain inherently subjective and struggle with complex, non-linear interactions (Jenifer and Jha 2017). Consequently, the field has shifted decisively toward data-driven Machine Learning (ML) algorithms, which demonstrate superior capability in modelling spatial heterogeneity (Rana et al. \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2025\u003c/span\u003e; Naghibi et al. \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2016\u003c/span\u003e; Saha and Pal \u003cspan citationid=\"CR37\" class=\"CitationRef\"\u003e2025\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eA primary deficiency in current literature is treating geospatial data as independent. Most studies utilize standard random k-fold cross-validation to tune and evaluate models. This approach violates the fundamental assumption of independent and identically distributed (i.i.d.) data, as geospatial variables are invariably governed by Tobler\u0026rsquo;s First Law of Geography (Tobler \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e1970\u003c/span\u003e), leading to spatial autocorrelation. Random splitting leads to \"spatial information leakage,\" where training and testing points are geographically adjacent and statistically dependent. This results in inflated, over-optimistic performance metrics that fail to generalize to unseen locations (Roberts et al. 2017).\u003c/p\u003e \u003cp\u003eSecondly, GWPZ mapping is often mis-framed as a binary classification task or simple regression, neglecting the ordinal nature of potential zones (e.g., distinguishing 'Excellent' from 'Good') required for prioritization. Few studies validate the preservation of this ordinal rank (e.g., using Spearman\u0026rsquo;s \u0026#120588;) in model outputs. Although recent studies have explored Deep Learning techniques such as Convolutional Neural Networks (CNNs) for spatial mapping (Panahi et al. \u003cspan citationid=\"CR30\" class=\"CitationRef\"\u003e2020\u003c/span\u003e), the potential of modern architectures to capture geological nuances via entity embeddings remains underutilized (Guo and Berkhahn \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). Most critically, standard ML outputs are deterministic, providing single predictions without quantifying epistemic uncertainty. In decision-making contexts, a prediction of 'High Potential' carries significantly different risks depending on whether the model is 90% confident or only 50% confident. The lack of spatially explicit uncertainty quantification renders many current GWPZ maps unsuitable for high-stakes resource allocation, where the cost of a dry borehole is substantial (Kosmidis et al. 2020). While recent studies have applied artificial intelligence (AI) and advanced machine learning techniques to hydrology, and others have successfully utilized optimization techniques such as Genetic Algorithms for broad regional mapping, these assessments often prioritize temporal trends over the rigorous quantification of spatial prediction uncertainty (Mussa et al. \u003cspan citationid=\"CR28\" class=\"CitationRef\"\u003e2024\u003c/span\u003e; Dikshit et al. \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Saha et al. \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e2024\u003c/span\u003e). Consequently, the simultaneous integration of spatial autocorrelation correction and uncertainty analysis remains largely unexplored in hard-rock terrains (Liu et al. \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2025\u003c/span\u003e; Vishwakarma et al. \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2021\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eTo address these gaps, this study presents a comparative machine learning framework for GWPZ mapping in the Khatra sub-division, West Bengal, driven by four methodological innovations. First, Spatial Cross-Validation was integrated with Optuna hyperparameter tuning to eliminate spatial autocorrelation bias and ensure robust model generalization. Second, model selection prioritized ordinal ranking fidelity (Spearman\u0026rsquo;s \u0026#120588;) alongside traditional binary metrics (AUC) to ensure proper hierarchical classification of potential zones. Third, deterministic predictions were transformed into risk-aware decision-support tools by generating spatially explicit uncertainty maps via Quantile LightGBM regression. Finally, the study evaluated advanced predictive architectures, including an entity-embedded Deep Learning network (MLP-DL) and a spatially augmented XGBoost (XGBoost-SF) designed to capture multi-scale topographic context.\u003c/p\u003e \u003cp\u003eBy benchmarking these advanced architectures alongside established ensemble methods (RF, CatBoost, SVM, and standard XGBoost), this work establishes a new standard for reliable, transparent, and uncertainty-aware groundwater modelling in complex hard-rock terrains.\u003c/p\u003e"},{"header":"2. Data Sources \u0026 Methodology","content":"\u003cdiv id=\"Sec3\" class=\"Section2\"\u003e \u003ch2\u003e2.1 Study Region\u003c/h2\u003e \u003cp\u003eThe present study centres on the Khatra Sub-division of Bankura District, West Bengal, India (Fig.\u0026nbsp;\u003cspan refid=\"Fig1\" class=\"InternalRef\"\u003e1\u003c/span\u003e), encompassing a total area of approximately 2,433 km\u0026sup2; between longitudes 86\u0026deg;36'52\" E to 87\u0026deg;19'04\" E and latitudes 22\u0026deg;37'61\" N to 23\u0026deg;16'56\" N. Administratively, the sub-division comprises eight Community Development (CD) blocks: Khatra, Indpur, Hirband, Ranibandh, Simlapal, Taldangra, Raipur, and Sarenga. The region marks the transition from the Chotanagpur Plateau to Gangetic alluvial plains.\u003c/p\u003e \u003cp\u003eThe area lies geologically on the erosional eastern fringe of the Chotanagpur Granite Gneiss Complex (CGGC), dominated by Archean crystalline rocks (granite gneiss, mica schist) with negligible primary porosity. Groundwater occurrence is strictly confined to secondary porosity elements like fractures and weather regolith (Central Ground Water Board [CGWB], 2022). Topographically, the region exhibits a westward-to-eastward gradient (400 m to 60 m Above Mean Sea Level or AMSL), featuring undulating pediments in the west that merge into gently sloping alluvial plains in the east (CGWB, 2022).\u003c/p\u003e \u003cp\u003eThe study area experiences a tropical sub-humid regime with an average annual rainfall of approximately 1,386 mm, the vast majority of which occurs during the southwest monsoon (June to September). The drainage architecture is controlled by two prominent river systems, the Kangsabati and the Shilabati. While the Kangsabati River is regulated by the Mukutmanipur Dam in the Khatra block to support canal irrigation, surface water availability remains critically insufficient for agricultural demands during the dry \u003cem\u003erabi\u003c/em\u003e and pre-\u003cem\u003ekharif\u003c/em\u003e seasons (CGWB, 2022). This scarcity, combined with high surface runoff facilitated by the undulating terrain, has driven an increased dependence on groundwater extraction via borewells, placing unsustainable stress on the limited storage capacity of the fractured aquifers.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec4\" class=\"Section2\"\u003e \u003ch2\u003e2.2 Geospatial Database and Ground Truth Inventory\u003c/h2\u003e \u003cp\u003eThe spatial database for this study, anchored by a Survey of India (SOI) topographic base map (1:50,000), integrates multi-source satellite imagery, ancillary maps, and hydrological field data. Thirteen thematic predictor layers were selected based on established literature (Naghibi et al., \u003cspan citationid=\"CR29\" class=\"CitationRef\"\u003e2016\u003c/span\u003e) and processed within a QGIS environment. Topographic and structural variables\u0026mdash;Elevation, Slope, Drainage Density, and Lineament Density\u0026mdash;were extracted from the ALOS PALSAR Digital Elevation Model (12.5 m resolution). Land Use/Land Cover (LULC) and the Normalized Difference Vegetation Index (NDVI) were derived from Landsat 9 OLI imagery (30 m, 2023) using the specified bands, with the LULC classification rigorously validated against ground control points via Google Earth Pro (Overall Accuracy: 90.8%, Kappa coefficient: 0.846). Geology and soil texture parameters were digitized from the Geological Survey of India (GSI) District Resource Map (1:250,000) and the National Bureau of Soil Survey and Land Use Planning (NBSS\u0026amp;LUP) soil texture map (1:50,000), while the geomorphological vector layer was sourced from the Bhukush GSI geoportal.\u003c/p\u003e \u003cp\u003eHydrometeorological forcing incorporated 2023 CHIRPS rainfall data (0.05\u0026deg;) and surface runoff, which was explicitly modelled using the NRCS-Curve Number (CN) method by integrating NASA HYSOG data, rainfall, and LULC. Finally, subsurface conditions, specifically Mean Water Level (2019\u0026ndash;2023) in the form of point data was interpolated via Inverse Distance Weighting (IDW) and the Aquifer Depth layer was digitised from aquifer deposition 2013 data from the India Water Resources Information System (WRIS).\u003c/p\u003e \u003cp\u003eTo establish the target variable, groundwater yield (Litres Per Minute, LPM) was derived from NRSC ISRO prospect maps based on extensive field observations (April\u0026ndash;May 2009 and January\u0026ndash;February 2010). Although dynamic surface predictors utilized recent 2023 imagery, this historical baseline remains hydrologically valid because groundwater accumulation in hard-rock terrains is primarily driven by static structural controls (e.g., geomorphology, geology) rather than short-term surface changes. From this dataset, 303 points were randomly sampled and categorized into five ordinal classes (1\u0026ndash;5 scale): Very Low (30\u0026ndash;50), Low (50\u0026ndash;100), Moderate (100\u0026ndash;200), High (200\u0026ndash;400), and Very High (400\u0026ndash;800 LPM), then partitioned into 80% training and 20% testing subsets (random seed: 42). To prevent spatial data leakage, a completely independent external validation dataset of 180 points was curated from the same source with strictly zero spatial overlap with the primary data\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec5\" class=\"Section2\"\u003e \u003ch2\u003e2.3 Hydrogeological Rationale of Thematic Predictors\u003c/h2\u003e \u003cp\u003eThe delineation of groundwater potential in hard-rock terrains depends on the non-linear interaction between surface boundary conditions and subsurface storage capacity. The thirteen selected thematic layers given in Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e and Table\u0026nbsp;\u003cspan refid=\"Tab1\" class=\"InternalRef\"\u003e1\u003c/span\u003e were categorized into three functional groups:\u003c/p\u003e \u003cp\u003eTopographic and Geomorphic Controls: Topography heavily dictates surface water retention. Lower elevations (ranging down to 1 m AMSL) facilitate water accumulation and longer residence times, favouring recharge, whereas higher elevations (up to 248 m AMSL) promote rapid runoff (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ea) (Todd \u0026amp; Mays, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2004\u003c/span\u003e). Similarly, slope exhibits an inverse relationship with infiltration; gentle slopes (\u0026lt;\u0026thinsp;5\u0026deg;) enhance residence time, while steep slopes (up to 32.53\u0026deg;) accelerate runoff (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eb). Geomorphologically, alluvial and flood plains offer high primary porosity, pediment complexes serve as shallow aquifers, and dissected hills act as primary runoff zones (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ek) (Machiwal et al., \u003cspan citationid=\"CR25\" class=\"CitationRef\"\u003e2011\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eHydrological Drivers: Rainfall (1168\u0026ndash;1360 mm) acts as the primary source of recharge, directly correlating with increased potential in permeable zones (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ec) (Adiat et al., \u003cspan citationid=\"CR1\" class=\"CitationRef\"\u003e2012\u003c/span\u003e). Conversely, surface runoff (1047\u0026ndash;1351 mm) and high drainage density (up to 3.14 km/km\u0026sup2;) serve as inverse indicators; high values signify impermeable lithology and low infiltration capacity (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ed, e) (Chowdary et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Yeh et al., \u003cspan citationid=\"CR49\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). Furthermore, Land Use/Land Cover (LULC) and NDVI (-0.75 to 0.88) dictate percolation efficiency; vegetated and agricultural lands reduce runoff velocity and increase soil porosity via root systems, whereas impervious built-up areas impede recharge (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003ef, j) (Sander, \u003cspan citationid=\"CR40\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Shaban et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2006\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eSubsurface and Structural Constraints: In crystalline Geological units like Granite Gneiss and Schist, groundwater movement relies entirely on secondary porosity, making high Lineament Density (up to 0.56 km/km\u0026sup2;) critical conduits for flow (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eg, m) (Mogaji et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Singhal \u0026amp; Gupta, \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Pedologically, coarse Soil Textures like sandy loam facilitate rapid recharge compared to clay-rich soils (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003el) (Hillel, \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). Subsurface capacity is indicated by Mean Water Level and Aquifer Depth; shallow water levels (\u0026lt;\u0026thinsp;5 mbgl) signify saturated valley fills, while thicker aquifers (90\u0026ndash;150 m) indicate larger storage reservoirs capable of sustaining extraction during dry periods (Fig.\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003eh, i) (Prasad et al., \u003cspan citationid=\"CR32\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Todd \u0026amp; Mays, \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2004\u003c/span\u003e).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig2\" class=\"InternalRef\"\u003e2\u003c/span\u003e Spatial distribution of input thematic layers: (a) Elevation, (b) Slope, (c) Rainfall, (d) Runoff, (e) Drainage Density, (f) NDVI, (g) Lineament Density, (h) Mean Water Level, and (a) Aquifer Depth, (j) LULC, (k) Geomorphology, (l) Soil Texture, and (m) Geology.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003cdiv class=\"gridtable\"\u003e\u003ctable float=\"Yes\" id=\"Tab1\" border=\"1\"\u003e \u003ccaption language=\"En\"\u003e \u003cdiv class=\"CaptionNumber\"\u003eTable 1\u003c/div\u003e \u003cdiv class=\"CaptionContent\"\u003e \u003cp\u003eData sources and characteristics of the thematic layers.\u003c/p\u003e \u003c/div\u003e \u003c/caption\u003e \u003ccolgroup cols=\"4\"\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c1\" colnum=\"1\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c2\" colnum=\"2\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c3\" colnum=\"3\"\u003e\u003c/div\u003e \u003cdiv align=\"left\" class=\"colspec\" colname=\"c4\" colnum=\"4\"\u003e\u003c/div\u003e \u003cthead\u003e \u003ctr\u003e \u003cth align=\"left\" colname=\"c1\"\u003e \u003cp\u003eParameters\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c2\"\u003e \u003cp\u003eSource\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c3\"\u003e \u003cp\u003eData Type\u003c/p\u003e \u003c/th\u003e \u003cth align=\"left\" colname=\"c4\"\u003e \u003cp\u003eNative Resolution/Scale\u003c/p\u003e \u003c/th\u003e \u003c/tr\u003e \u003c/thead\u003e \u003ctbody\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eElevation, Slope, Drainage Density, Lineament Density\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eALOS PALSAR from Alaska Satellite Facility (ASF) and Japan Aerospace Exploration Agency (JAXA). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://datapool.asf.alaska.edu\u003c/span\u003e\u003cspan address=\"https://datapool.asf.alaska.edu\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eRaster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e12.5 m\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eRainfall\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eClimate Hazards Group Infrared Precipitation with Stations (CHIRPS). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://data.chc.ucsb.edu\u003c/span\u003e\u003cspan address=\"https://data.chc.ucsb.edu\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eRaster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e0.05\u0026deg; x 0.05\u0026deg;\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eRunoff\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eObtained using NRCS-CN Method\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eRaster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e12.5 m\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eNDVI, LULC\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eLandsat 9 OLI (30 m). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://earthexplorer.usgs.gov\u003c/span\u003e\u003cspan address=\"https://earthexplorer.usgs.gov\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eRaster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e30 m\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eGeomorphology, Geology\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eGeological Survey of India (GSI) (1:250k). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://gsi.gov.in/home\u003c/span\u003e\u003cspan address=\"https://gsi.gov.in/home\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eVector and Raster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1:250k\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eSoil Texture\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eNational Bureau of Soil Survey \u0026amp; Land Use Planning (NBSS\u0026amp;LUP). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://icar-nbsslup.org.in\u003c/span\u003e\u003cspan address=\"https://icar-nbsslup.org.in\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eVector\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1:50k\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eMean Water Level, Aquifer Depth\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eIndia Water Resource Information System (WRIS). (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://indiawris.gov.in\u003c/span\u003e\u003cspan address=\"https://indiawris.gov.in\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eVector and Raster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003ePoints, 1:250k\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003ctr\u003e \u003ctd align=\"left\" colname=\"c1\"\u003e \u003cp\u003eGroundwater Prospect\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c2\"\u003e \u003cp\u003eThe National Remote Sensing Centre (NRSC) and ISRO. (\u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://maps.wbphed.gov.in\u003c/span\u003e\u003cspan address=\"https://maps.wbphed.gov.in\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e)\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c3\"\u003e \u003cp\u003eRaster\u003c/p\u003e \u003c/td\u003e \u003ctd align=\"left\" colname=\"c4\"\u003e \u003cp\u003e1:50k\u003c/p\u003e \u003c/td\u003e \u003c/tr\u003e \u003c/tbody\u003e \u003c/colgroup\u003e \u003c/table\u003e\u003c/div\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec6\" class=\"Section2\"\u003e \u003ch2\u003e2.4 Data Preprocessing and Algorithm Architectures\u003c/h2\u003e \u003cp\u003eTo ensure spatial consistency across the multi-source database (Fig.\u0026nbsp;\u003cspan refid=\"Fig3\" class=\"InternalRef\"\u003e3\u003c/span\u003e), all thematic layers were resampled to a unified 12.5 m resolution, using the ALOS PALSAR DEM as the reference master layer. Nearest-neighbour interpolation was employed for categorical layers to preserve discrete class fidelity, while bilinear interpolation was applied to continuous variables to maintain surface gradients in QGIS software. Following the data resampling and interpolation step, six distinct architectures were evaluated. Preprocessing pipelines and feature engineering were tailored specifically to the mathematical requirements of each algorithm.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cdiv id=\"Sec7\" class=\"Section3\"\u003e \u003ch2\u003e2.4.1 Random Forest (RF)\u003c/h2\u003e \u003cp\u003eRandom Forest, developed by Leo Breiman (\u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2001\u003c/span\u003e), is an ensemble learning method that constructs a multitude of decision trees at training time. It operates on the principle of bagging (bootstrap aggregating), where each tree is trained on a random subset of the data sampled with replacement. Node splitting within the individual trees is governed by impurity metrics such as the Gini index (Eq.\u0026nbsp;\u003cspan refid=\"Equ1\" class=\"InternalRef\"\u003e1\u003c/span\u003e):\u003cdiv id=\"Equ1\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ1\" name=\"EquationSource\"\u003e\n$$\\:Gini=1-{\\sum\\:}_{i=1}^{C}{\\left({p}_{i}\\right)}^{2}$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e1\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003eWhere \u003cem\u003eC\u003c/em\u003e is the number of classes and \u003cem\u003ep\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e is the proportion of class \u003cem\u003ei\u003c/em\u003e. RF introduces a second layer of randomness by selecting a random subset of features at each split point (controlled by max_features). This decorrelates the trees, significantly reducing the variance of the model without increasing bias, making it highly robust to noise and overfitting. For this, in RF and the other tree-based ensembles, Min-Max scaling [0, 1] is used for continuous features, while categorical variables are integer-encoded. Because the model is dealing with a continuous regression task, the algorithm must aggregate the independent outputs of the forest into a single continuous value. It achieves this by calculating the arithmetic mean of the predictions from all \u003cem\u003eN\u003c/em\u003e individual trees (\u003cem\u003ef\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e), yielding the final groundwater potential prediction (\u003cem\u003eF\u003c/em\u003e(\u003cem\u003ex\u003c/em\u003e)) for any given point (\u003cem\u003ex\u003c/em\u003e), using Eq.\u0026nbsp;\u003cspan refid=\"Equ2\" class=\"InternalRef\"\u003e2\u003c/span\u003e:\u003cdiv id=\"Equ2\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ2\" name=\"EquationSource\"\u003e\n$$\\:F\\left(x\\right)=\\frac{1}{N}{\\sum\\:}_{i=1}^{N}{f}_{i}\\left(x\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e2\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec8\" class=\"Section3\"\u003e \u003ch2\u003e2.4.2 Extreme Gradient Boosting (Standard XGBoost)\u003c/h2\u003e \u003cp\u003eUnlike RF, which builds independent trees in parallel, Gradient Boosting builds an ensemble sequentially. Each new tree is trained to predict the pseudo-residuals (errors) of the preceding ensemble, effectively converting weak learners into a strong learner (Friedman, \u003cspan citationid=\"CR13\" class=\"CitationRef\"\u003e2001\u003c/span\u003e). eXtreme Gradient Boosting (XGBoost) is an advanced implementation optimized for speed and performance (Chen, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2016\u003c/span\u003e), minimizing a regularized objective function (\u003cem\u003eObj\u003c/em\u003e) that combines a convex loss function (\u003cem\u003eL\u003c/em\u003e) with a penalty term (\u0026#120570;) for model complexity (Eq.\u0026nbsp;\u003cspan refid=\"Equ3\" class=\"InternalRef\"\u003e3\u003c/span\u003e):\u003cdiv id=\"Equ3\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ3\" name=\"EquationSource\"\u003e\n$$\\:Obj={\\sum\\:}_{i=1}^{n}L\\left({y}_{i},\\widehat{{y}_{i}}\\right)+{\\sum\\:}_{m=1}^{M}{\\Omega\\:}\\left({f}_{m}\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e3\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cem\u003en\u003c/em\u003e represents the total number of training instances, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\text{y}}_{\\text{i}}\\)\u003c/span\u003e\u003c/span\u003e is the actual observed groundwater yield, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\widehat{{\\text{y}}_{\\text{i}}}\\)\u003c/span\u003e\u003c/span\u003e is the model's prediction. The first term, \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\text{L}\\left({\\text{y}}_{\\text{i}},\\widehat{{\\text{y}}_{\\text{i}}}\\right)\\)\u003c/span\u003e\u003c/span\u003e, calculates the training loss by measuring the difference between the actual and predicted values. The second term introduces the regularization penalty \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Omega\\:}\\)\u003c/span\u003e\u003c/span\u003e applied across all \u003cem\u003eM\u003c/em\u003e sequential trees \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{f}_{m}\\)\u003c/span\u003e\u003c/span\u003e in the ensemble. Specifically, the penalty is defined as \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\Omega\\:}\\left({f}_{m}\\right)={\\gamma\\:}\\text{T}+\\frac{1}{2}{\\lambda\\:}{\\left|\\left|\\text{w}\\right|\\right|}^{2}\\)\u003c/span\u003e\u003c/span\u003e, where \u003cem\u003eT\u003c/em\u003e represents the total number of leaves in the tree, \u003cem\u003ew\u003c/em\u003e represents the leaf weights (scores), and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\gamma\\:}\\)\u003c/span\u003e\u003c/span\u003e and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\lambda\\:}\\)\u003c/span\u003e\u003c/span\u003e are parameters that control the severity of the regularization. This explicit penalization is critical for preventing the model from overfitting on limited, high-dimensional training data. This study configured the model with tree_method='hist' and enable_categorical=True to natively handle the integer-coded categorical features with efficiency.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec9\" class=\"Section3\"\u003e \u003ch2\u003e2.4.3 Spatially-Augmented XGBoost\u003c/h2\u003e \u003cp\u003eStandard machine learning algorithms remain inherently \"point-based,\" often failing to account for spatial autocorrelation as described by Tobler\u0026rsquo;s First Law of Geography (Tobler, \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e1970\u003c/span\u003e). To address this limitation, a spatially-augmented variant (XGBoost-SF) was developed. It shares the identical mathematical foundation as Standard XGBoost but incorporates engineered spatial textures. Three features were derived from the Elevation raster: a Focal Mean (5 x 5-pixel window) to smooth noise and highlight regional trends; a Focal Standard Deviation to quantify local topographic roughness; and a Sobel Gradient (3 x 3 filter) to assess multi-directional slope steepness. By expanding the input matrix with these textures, the gradient boosting objective function explicitly learns regional topographic transitions rather than treating each pixel as an isolated entity.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec10\" class=\"Section3\"\u003e \u003ch2\u003e2.4.4 CatBoost\u003c/h2\u003e \u003cp\u003eCategorical Boosting (CatBoost) is a GBDT algorithm specifically engineered to handle categorical data (Prokhorenkova et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2018\u003c/span\u003e). Its primary innovation is \"ordered target statistics\"- a permutation-based strategy for encoding categorical features that prevents target leakage and allows for the full, unbiased utilization of high-cardinality discrete features. It also grows \"symmetric\" or \"oblivious\" trees, which act as a form of regularization and significantly speed up inference. This architectural symmetry was deemed highly relevant given the established importance of categorical layers like Geomorphology and Geology in hydrogeological modelling.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec11\" class=\"Section3\"\u003e \u003ch2\u003e2.4.5 Support Vector Machine (SVM)\u003c/h2\u003e \u003cp\u003eSupport Vector Regression (SVR) is a kernel-based method rooted in statistical learning theory (Vapnik, 1995). Unlike traditional least-squares regression, SVR seeks a flat function \u003cem\u003ef\u003c/em\u003e (\u003cem\u003ex\u003c/em\u003e) deviating from the target \u003cem\u003ey\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e by at most a margin \u0026#120656; (Drucker et al., 1997). This \u0026#120656;-insensitive loss function only penalizes errors falling outside this specific margin. To model non-linear spatial relationships, data is mapped to a higher-dimensional feature space using the Radial Basis Function (RBF) kernel (Eq.\u0026nbsp;\u003cspan refid=\"Equ4\" class=\"InternalRef\"\u003e4\u003c/span\u003e):\u003cdiv id=\"Equ4\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ4\" name=\"EquationSource\"\u003e\n$$\\:\\left({x}_{i},{x}_{j}\\right)=\\text{exp}\\left(-\\gamma\\:{\\left|\\left|{x}_{i}-{x}_{j}\\right|\\right|}^{2}\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e4\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{x}_{i}\\:\\)\u003c/span\u003e\u003c/span\u003eand \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{x}_{j}\\)\u003c/span\u003e\u003c/span\u003e are distinct geographical feature vectors, and \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\left|\\left|{x}_{i}-{x}_{j}\\right|\\right|}^{2}\\)\u003c/span\u003e\u003c/span\u003e calculates the squared Euclidean distance. The kernel coefficient \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{\\gamma\\:}\\:\\)\u003c/span\u003e\u003c/span\u003e(gamma) defines the similarity radius of individual support vectors. Together with the regularization penalty \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:C\\)\u003c/span\u003e\u003c/span\u003e (for margin violations) and \u0026#120656;, it dictates the model's complexity and ability to generalize. Unlike tree ensembles, SVM strictly requires One-Hot Encoding for categorical inputs and Standard Scaling (zero mean, unit variance) to prevent magnitude bias in distance calculations (Hastie et al., \u003cspan citationid=\"CR17\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; G\u0026eacute;ron, \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2022\u003c/span\u003e). Due to its \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:\\mathcal{O}\\left({\\text{n}}^{3}\\right)\\)\u003c/span\u003e\u003c/span\u003e computational complexity, the SVR was trained on a random sub-sample of 8,000 points.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec12\" class=\"Section3\"\u003e \u003ch2\u003e2.4.6 Multi-Layer Perceptron Deep Learning (MLP-DL)\u003c/h2\u003e \u003cp\u003eA custom tabular Deep Learning (DL) model was implemented in PyTorch to process mixed data types. To prevent the sparsity of one-hot encoding, high-cardinality categorical variables (e.g., Geology) were passed through nn.Embedding layers (Guo \u0026amp; Berkhahn, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). This generated dense, low-dimensional vector representations capturing semantic class similarities. These embeddings were concatenated with continuous inputs and routed through a simplified Self-Attention mechanism (Vaswani et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2017\u003c/span\u003e) to dynamically weigh feature importance. The core backbone consisted of a Residual MLP with four hidden layers ([256, 256, 128, 64]) utilizing LeakyReLU activation (Eq.\u0026nbsp;\u003cspan refid=\"Equ5\" class=\"InternalRef\"\u003e5\u003c/span\u003e):\u003cdiv id=\"Equ5\" class=\"Equation\"\u003e\u003cdiv format=\"TEX\" class=\"mathdisplay\" id=\"FileID_Equ5\" name=\"EquationSource\"\u003e\n$$\\:f\\left(x\\right)={max}\\left(0.01x,x\\right)$$\u003c/div\u003e\u003cdiv class=\"EquationNumber\"\u003e5\u003c/div\u003e\u003c/div\u003e\u003c/p\u003e \u003cp\u003ewhere \u003cem\u003ex\u003c/em\u003e represents the input signal to the neuron. Unlike standard ReLU, LeakyReLU applies a small positive slope (0.01) to negative values, ensuring continuous gradient flow and mitigating the vanishing gradient problem. Finally, the activated data passed through Dropout (0.1) layers for regularization.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec13\" class=\"Section2\"\u003e \u003ch2\u003e2.5 Hyperparameter Optimization and Spatial Cross-Validation\u003c/h2\u003e \u003cp\u003eExcept for the MLP-DL, all architectures were optimized via Optuna using a 50-trial Tree-Structured Parzen Estimator (TPE) (Akiba et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Online Resource 1). To prevent the optimistic bias and spatial data leakage typical of random sampling, Spatial Cross-Validation (SpatialKFold) was integrated directly into the objective function. By partitioning the study area into geometrically disjoint blocks (2x2 geographic grid) via coordinate k-means clustering (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), models were forced to learn generalized hydrogeological rules rather than memorizing local coordinates. This realistically simulates exploration in uncharted terrains. Additionally, the RF and SVM objective functions incorporated a penalty for train-validation divergence to explicitly enforce algorithmic generalization.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003eFigure\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e Spatial Cross-Validation blocks (coloured clusters).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec14\" class=\"Section2\"\u003e \u003ch2\u003e2.6 Feature Diagnostics\u003c/h2\u003e \u003cp\u003eTo ensure the model relies on physically valid hydro-climatic relationships rather than spurious correlations, we implemented a dual-stage diagnostic framework.\u003c/p\u003e \u003cp\u003eThe first stage focused on multicollinearity screening, specifically restricted to continuous climatic and topographic predictors (e.g., Rainfall, Slope, Elevation). Variance Inflation Factor (VIF) analysis, a linear diagnostic derived from Ordinary Least Squares (OLS) regression, was utilized for this purpose. Applying VIF to high-cardinality categorical features like Geology would require one-hot encoding, which introduces artificial multicollinearity (the 'dummy variable trap') without providing meaningful insight. By restricting VIF solely to continuous raster, physical redundancies\u0026mdash;such as those between Elevation and Slope\u0026mdash;were successfully identified without introducing numerical artefacts.\u003c/p\u003e \u003cp\u003eThe second stage evaluated global feature importance to assess the predictive power of the full mixed-type feature set. For this, model-agnostic techniques, specifically SHAP (SHapley Additive exPlanations) and Permutation Importance, were employed. Unlike VIF, these methods quantify the marginal contribution of a feature to the model's decision boundary rather than its linear correlation. This approach allows for a direct, unbiased comparison between continuous drivers (e.g., NDVI) and categorical controls (e.g., Geomorphology) on a unified scale, effectively capturing both linear and non-linear interactions that traditional diagnostics typically miss (Lundberg \u0026amp; Lee, \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Fisher et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2019\u003c/span\u003e).\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec15\" class=\"Section2\"\u003e \u003ch2\u003e2.7 Model Interpretability\u003c/h2\u003e \u003cp\u003eA multi-stage validation workflow was implemented to assess model robustness, stability, and feature logic beyond simple accuracy metrics. To establish baseline predictive performance, internal validation on the 20% hold-out set utilized the Coefficient of Determination (\u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e), Root Mean Square Error (RMSE), and Mean Absolute Error (MAE). Generalization was then rigorously tested against an independent external ground truth dataset. This external phase employed the Area Under the Receiver Operating Characteristic Curve (AUC-ROC) to evaluate binary discrimination capability, alongside Spearman's Rank Correlation (\u0026#120588;) to assess the models' ordinal fidelity across the discrete GWPZ classes.\u003c/p\u003e \u003cp\u003eBeyond standard error metrics, the workflow explicitly quantified algorithmic uncertainty and input sensitivity. Model stability was evaluated via a 60-iteration bootstrap analysis (Efron \u0026amp; Tibshirani, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e1994\u003c/span\u003e), which calculated the mean per-point standard deviation to establish the deterministic consensus of each architecture. Finally, non-linear relationships were mathematically decomposed using Sobol Global Sensitivity Analysis (Soboĺ, \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e1993\u003c/span\u003e; Saltelli et al., \u003cspan citationid=\"CR39\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). By calculating First-Order (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{i}\\)\u003c/span\u003e\u003c/span\u003e) and Total-Order (\u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:{S}_{T}\\)\u003c/span\u003e\u003c/span\u003e) indices for the continuous features, this analysis successfully isolated the variance contribution attributable strictly to individual inputs versus complex multi-feature interactions.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec16\" class=\"Section2\"\u003e \u003ch2\u003e2.8 Spatial Autocorrelation Analysis\u003c/h2\u003e \u003cp\u003eTo evaluate spatial independence\u0026mdash;a fundamental assumption in geospatial modelling\u0026mdash;autocorrelation was assessed using the esda package (Rey \u0026amp; Anselin, \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2007\u003c/span\u003e). A spatial weights matrix was defined via \u003cem\u003ek\u003c/em\u003e-Nearest Neighbours \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:(k=5)\\)\u003c/span\u003e\u003c/span\u003e with significance tested across 999 simulations. The analysis was two-fold: (a) Input Data: Global Moran's I (Moran, \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e1950\u003c/span\u003e) was calculated for the raw GWPZ values. A significant clustering \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:(p\u0026lt;0.05)\\)\u003c/span\u003e\u003c/span\u003e confirmed the necessity of the Spatial Cross-Validation strategy. (b) Residual Analysis: Post-modelling, both Global Moran's I and Geary's C (Geary, \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e1954\u003c/span\u003e) were computed on the test set residuals (\u003cem\u003ee\u003c/em\u003e\u003csub\u003e\u003cem\u003ei\u003c/em\u003e\u003c/sub\u003e). A non-significant result \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:(p\u0026gt;0.05)\\)\u003c/span\u003e\u003c/span\u003e was required to confirm that the model successfully captured the underlying spatial structure without leaving systematic, clustered errors.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec17\" class=\"Section2\"\u003e \u003ch2\u003e2.9 Spatial Uncertainty Analysis\u003c/h2\u003e \u003cp\u003eTo move beyond a single mean prediction and quantify epistemic (model) uncertainty, a quantile regression workflow was implemented. Three LightGBM models (Ke et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2017\u003c/span\u003e) were trained using the quantile objective to predict the 5th, 50th, and 95th percentiles (Koenker \u0026amp; Bassett Jr, \u003cspan citationid=\"CR22\" class=\"CitationRef\"\u003e1978\u003c/span\u003e). LightGBM was selected for its computational efficiency in estimating conditional quantiles. The final uncertainty map was generated by calculating the 90% Prediction Interval (PI) width \u003cspan class=\"InlineEquation\"\u003e\u003cspan class=\"mathinline\"\u003e\\(\\:({\\mathcal{Q}}_{95}-{\\mathcal{Q}}_{05})\\)\u003c/span\u003e\u003c/span\u003e for each pixel, allowing stakeholders to distinguish high-confidence zones from areas requiring further field verification spatially. We extracted the 95% Prediction Intervals (PI) for the independent test set to assess the coverage probability, comparing the predicted mean yield against the width of the uncertainty bound.\u003c/p\u003e \u003c/div\u003e"},{"header":"3. Results and Discussion","content":"\u003cp\u003eThis study rigorously evaluated six distinct machine learning architectures for Groundwater Potential Zonation (GWPZ): RF, CatBoost, SVM, MLP-DL, Standard XGBoost, and XGBoost-SF. All models underwent identical Optuna hyperparameter tuning and spatial cross-validation (CV) to mitigate spatial autocorrelation. The assessment covers spatial analysis, internal predictive mechanics, and definitive external validation.\u003c/p\u003e \u003cdiv id=\"Sec19\" class=\"Section2\"\u003e \u003ch2\u003e3.1 Internal Predictive Mechanisms\u003c/h2\u003e \u003cp\u003eModel predictive capability was initially evaluated on the 20% hold-out test set to assess the trade-off between pure accuracy and algorithmic generalization. The spatially-augmented XGBoost-SF demonstrated superior internal accuracy, yielding the highest Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e (0.853) and lowest errors (RMSE\u0026thinsp;=\u0026thinsp;0.508, MAE\u0026thinsp;=\u0026thinsp;0.403) with a minimal overfitting gap of 1.1% compared to its Training \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e (0.864). This outperformed the standard XGBoost (Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.846, Train \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.876, gap\u0026thinsp;=\u0026thinsp;3.0%, RMSE\u0026thinsp;=\u0026thinsp;0.520, MAE\u0026thinsp;=\u0026thinsp;0.425), confirming that the integration of the engineered focal mean elevation feature provided highly valuable spatial context. Visualized via a radar chart (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e), XGBoost-SF achieved the optimal balance of high accuracy and low error.\u003c/p\u003e \u003cp\u003eThe remaining tree-based ensembles also exhibited strong predictive capabilities and small generalization gaps, indicating that regularization parameters optimized via Optuna successfully prevented memorization. CatBoost achieved a Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e of 0.834 (Train \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.863, gap\u0026thinsp;=\u0026thinsp;2.9%, RMSE\u0026thinsp;=\u0026thinsp;0.540, MAE\u0026thinsp;=\u0026thinsp;0.419), followed closely by Random Forest with a Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e of 0.823 (Train \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.865, gap\u0026thinsp;=\u0026thinsp;4.2%, RMSE\u0026thinsp;=\u0026thinsp;0.558, MAE\u0026thinsp;=\u0026thinsp;0.432).\u003c/p\u003e \u003cp\u003eConversely, the MLP-DL and SVM architectures displayed highly constrained learning, avoiding overfitting almost entirely. The MLP-DL achieved a Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e of 0.831 (RMSE\u0026thinsp;=\u0026thinsp;0.544, MAE\u0026thinsp;=\u0026thinsp;0.423) with a negligible overfitting gap of 0.5% against its Training \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e (0.836), suggesting the Dropout layers and entity embeddings allowed the network to learn a highly generalized representation. Remarkably, SVM displayed a negative overfitting gap (-2.0%), performing better on the test set (Test \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.786, RMSE\u0026thinsp;=\u0026thinsp;0.613, MAE\u0026thinsp;=\u0026thinsp;0.521) than the training sub-sample (Train \u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.766). While this testifies to the powerful global regularization of the SVM\u0026rsquo;s objective function\u0026mdash;establishing it as a robust baseline for data-scarce regions\u0026mdash;it resulted in lower overall peak accuracy compared to the boosting models.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec20\" class=\"Section2\"\u003e \u003ch2\u003e3.2 Hydrogeological Drivers and Feature Logic\u003c/h2\u003e \u003cp\u003eFeature importance was triangulated using Variance Inflation Factor (VIF), SHAP (SHapley Additive exPlanations), Permutation Importance, and Sobol Sensitivity Analysis to ensure physical validity and algorithm transparency.\u003c/p\u003e \u003cdiv id=\"Sec21\" class=\"Section3\"\u003e \u003ch2\u003e3.2.1 Multicollinearity Test\u003c/h2\u003e \u003cp\u003eVIF analysis (Fig.\u0026nbsp;\u003cspan refid=\"Fig5\" class=\"InternalRef\"\u003e5\u003c/span\u003e; Online Resource 1) confirmed the physical independence of the base thematic predictors, with all continuous variables\u0026mdash;including Runoff (VIF\u0026thinsp;=\u0026thinsp;5.25) and Rainfall (VIF\u0026thinsp;=\u0026thinsp;4.56)\u0026mdash;registering well below the severe multicollinearity threshold of 10. Conversely, the introduction of engineered spatial textures in the XGBoost-SF model inherently generated extreme mathematical redundancy, specifically between Slope and GWPZ Gradient (VIF\u0026thinsp;\u0026gt;\u0026thinsp;71,000) and Elevation and GWPZ Mean (VIF\u0026thinsp;\u0026gt;\u0026thinsp;8,000).\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec22\" class=\"Section3\"\u003e \u003ch2\u003e3.2.2 Feature Importance and Sensitivity\u003c/h2\u003e \u003cp\u003eSynthesizing global SHAP values (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003e) and Permutation Mean \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e Drop (Fig.\u0026nbsp;\u003cspan refid=\"Fig7\" class=\"InternalRef\"\u003e7\u003c/span\u003e) reveals how algorithms navigated collinearity and classified predictors into hydrogeological tiers. First, tree-based models successfully managed dataset redundancies without manual feature elimination by consistently prioritizing Runoff over Rainfall (mean SHAP\u0026thinsp;\u0026lt;\u0026thinsp;0.012) (Breiman, \u003cspan citationid=\"CR3\" class=\"CitationRef\"\u003e2001\u003c/span\u003e; Chen, \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). This is physically sound, as Runoff acts as a localized synthesis of soil texture, land use, and slope, whereas rainfall provides a uniform regional input. Beyond resolving collinearity, a dominant triad of Geomorphology, LULC, and NDVI emerged universally. Geomorphology acted as the absolute primary control (SHAP: 0.681\u0026ndash;0.804; Mean \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e Drop\u0026thinsp;\u0026gt;\u0026thinsp;0.95), dictating static aquifer boundary conditions, while LULC and NDVI served as critical proxies for dynamic surface permeability.\u003c/p\u003e \u003cp\u003eBecause Geomorphology dominates this discrete space, Sobol Sensitivity Analysis (Fig.\u0026nbsp;\u003cspan refid=\"Fig8\" class=\"InternalRef\"\u003e8\u003c/span\u003e) was utilized to isolate the variance from continuous physical variables. The resulting first-order indices (\u003cem\u003eS\u003c/em\u003e\u003csub\u003e1\u003c/sub\u003e) identified NDVI as the primary contributor (\u003cem\u003eS\u003c/em\u003e\u003csub\u003e1\u003c/sub\u003e\u0026thinsp;=\u0026thinsp;0.6667), followed by Drainage Density (\u003cem\u003eS\u003c/em\u003e\u003csub\u003e1\u003c/sub\u003e\u0026thinsp;=\u0026thinsp;0.0962) and Elevation (\u003cem\u003eS\u003c/em\u003e\u003csub\u003e1\u003c/sub\u003e\u0026thinsp;=\u0026thinsp;0.0517). Thus, while Geomorphology sets the structural container, vegetation health (NDVI) fine-tunes near-surface soil moisture and shallow groundwater availability during dry seasons. Furthermore, analysis of the spatially augmented XGBoost-SF (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003eb) validated the spatial feature engineering methodology. The focal mean texture feature (GWPZ_MEAN) emerged as the fourth most important predictor (SHAP\u0026thinsp;=\u0026thinsp;0.077), demonstrating that providing the gradient boosting algorithm with 5x5-pixel regional topographic awareness significantly influences its decision tree splits.\u003c/p\u003e \u003cp\u003eNon-tree algorithms exhibited unique inductive biases. The MLP-DL model (Fig.\u0026nbsp;\u003cspan refid=\"Fig6\" class=\"InternalRef\"\u003e6\u003c/span\u003ef) elevated Geology (SHAP\u0026thinsp;=\u0026thinsp;0.150) and Soil Texture to much greater prominence than the tree-based models, where Geology caused an anomalous 0.094 drop in \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e. This is directly attributed to the nn.Embedding layers (Guo \u0026amp; Berkhahn, \u003cspan citationid=\"CR16\" class=\"CitationRef\"\u003e2016\u003c/span\u003e), which allowed the neural network to learn dense, multi-dimensional vector representations and exploit complex semantic relationships that are missed by standard integer-encoding. Conversely, the SVM demonstrated a distinct kernel bias towards hydrologic flow. Because SVM required one-hot encoding, the categorical dominance of Geomorphology was fragmented across multiple dummy variables (e.g., Geomorphology_7, Geomorphology_1). This fragmentation allowed Runoff to become a highly influential feature (SHAP\u0026thinsp;=\u0026thinsp;0.205; Permutation \u003cem\u003eR\u003c/em\u003e\u003csup\u003e2\u003c/sup\u003e Drop\u0026thinsp;=\u0026thinsp;0.112), demonstrating the RBF kernel's sensitivity to continuous, scaled hydrologic gradients.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec23\" class=\"Section2\"\u003e \u003ch2\u003e3.3 Spatial Zonation Analysis\u003c/h2\u003e \u003cp\u003eTo evaluate how the statistical differences identified in the previous section translate into physical zoning, the continuous Groundwater Potential Index (GWPI) generated by all six optimized models was classified using the Natural Breaks (Jenks) method to minimize intra-class variance (Jenks, \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e1967\u003c/span\u003e). This aligned with the ordinal-rank validation strategy, classifying the study area into five ordinal zones: \"Very Low,\" \"Low,\" \"Moderate,\" \"High,\" and \"Very High.\" Visual inspection of the classified maps (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003e) alongside their comparative areal distributions (Fig.\u0026nbsp;\u003cspan refid=\"Fig10\" class=\"InternalRef\"\u003e10\u003c/span\u003e) reveals critical insights into model stability, geomorphic adherence, and algorithmic inductive bias.\u003c/p\u003e \u003cp\u003eXGBoost (Standard \u0026amp; SF): Both the Standard XGBoost (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ea) and the XGBoost-SF (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003eb) produced hydro-geologically plausible maps with remarkably balanced areal distributions. The combined 'Very High' and \u0026lsquo;High\u0026rsquo; zones (4.9% in Standard, 5.2% in SF) were clearly delineated along primary river valleys characterized by gentle slopes and saturated sediments, creating smooth, logical transitions to 'Moderate' weathered pediments. The \u0026lsquo;Moderate\u0026rsquo; zone constitutes approximately 34.2% (Standard) to 36.3% (SF) of the area. Hydrologically, this aligns perfectly with the reality of hard-rock terrains, where true high-yield aquifers are rare, but a significant portion of the landscape consists of weathered mantles offering moderate, seasonally variable yields. The remaining area falls under the 'Low' (31% in Standard, 29.8% in SF) and 'Very Low' (29.9% in Standard, 28.7% in SF) zones, which are heavily concentrated in the south-western and south-to-southeastern regions. Notably, the XGB-SF map exhibited slightly expanded 'Moderate' zones in upland plateaus due to the engineered focal mean feature.\u003c/p\u003e \u003cp\u003eCatBoost: This model generated a highly conservative landscape dominated by \u0026lsquo;Low\u0026rsquo; (63.6%) and \u0026lsquo;Very Low\u0026rsquo; (26.5%) zones (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ec). High potential was strictly confined to narrow river corridors, with \u0026lsquo;High\u0026rsquo; and \u0026lsquo;Very High\u0026rsquo; zones representing a negligible 5% combined. The \u0026lsquo;Moderate\u0026rsquo; potential zone was exceptionally limited (4.9%), restricted to areas in proximity to the southern river sections.\u003c/p\u003e \u003cp\u003eMLP-DL: Conversely, MLP-DL demonstrated the most optimistic areal distribution, shifting the spectral curve toward higher potentiality (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ef). It identified a combined 27.2% of the region as \u0026lsquo;High\u0026rsquo; to \u0026lsquo;Very High\u0026rsquo; potential and allocated the largest proportion of any model (40.8%) to the \u0026lsquo;Moderate\u0026rsquo; class. The remaining elevated forested land was classified as \u0026lsquo;Low\u0026rsquo; (15.9%) and \u0026lsquo;Very Low\u0026rsquo; (16%). This tendency to over-predict in marginal areas and produce fragmented outputs reflects the stochastic training nature of neural networks on limited datasets, which often fail to suppress noisy signals.\u003c/p\u003e \u003cp\u003eRF: The RF model allocated 5% and 12.3% of the study area to the \u0026lsquo;Very High\u0026rsquo; and \u0026lsquo;High\u0026rsquo; zones, respectively (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ed). While these were generally situated close to waterbodies, the model anomalously predicted high potential across the south-western hilly tract. Almost the entire northern section was identified as \u0026lsquo;Moderate\u0026rsquo; (20.9%), whereas a vast area covering the central, eastern, and select western and southern regions was classified as \u0026lsquo;Low\u0026rsquo; (38.6%) and \u0026lsquo;Very Low\u0026rsquo; (23.2%).\u003c/p\u003e \u003cp\u003eSVM: Displaying large, contiguous blocks of potential, the SVM map (Fig.\u0026nbsp;\u003cspan refid=\"Fig9\" class=\"InternalRef\"\u003e9\u003c/span\u003ee) reflects the \"gap penalty\" and RBF kernel properties of the algorithm (Vapnik, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). While this generalization creates broad zoning, it glosses over the local hydrogeological nuances required for site-specific exploration. It predicted 4.6% and 8.5% of the area as \u0026lsquo;Very High\u0026rsquo; and \u0026lsquo;High\u0026rsquo; potential, primarily within southern, western, and eastern river valleys. Adjoining areas were classified as \u0026lsquo;Moderate\u0026rsquo; (28.7%), while rocky uplands and elevated hills defaulted to the \u0026lsquo;Low\u0026rsquo; (34.2%) and \u0026lsquo;Very Low\u0026rsquo; (24.0%) potential zones.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec24\" class=\"Section2\"\u003e \u003ch2\u003e3.4 Spatial Robustness and Uncertainty\u003c/h2\u003e \u003cp\u003eThe integration of spatial diagnostics and bootstrap resampling allowed for a rigorous evaluation of model reliability beyond standard error metrics.\u003c/p\u003e \u003cdiv id=\"Sec25\" class=\"Section3\"\u003e \u003ch2\u003e3.4.1 Spatial Independence\u003c/h2\u003e \u003cp\u003eTo confirm the mitigation of spatial autocorrelation bias, the distribution of spatial cross-validation blocks was visualised (Fig.\u0026nbsp;\u003cspan refid=\"Fig4\" class=\"InternalRef\"\u003e4\u003c/span\u003e). Assigning ground truth points to geographically distinct 2 x 2 blocks ensures that the training and validation subsets are physically separated by a sufficient distance to minimise information leakage. As indicated by the non-significant Global Moran\u0026rsquo;s \u003cem\u003eI\u003c/em\u003e (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.218) and Geary\u0026rsquo;s \u003cem\u003eC\u003c/em\u003e (\u003cem\u003ep\u003c/em\u003e\u0026thinsp;=\u0026thinsp;0.106) with values 0.036 and 0.899 for the test set residuals, this blocking strategy successfully compelled the model to learn regional hydrogeological drivers rather than local spatial dependencies. Consequently, the reported predictive accuracy is a robust measure of the model's ability to generalize to unseen terrains within the crystalline basement of the Khatra region.\u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec26\" class=\"Section3\"\u003e \u003ch2\u003e3.4.2 Model Stability and Confidence\u003c/h2\u003e \u003cp\u003eThe deterministic consensus of the optimal models was interrogated through a 60-iteration bootstrap analysis (Efron \u0026amp; Tibshirani, \u003cspan citationid=\"CR11\" class=\"CitationRef\"\u003e1994\u003c/span\u003e). A clear trade-off emerged between predictive performance and algorithmic stability. Regarding algorithmic mechanics, SVM proved the most deterministic (Bootstrap Uncertainty STD\u0026thinsp;=\u0026thinsp;0.119), a result attributed to its global convex optimization which guarantees a unique solution (Vapnik, \u003cspan citationid=\"CR46\" class=\"CitationRef\"\u003e2013\u003c/span\u003e). This was closely followed by the highly stable tree-based ensembles: Standard XGBoost (STD\u0026thinsp;=\u0026thinsp;0.129), Random Forest (STD\u0026thinsp;=\u0026thinsp;0.135), and the spatially-augmented XGBoost-SF (STD\u0026thinsp;=\u0026thinsp;0.141). CatBoost demonstrated moderate stability (STD\u0026thinsp;=\u0026thinsp;0.197). In contrast, the MLP-DL exhibited the highest instability (STD\u0026thinsp;=\u0026thinsp;0.340) due to its non-convex loss landscape (Paszke et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2019\u003c/span\u003e), manifesting as high-frequency noise and stochastic variance in the spatial outputs. This suggests that for smaller geospatial datasets (\u0026lt;\u0026thinsp;10,000 points), the complexity of deep learning may introduce unnecessary stochasticity compared to gradient boosting or kernel methods. Ultimately, this algorithmic stability metric serves as a crucial prerequisite; however, the final model selection must also satisfy rigorous independent testing discussed underneath.\u003c/p\u003e \u003c/div\u003e \u003c/div\u003e \u003cdiv id=\"Sec27\" class=\"Section2\"\u003e \u003ch2\u003e3.5 Model Validation and Selection\u003c/h2\u003e \u003cp\u003eTo evaluate true hydrogeological applicability, the six models were tested against an independent external dataset using AUC-ROC for binary discrimination and Spearman's Rank Correlation (\u0026#120588;) for ordinal fidelity. While all models exhibited exceptional binary discrimination (Fig.\u0026nbsp;\u003cspan refid=\"Fig11\" class=\"InternalRef\"\u003e11\u003c/span\u003e)\u0026mdash;with CatBoost achieving the highest AUC (0.9997), followed by MLP-DL (0.9961), XGBoost-SF (0.9919), SVM (0.9919), Standard XGBoost (0.9892), and RF (0.9856)\u0026mdash;relying solely on AUC for five-class zonation is misleading. It often fails to capture subtle, transitional 'Moderate' zones critical for regional planning, excelling primarily at separating absolute extremes ('Very High' vs. 'Very Low'). Conversely, Spearman\u0026rsquo;s Rho (\u0026#120588;) revealed the true predictive hierarchy (Fig.\u0026nbsp;\u003cspan refid=\"Fig12\" class=\"InternalRef\"\u003e12\u003c/span\u003e). Standard XGBoost achieved the highest rank fidelity (\u0026#120588; = 0.8272), followed closely by CatBoost (\u0026#120588; = 0.8244) and XGBoost-SF (\u0026#120588; = 0.8194). MLP-DL (\u0026#120588; = 0.7596) and SVM (\u0026#120588; = 0.7071) showed moderate performance, while Random Forest recorded a significant drop (\u0026#120588; = 0.6014).\u003c/p\u003e \u003cp\u003eCatBoost\u0026rsquo;s near-perfect AUC but lower \u0026#120588; reflects a binary nature that struggles with intermediate classes, whereas Standard XGBoost accurately captured continuous hydrogeological gradients. Consequently, Standard XGBoost was selected as the optimal architecture. It delivered the highest ordinal fidelity, coherent spatial logic, and high bootstrap stability (STD\u0026thinsp;=\u0026thinsp;0.129), all without the computational overhead of the XGBoost-SF variant. With Standard XGBoost established for definitive point-estimates, the final critical step is quantifying the model's spatial confidence boundaries, evaluated in the upcoming section.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec28\" class=\"Section2\"\u003e \u003ch2\u003e3.6 Spatial Uncertainty Mapping\u003c/h2\u003e \u003cp\u003eWhile Standard XGBoost provided the optimal deterministic GWPZ scores, the resulting 90% Prediction Interval (PI) map (Fig.\u0026nbsp;\u003cspan refid=\"Fig13\" class=\"InternalRef\"\u003e13\u003c/span\u003e) fundamentally alters practical decision support by visually representing the model's confidence levels. Unlike standard models that provide a single baseline estimate, this approach maps the variance in predictions. The output reveals narrow PI zones that indicate high model consensus, identifying safe, low-risk targets suitable for immediate investment. Conversely, the wide PI zones highlight regions of high epistemic uncertainty\u0026mdash;areas where the model exhibits lower confidence due to conflicting features or sparse data. By explicitly mapping these confidence boundaries, stakeholders are cautioned against blind drilling; instead, these high-uncertainty areas can be designated for mandatory geophysical verification, effectively bridging the gap between predictive modelling and field-level risk management.\u003c/p\u003e \u003cp\u003e \u003c/p\u003e \u003c/div\u003e \u003cdiv id=\"Sec29\" class=\"Section2\"\u003e \u003ch2\u003e3.7 Limitations and Future Scope\u003c/h2\u003e \u003cp\u003eWhile the proposed Uncertainty-Quantified framework demonstrates high algorithmic stability and successfully mitigates spatial autocorrelation, a minor epistemic uncertainty is introduced by the decadal temporal gap between the static ground truth yield points (2009\u0026ndash;2010) and the dynamic surface predictors (2023 LULC/NDVI). Consequently, the primary contribution of this research is the demonstration of a spatially robust machine learning pipeline\u0026mdash;integrating Spatial CV with Quantile Regression\u0026mdash;rather than a contemporary operational audit of the Khatra sub-division. Future studies will focus on deploying this validated ML architecture across real-time, high-density piezometric monitoring networks to further operationalize and refine its predictive capabilities.\u003c/p\u003e \u003c/div\u003e"},{"header":"4. Conclusion","content":"\u003cp\u003eThis study established a rigorous comparative framework to evaluate the trade-offs between predictive accuracy, stability, and hydrogeological interpretability for Groundwater Potential Zonation (GWPZ) in hard-rock terrains. Moving beyond single-algorithm applications, the Standard XGBoost model emerged as the superior architecture for definitive spatial mapping. While CatBoost excelled in binary discrimination and XGBoost-SF in internal precision, Standard XGBoost achieved the optimal synthesis of performance and ordinal rank fidelity (\u0026#120588; = 0.827).\u003c/p\u003e \u003cp\u003eMethodologically, this study demonstrates that standard random validation is insufficient for heterogeneous geospatial data. Implementing Spatial Cross-Validation revealed that conventional metrics inherently overestimate model reliability. Enforcing spatial independence provides a conservative, physically grounded benchmark reflecting true exploration challenges in uncharted aquifers. Regarding algorithmic dynamics, SVMs proved to be robust baselines for conservative planning, whereas MLP-DL exhibited high stochastic variance, indicating a need for further regularization in data-scarce domains. Physically, all models consistently prioritized Geomorphology and Runoff, confirming that landform geometry and surface permeability remain the primary controls on local aquifer storage.\u003c/p\u003e \u003cp\u003eCrucially, this work advances the paradigm from deterministic mapping to risk-aware spatial planning via a novel dual-algorithm framework. By coupling Standard XGBoost for primary GWPZ estimation with Quantile LightGBM regression for spatial uncertainty mapping, this study provides stakeholders with quantifiable confidence bounds. This clearly delineates high-confidence valley fills from high-uncertainty upland transitions, offering actionable guidance for regional planners to prioritize development in 'High Potential/Low Uncertainty' zones while mandating geophysical investigations in 'High Uncertainty' areas. Furthermore, the strong influence of land-use features suggests moderate-potential pediment plains are optimal for Managed Aquifer Recharge (MAR) interventions. Future research should integrate dynamic temporal variables, such as seasonal water level fluctuations, to refine operational yield estimates. Ultimately, this study facilitates a vital shift from static 'map-making' to dynamic 'risk-management,' empowering water managers to reject risky high-yield predictions in favour of stable, verified alternatives.\u003c/p\u003e"},{"header":"Declarations","content":" \u003cp\u003e \u003cstrong\u003eEthics approval and consent to participate:\u003c/strong\u003e \u003cp\u003eNot applicable. This study does not involve human participants, human tissue, or animal subjects.\u003c/p\u003e \u003c/p\u003e \u003cp\u003e \u003cstrong\u003eConsent for publication:\u003c/strong\u003e \u003cp\u003eNot applicable. The manuscript does not contain data from any person.\u003c/p\u003e \u003c/p\u003e\u003cp\u003e \u003ch2\u003eCompeting Interests:\u003c/h2\u003e \u003cp\u003eThe authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.\u003c/p\u003e \u003c/p\u003e\u003ch2\u003eFunding:\u003c/h2\u003e \u003cp\u003eThe first author (Suraj Dule) acknowledges the financial support provided by the University Grants Commission (UGC), Government of India, in the form of a Junior Research Fellowship (Ref. No. 210510010833).\u003c/p\u003e\u003ch2\u003eAuthor Contribution\u003c/h2\u003e\u003cp\u003eSuraj Dule: Methodology, Validation, Resources, Software, Formal analysis, Investigation, Data curation, Writing \u0026ndash; original draft, Visualization. Arabinda Sharma: Conceptualization, Methodology, Resources, Supervision, Writing \u0026ndash; review \u0026amp; editing. All authors read and approved the final manuscript.\u003c/p\u003e\u003ch2\u003eAcknowledgement\u003c/h2\u003e\u003cp\u003eThe first author acknowledges the UGC for providing financial assistance through the Junior Research Fellowship. The authors gratefully acknowledge the School of Geography, Gangadhar Meher University, for providing the necessary computational infrastructure. We also thank the Geological Survey of India (GSI) and Central Ground Water Board (CGWB) for making the ancillary geospatial data publicly available.\u003c/p\u003e\u003ch2\u003eData Availability\u003c/h2\u003e\u003cp\u003eThe raw spatial, geological, and hydro-meteorological datasets analysed in this study are publicly available from their respective institutional repositories (e.g., Geological Survey of India, USGS EarthExplorer, Climate Hazards Center, and India WRIS) as detailed in Table 1. The processed geospatial matrices, extracted validation point data, and generated machine learning model artifacts supporting the findings of this study are available from the corresponding author upon reasonable request.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\u003cli\u003e\u003cspan\u003eAdiat, K. A. N., Nawawi, M. N. M., \u0026amp; Abdullah, K. (2012). Assessing the accuracy of GIS-based elementary multi criteria decision analysis as a spatial prediction tool\u0026ndash;a case of predicting potential zones of sustainable groundwater resources. \u003cem\u003eJournal of hydrology\u003c/em\u003e, \u003cem\u003e440\u003c/em\u003e, 75\u0026ndash;89.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eAkiba, T., Sano, S., Yanase, T., Ohta, T., \u0026amp; Koyama, M. (2019, July). Optuna: A next-generation hyperparameter optimization framework. In \u003cem\u003eProceedings of the 25th ACM SIGKDD international conference on knowledge discovery \u0026amp; data mining\u003c/em\u003e (pp. 2623\u0026ndash;2631).\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBreiman, L. (2001). Random forests. \u003cem\u003eMachine learning\u003c/em\u003e, \u003cem\u003e45\u003c/em\u003e(1), 5\u0026ndash;32.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eBrenning, A. (2012, July). Spatial cross-validation and bootstrap for the assessment of prediction rules in remote sensing: The R package sperrorest. In \u003cem\u003e2012 IEEE international geoscience and remote sensing symposium\u003c/em\u003e (pp. 5372\u0026ndash;5375). IEEE.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eCentral Ground Water Board. (2022). \u003cem\u003eReport on national aquifer mapping and management plan of Bankura district, West Bengal\u003c/em\u003e. Department of Water Resources, River Development and Ganga Rejuvenation, Government of India. \u003cspan class=\"ExternalRef\"\u003e\u003cspan class=\"RefSource\"\u003ehttps://cgwb.gov.in/cgwbpnm/download/897\u003c/span\u003e\u003cspan address=\"https://cgwb.gov.in/cgwbpnm/download/897\" targettype=\"URL\" class=\"RefTarget\"\u003e\u003c/span\u003e\u003c/span\u003e\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChen, T. (2016). XGBoost: A Scalable Tree Boosting System. \u003cem\u003eCornell University\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eChowdary, V. M., Chakraborthy, D., \u0026amp; Jeyaram, A. (2009). Multi-criteria decision-making approach for watershed management using remote sensing and GIS. \u003cem\u003eWater Resources Management\u003c/em\u003e, 23(8), 1553\u0026ndash;1573.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDikshit, A., Pradhan, B., \u0026amp; Alamri, A. M. (2021). Pathways and challenges of the application of artificial intelligence to geohazards modelling. \u003cem\u003eGondwana Research\u003c/em\u003e, \u003cem\u003e100\u003c/em\u003e, 290\u0026ndash;301.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDrucker, H., Burges, C. J., Kaufman, L., Smola, A., \u0026amp; Vapnik, V. (1996). Support vector regression machines. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e9\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eDule, S., \u0026amp; Sharma, A. (2025). Deciphering groundwater potential zones of Sarenga block, West Bengal, India using Geographic Information System and Multi-Criteria Decision-Analysis methods. \u003cem\u003eInternational Journal of Engineering and Geosciences\u003c/em\u003e, \u003cem\u003e11\u003c/em\u003e(2), 301\u0026ndash;320.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eEfron, B., \u0026amp; Tibshirani, R. J. (1994). \u003cem\u003eAn introduction to the bootstrap\u003c/em\u003e. Chapman and Hall/CRC.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFisher, A., Rudin, C., \u0026amp; Dominici, F. (2019). All models are wrong, but many are useful: Learning a variable's importance by studying an entire class of prediction models simultaneously. \u003cem\u003eJournal of Machine Learning Research\u003c/em\u003e, \u003cem\u003e20\u003c/em\u003e(177), 1\u0026ndash;81.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eFriedman, J. H. (2001). Greedy function approximation: a gradient boosting machine. \u003cem\u003eAnnals of statistics\u003c/em\u003e, 1189\u0026ndash;1232.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGeary, R. C. (1954). The contiguity ratio and statistical mapping. \u003cem\u003eThe incorporated statistician\u003c/em\u003e, \u003cem\u003e5\u003c/em\u003e(3), 115\u0026ndash;146.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eG\u0026eacute;ron, A. (2022). \u003cem\u003eHands-on machine learning with Scikit-Learn, Keras, and TensorFlow\u003c/em\u003e. \" O'Reilly Media, Inc.\".\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eGuo, C., \u0026amp; Berkhahn, F. (2016). Entity embeddings of categorical variables. \u003cem\u003earXiv preprint arXiv:1604.06737\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHastie, T., Tibshirani, R., \u0026amp; Friedman, J. (2009). The elements of statistical learning: Data mining, inference, and prediction (2nd ed.). Springer.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eHillel, D. (2013). \u003cem\u003eIntroduction to soil physics\u003c/em\u003e. Academic press.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eJenks, G. F. (1967). The Data Model Concept in Statistical Mapping. \u003cem\u003eInternational Yearbook of Cartography\u003c/em\u003e, \u003cem\u003e7\u003c/em\u003e, 186\u0026ndash;190.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKar, S., Barman, B., \u0026amp; Das, J. (2025). Identification of Groundwater Potential Zones in North-Western Bankura, West Bengal, India, Using GIS Based MCDA Techniques TOPSIS and CoCoSo. \u003cem\u003eJournal of the Indian Society of Remote Sensing\u003c/em\u003e, 1\u0026ndash;17.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKe, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., \u0026hellip; Liu, T. Y. (2017). Lightgbm: A highly efficient gradient boosting decision tree. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e30\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eKoenker, R., \u0026amp; Bassett Jr, G. (1978). Regression quantiles. \u003cem\u003eEconometrica: journal of the Econometric Society\u003c/em\u003e, 33\u0026ndash;50.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLiu, R., Gou, J., Han, J., Cao, K., Pan, X., Wang, D., \u0026amp; Wei, L. (2025). Optimizing groundwater potential assessment: uncertainty reduction through sample balancing and enhanced hybrid modeling. \u003cem\u003eStochastic Environmental Research and Risk Assessment\u003c/em\u003e, 1\u0026ndash;21.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eLundberg, S. M., \u0026amp; Lee, S. I. (2017). A unified approach to interpreting model predictions. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e30\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMachiwal, D., Jha, M. K., \u0026amp; Mal, B. C. (2011). Assessment of groundwater potential in a semi-arid region of India using remote sensing, GIS and MCDM techniques. \u003cem\u003eWater resources management\u003c/em\u003e, \u003cem\u003e25\u003c/em\u003e(5), 1359\u0026ndash;1386.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMogaji, K. A., Lim, H. S., \u0026amp; Abdullah, K. (2015). Regional prediction of groundwater potential using GIS-based Dempster\u0026ndash;Shafer theory and frequency ratio models. \u003cem\u003eEnvironmental Earth Sciences\u003c/em\u003e, 73(1), 327\u0026ndash;346.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMoran, P. A. (1950). Notes on continuous stochastic phenomena. \u003cem\u003eBiometrika\u003c/em\u003e, \u003cem\u003e37\u003c/em\u003e(1/2), 17\u0026ndash;23.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eMussa, M. M., Lohani, T. K., \u0026amp; Eshete, A. A. (2024). Evaluation of Groundwater Potential Zones Using GIS-Based Machine Learning Ensemble Models in the Gidabo Watershed, Ethiopia. \u003cem\u003eGlobal Challenges\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(12), 2400137.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eNaghibi, S. A., Pourghasemi, H. R., \u0026amp; Dixon, B. (2016). GIS-based groundwater potential mapping using boosted regression tree, classification and regression tree, and random forest machine learning models in Iran. \u003cem\u003eEnvironmental monitoring and assessment\u003c/em\u003e, \u003cem\u003e188\u003c/em\u003e(1), 44.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePanahi, M., Sadhasivam, N., Pourghasemi, H. R., Rezaie, F., \u0026amp; Lee, S. (2020). Spatial prediction of groundwater potential mapping based on convolutional neural network (CNN) and support vector regression (SVR). \u003cem\u003eJournal of Hydrology\u003c/em\u003e, \u003cem\u003e588\u003c/em\u003e, 125033.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePaszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., \u0026hellip; Chintala, S. (2019). Pytorch: An imperative style, high-performance deep learning library. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e32\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003ePrasad, R. K., Mondal, N. C., Banerjee, P., Nandakumar, M. V., \u0026amp; Singh, V. S. (2008). Deciphering potential groundwater zone in hard rock through the application of GIS. \u003cem\u003eEnvironmental geology\u003c/em\u003e, \u003cem\u003e55\u003c/em\u003e(3), 467\u0026ndash;475.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eProkhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., \u0026amp; Gulin, A. (2018). CatBoost: unbiased boosting with categorical features. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e31\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRana, M. S. P., Rahman, M. T., \u0026amp; Hassan, M. F. (2025). Mapping groundwater potential zone by robust machine learning algorithms \u0026amp; remote sensing techniques in agriculture dominated area, Bangladesh. \u003cem\u003eCleaner Water\u003c/em\u003e, \u003cem\u003e3\u003c/em\u003e, 100064.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRay, S. K. (2025). Unveiling groundwater gems: A GIS-powered fusion of AHP and TOPSIS for mapping groundwater potential zones. \u003cem\u003eGroundwater for Sustainable Development\u003c/em\u003e, \u003cem\u003e29\u003c/em\u003e, 101431.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eRey, S. J., \u0026amp; Anselin, L. (2007). PySAL: A Python library of spatial analytical methods. \u003cem\u003eThe Review of Regional Studies\u003c/em\u003e, 37(1), 5\u0026ndash;27.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSaha, A., \u0026amp; Pal, S. C. (2025). Delineation of groundwater potential zones using data-driven approaches: towards achieving sustainable groundwater management in drought-prone region of Eastern India. \u003cem\u003eEnvironmental Monitoring and Assessment\u003c/em\u003e, \u003cem\u003e197\u003c/em\u003e(10), 1\u0026ndash;22.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSaha, R., Chiravuri, S. S., Das, I. C., Kandrika, S., Kumranchat, V. K., Chauhan, P., \u0026amp; Chitikela, V. L. (2024). Urban aquifer health assessment and its management for sustainable water supply: an innovative approach using machine learning techniques. \u003cem\u003eGroundwater for Sustainable Development\u003c/em\u003e, \u003cem\u003e25\u003c/em\u003e, 101130.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSaltelli, A., Annoni, P., Azzini, I., Campolongo, F., Ratto, M., \u0026amp; Tarantola, S. (2010). Variance based sensitivity analysis of model output. Design and estimator for the total sensitivity index. \u003cem\u003eComputer physics communications\u003c/em\u003e, \u003cem\u003e181\u003c/em\u003e(2), 259\u0026ndash;270.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSander, P. (2007). Lineaments in groundwater exploration: a review of applications and limitations. \u003cem\u003eHydrogeology journal\u003c/em\u003e, \u003cem\u003e15\u003c/em\u003e(1), 71\u0026ndash;74.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eShaban, A., Khawlie, M., \u0026amp; Abdallah, C. (2006). Use of remote sensing and GIS to determine recharge potential zones: the case of Occidental Lebanon. \u003cem\u003eHydrogeology Journal\u003c/em\u003e, \u003cem\u003e14\u003c/em\u003e(4), 433\u0026ndash;443.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSinghal, B. B. S., \u0026amp; Gupta, R. P. (2010). \u003cem\u003eApplied hydrogeology of fractured rocks\u003c/em\u003e. Springer Science \u0026amp; Business Media.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eSoboĺ, I. M. (1993). Sensitivity estimates for nonlinear mathematical models. \u003cem\u003eMath. Model. Comput. Exp.\u003c/em\u003e, \u003cem\u003e1\u003c/em\u003e, 407.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTobler, W. R. (1970). A computer movie simulating urban growth in the Detroit region. \u003cem\u003eEconomic geography\u003c/em\u003e, \u003cem\u003e46\u003c/em\u003e(sup1), 234\u0026ndash;240.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eTodd, D. K., \u0026amp; Mays, L. W. (2004). \u003cem\u003eGroundwater hydrology\u003c/em\u003e. John Wiley \u0026amp; Sons.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVapnik, V. (2013). \u003cem\u003eThe nature of statistical learning theory\u003c/em\u003e. Springer science \u0026amp; business media.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVaswani, A., Shazeer, N., Parmar, N., Uszkoreit, J., Jones, L., Gomez, A. N., \u0026hellip; Polosukhin, I. (2017). Attention is all you need. \u003cem\u003eAdvances in neural information processing systems\u003c/em\u003e, \u003cem\u003e30\u003c/em\u003e.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eVishwakarma, B. D., Zhang, J., \u0026amp; Sneeuw, N. (2021). Downscaling GRACE total water storage change using partial least squares regression. \u003cem\u003eScientific data\u003c/em\u003e, \u003cem\u003e8\u003c/em\u003e(1), 95.\u003c/span\u003e\u003c/li\u003e \u003cli\u003e\u003cspan\u003eYeh, H. F., Cheng, Y. S., Lin, H. I., \u0026amp; Lee, C. H. (2016). Mapping groundwater recharge potential zone using a GIS approach in Hualian River, Taiwan. \u003cem\u003eSustainable Environment Research\u003c/em\u003e, \u003cem\u003e26\u003c/em\u003e(1), 33\u0026ndash;43.\u003c/span\u003e\u003c/li\u003e\u003c/ol\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Groundwater Potential Mapping, Machine Learning, Deep Learning, Spatial Cross-Validation, Spatially Explicit Uncertainty, Hard-Rock Aquifers","lastPublishedDoi":"10.21203/rs.3.rs-9249309/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-9249309/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eSustainable groundwater management in complex hard-rock terrains requires the precise delineation of potential zones. Traditional learning approaches used for groundwater potential mapping often fail to account for the spatial dependence of hydrogeological data and the inherent uncertainty in predictions. This study presents a rigorous comparative framework for Groundwater Potential Zonation (GWPZ) evaluating six architectures: Random Forest, CatBoost, Support Vector Machines (SVM), standard XGBoost, spatially augmented XGBoost (XGB-SF), and Deep Learning Multi-Layer Perceptron (MLP-DL). To mitigate the spatial autocorrelation bias common in conventional ML/DL applications, the study proposes an innovative spatial cross-validation strategy integrated within an Optuna hyperparameter tuning workflow. Results indicated that while the XGB-SF model achieved the highest internal accuracy (\u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.853), the standard XGBoost model proved optimal for operational mapping. It delivered the most suitable balance of predictive accuracy (\u003cem\u003eR\u003c/em\u003e\u003csup\u003e\u003cem\u003e2\u003c/em\u003e\u003c/sup\u003e\u0026thinsp;=\u0026thinsp;0.846), algorithmic stability (Bootstrap STD\u0026thinsp;=\u0026thinsp;0.129), and ordinal ranking fidelity for groundwater potential (Spearman\u0026rsquo;s \u0026#120588; = 0.827). Furthermore, feature importance analysis revealed that tree-based ensembles successfully navigated severe multicollinearity without manual feature elimination, consistently identifying geomorphology and surface runoff as the dominant physical controls on hard-rock aquifer storage. Conversely, SVM demonstrated exceptional stability with zero overfitting, serving as a robust baseline for data-scarce regions. Crucially, this research advances decision support by generating spatially explicit uncertainty maps via quantile regression. This enables planners to distinguish high-confidence zones suitable for immediate exploitation from high-uncertainty areas requiring geophysical verification, providing a scientifically defensible framework for exploration in heterogeneous aquifer systems.\u003c/p\u003e","manuscriptTitle":"Modelling Groundwater Potential of Tropical Hard-Rock Aquifers using Spatial Cross-Validation augmented with Uncertainty and Sensitivity Analysis","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2026-04-20 08:38:32","doi":"10.21203/rs.3.rs-9249309/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"5a4a181d-90f3-41fe-93ee-e5f6501fc172","owner":[],"postedDate":"April 20th, 2026","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2026-04-20T08:38:32+00:00","versionOfRecord":[],"versionCreatedAt":"2026-04-20 08:38:32","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-9249309","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-9249309","identity":"rs-9249309","version":["v1"]},"buildId":"XKTyCvWXoU3ODBz1xrDgd","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}