Unraveling Groundwater Stress: A Model Generalizable Across Space and Time Reveals Co-Declining Water Tables and Recharge Across Bangladesh

Saleh Sakib Ahmed; Rashed Uz Zzaman; Saifur Rahman Jony; Faizur Rahman Himel; Afroza Sharmin; A. H. M. Khalequr Rahman; M. Sohel Rahman; Sara Nowreen

Figure 1: Methodology: a. Phase 1: Teacher Model predicting maximum groundwater level (bgl) using 17 HGFs and year as inputs. b. Phase 2: Teacher Model predicting minimum groundwater level (bgl) conditioned on maximum groundwater level (bgl) to ensure consistency, with inputs including 17 HGFs, year, and maximum groundwater level (bgl). c. Upsampling Student Model combining grid-level low-resolution information (annual maximum and minimum groundwater storage along with 17 HGFs of grid center-point) and point-level high-resolution data (17 HGFs of 2 km), merged by GLDAS grid ID and year, to predict high-resolution annual maximum and minimum groundwater level (bgl). The target for the Student is the Teacher’s pseudo-labels and the original in-situ data

A. Year-Specific Architecture Comparison (Ablation Study)
Category / Experiment	Model / Framework	$\mathbf{R^{2}}$	MSE	RMSE	MAE	MAPE (%)
Max GWL	Random Forest	0.3115	22.8765	4.7068	2.3979	34.20
	Neural Network	0.1993	26.5521	5.0912	2.7973	41.13
	XGBoost	0.2172	25.9802	5.0143	2.6763	39.29
	LightGBM	0.2688	24.3531	4.8523	2.6511	38.78
Min GWL	Random Forest	0.1082	33.1914	5.6804	3.0257	82.76

Table 1: Comprehensive Performance Comparison of Teacher Model Candidates for Task 1 (Max and Min GWL Prediction). Section A evaluates model performance under year-specific training constraints. Section B compares combined-year training architectures, identifying Random Forest as the top-performing ML backbone. Section C benchmarks the proposed Groundhog framework (Combined-Year Random Forest) against traditional spatial interpolation baselines.

Figure 2: Performance Metrics: a. Radar plot for different metrics of IDW, GroundHog and Year-Specific Random Forest. b.

R^{2}

scores of IDW, GroundHog and Year-Specific Random Forest. c. Leave-One-Year-Out (LOYO) validation for the Upsampling Model.

Target Group	Metric	Random Forest	Neural Network	XGBoost	LightGBM
Overall	$R^{2}$	0.9344	0.7125	0.8355	0.7487
	MSE	0.6466	2.9228	1.6891	2.6023
	RMSE	0.8041	1.7096	1.2997	1.6132
	MAE	0.3769	1.0707	0.7714	0.9432
Min GWL	$R^{2}$	0.9191	0.6803	0.8233	0.7383
	MSE	0.6322	2.4991	1.3812	2.0452

Table 2: Performance Comparison of Task 2 Student Model Candidates Across Multi-Scale Upsampling Targets. Bold values indicate top-performing models, while underlined values denote second-best performers.

Figure 3: Geospatial model validation, epistemic uncertainty, and Monte Carlo trend robustness suite for groundwater recharge across Bangladesh. (A) District-level Leave-One-District-Out (LODO) cross-validation error, quantified via Mean Absolute Percentage Error ( MAPE %), highlighting regions where localized anthropogenic extraction or unique physical boundary conditions limit out-of-fold generalization. (B) Random Forest ensemble model uncertainty, measured as the mean standard deviation across individual decision trees (

m

). (C) Binary classification of point-wise trend robustness derived from

1,000

Monte Carlo simulation iterations, isolating statistically robust trends (deep teal; empirical

95\%

confidence interval strictly excluding zero) from non-robust signals (soft coral; error envelope spanning zero). (D) Point-wise mean annual groundwater recharge trend slopes ( mm/year ), using a symmetric diverging colormap to differentiate localized recharge accretion (blue) from temporal depletion (red).

Figure 4: Comparison of sparse observations, IDW interpolation and pseudo-ground-truth groundwater-level maps for 2008. Maps show depth below ground level (m), with maximum and minimum groundwater levels shown in the bottom and top rows, respectively. (a,d) Sparse in-situ observations; (b,e) IDW interpolation, which produces physically inconsistent regions where minimum depth exceeds maximum depth (circled); and (c,f) GroundHog Teacher predictions, which provide continuous spatial fields while preserving the physical ordering of groundwater-level extremes.

Figure 5: Task 2: GLDAS groundwater storage and upsampled groundwater levels for 2008. (a,c) GLDAS minimum and maximum groundwater storage at 25 km resolution (mm); (b,d) corresponding upsampled annual minimum and maximum groundwater levels (m below ground level).

Figure 6: Spatiotemporal changes in national-average groundwater recharge and water-table depth, 2003–2024. (A) Linear trend in national-average annual recharge, showing a statistically significant decline of 0.110 cm yr -1 (

R^{2}=0.654

,

p<0.001

); the 95% confidence interval is shaded. (B) Distribution of annual recharge across 5-year epochs; boxes show interquartile ranges and red points show means

\pm

1 standard deviation. (C) Trajectory of national-average annual recharge (mm) against seasonal groundwater-level extremes (m bgl), with points colored by year. Regression bounds and dimension arrows illustrate a differential deepening effect: the annual minimum (shallower) groundwater level deepens at a faster rate than the maximum, leading to a progressively shrinking seasonal fluctuation gap. (D) Latitudinal variation in groundwater depth through time.

Figure 7: Temporal comparison of upsampled results: Spatial distribution of upsampled values in 2003 (top row) and 2022 (bottom row) for a,b. annual minimum groundwater depth (m bgl), c,d. annual maximum groundwater depth (m bgl), and e,f. annual recharge (cm).

Figure 8: Spatial trends in groundwater levels and recharge: Sen’s slope (top row; mm yr -1 for recharge, cm yr -1 for groundwater depth) and autocorrelation-adjusted Modified Mann–Kendall significance (

p^{*}<0.05

; bottom row) for a,d. minimum groundwater depth, b,e. maximum groundwater depth, and c,f. recharge. For recharge, Sen’s slopes are grouped into six categories: negligible trends (

-0.5

to

0.5

mm yr -1 ; light grey), moderate decline (

-3

to

-0.5

mm yr -1 ; yellow), strong decline (

-5

to

-3

mm yr -1 ; orange), severe decline (

-10

to

-5

mm yr -1 ; dark red), and positive accretion (

0.5

to

10

mm yr -1 ; blue shades).

Figure S2 : Geographic reference map of selected Bangladesh districts evaluated during Leave-One-District-Out (LODO) cross-validation.

Figure S3 : Distribution of key Agro-Ecological Zones (AEZs) and major river networks across Bangladesh.

Figure S4 : Effects of Variable Changes on Groundwater Level: a. Partial Dependence Plot (PDP) for the maximum groundwater level regressor illustrates feature impacts while holding others constant. b. Feature correlation heatmap highlights key relationships. c. 3D Partial Dependence Plot of Stream Power Index (SPI) and Topographic Wetness Index (TWI).

Figure S5 : Feature importance and SHAP values: In the SHAP plot, red values indicate higher feature values, while blue values indicate lower feature values. The X-axis represents SHAP values, showing positive (high) or negative (low) influence on the output. a, c show feature importance for the Max and Min GWL Models, respectively. b, d present SHAP values for the Max and Min GWL Models, respectively.

Figure S6 : Feature Importance of Student Upsampling model

Figure S7 : SHAP score of the Max GWL output of Student Upsampling model

Figure S8 : SHAP score of the Min GWL output of Student Upsampling model

Figure S9 : Future Prediction for 2024: Leveraging our Upsampling Model, we have upsampled the 2024 GLDAS groundwater storage data to obtain high-resolution groundwater level estimates even before in-situ measurements become available. The upsampled data is shown for sites where various organizations in Bangladesh monitor groundwater levels.

Figure S10 : Spatial trend significance using standard Mann–Kendall testing: Conventional Mann–Kendall trend significance across Bangladesh for a. minimum groundwater depth, b. maximum groundwater depth, and c. groundwater recharge. Categories indicate statistically significant negative trends (red), significant positive trends (blue), and non-significant trends (grey).

Figure S11 : Data Distribution: A, B. Annual maximum and minimum groundwater levels. C. Combined multi-year data showing maximum depths exceed minimums (below ground level).

Figure S12 : Reference Map: This map highlights various locations in Bangladesh that are mentioned in the manuscript.

Figure S13 : Hydro-Geological Feature Map: Distribution of 17 hydro-geological feature (HGF) values across the country.

Features (data type)	Characteristics	Equation/ Tool name	Data Source
Slope (Numerical)	The steepness of the terrain, derived from the DEM. Slope influences surface runoff, soil erosion, and infiltration Gao et al. (2024) .	$\text{Slope\%}=100\times\frac{\sqrt{AX^{2}+BY^{2}}}{\text{Pixel Size (DEM)}}$ Where: • $AX$ = Filtered DEM with x-gradient filter • $BY$ = Filtered DEM with y-gradient filter • Pixel Size (DEM) = Pixel size of the DEM Ritter (2007) Note: Derived from slope of Spatial Analyst Tool in ArcGIS Pro platform	Derived from Shuttle Radar Topography Mission (SRTM) NASA Earthdata (2025) based digital elevation model (SRTM-DEM) of 30 m resolution
Drainage Density (Numerical)	The total length of streams per unit area of the watershed. Higher drainage density indicates better drainage and less water retention in the soil Anderson and Anderson (2010) .	$D_{d}=\frac{\sum_{i=1}^{n}D_{i}}{A}$ Where: • $D_{i}$ : Total length of streams • $A$ : Area of a grid • $D_{d}$ : Drainage Density Horton (1932) Note: Derived from flow accumulation, stream order, fill, etc. of Spatial Analyst Tool in ArcGIS Pro platform	Derived from SRTM-DEM of 30 m resolution
Elevation (Numerical)	The height above sea level, derived from the DEM. Elevation influences climatic and hydrological conditions Kenward et al. (2000) .		Derived from SRTM-DEM of 30 m resolution
Distance from Stream (Numerical)	The proximity of a location to the nearest stream, which affects the availability of water and the potential for flooding Suwanno et al. (2023) .	$d_{ij}=\sqrt{\sum_{k=1}^{n}(X_{ik}-Y_{jk})^{2}}$ Where: • $d_{ij}$ : Distance from stream for $i,j$ locations • $k$ : Represents features Note: Used Euclidean distance of Spatial Analyst Tool to calculate in ArcGIS Pro platform.	Derived from SRTM-DEM of 30 m resolution
Topographic Wetness Index (TWI) (Numerical)	A measure of topographic control on hydrological processes, calculated using slope and upstream contributing area. Higher values indicate wetter areas.	$TWI=\ln\left(\frac{A_{s}}{\beta}\right)$ Where: • $A_{s}$ : Catchment area per unit • $\beta$ : Slope angle in degrees • $TWI$ : Topographic Wetness Index Note: Calculated from flow accumulation and slope using the Spatial Analyst Tool in ArcGIS Pro platform.
Terrain Ruggedness Index (TRI) (Numerical)	Quantifies the terrain’s ruggedness by measuring elevation differences between adjacent grid cells. Rugged terrain affects runoff and groundwater recharge Grinevskii (2014) ; Verma and Mirajkar (2024) .	$R=\sqrt{\text{max}^{2}-\text{min}^{2}}$ Where: • max : The highest cell value in the DEM • min : The lowest cell value in the DEM Note: Calculated using ArcGIS Pro platform.	Derived from SRTM-DEM of 30 m resolution

Table S2 : Summary of Numerical and Categorical Factors: First 17 features are referred to as Hydro-Geological Features (HGFs) in the text.

Unraveling Groundwater Stress: A Model Generalizable Across Space and Time Reveals Co-Declining Water Tables and Recharge Across Bangladesh

Abstract

Figures & tables

Hydrology

Pure and physics-guided deep learning approaches for spatio-temporal groundwater level prediction

Unsupervised Detection of Groundwater Storage Anomalies in Ghana Using GRACE Satellite Data

Reconstructing GRACE Terrestrial Water Storage with Spatio-Temporal Graph Neural Networks: An Application to South America

Abstract

Hydrology

Explore similar work

Pure and physics-guided deep learning approaches for spatio-temporal groundwater level prediction

Unsupervised Detection of Groundwater Storage Anomalies in Ghana Using GRACE Satellite Data

Reconstructing GRACE Terrestrial Water Storage with Spatio-Temporal Graph Neural Networks: An Application to South America