Inferring groundwater overdraft in data-scarce arid agro-ecosystems: a semi-empirical remote sensing framework applied to the Elfeija watershed, Morocco

Abstract

Introduction:

Small-scale irrigated agriculture in arid regions relies heavily on unmetered groundwater, creating an “invisible pumping” threat to aquifer sustainability. Monitoring these diffuse withdrawals remains a critical challenge for water governance. This study proposes a synergistic multi-sensor remote sensing workflow to infer groundwater abstraction and assess its impact on local water stress using a semi-empirical thermal-phenological model.

Methods:

We leveraged high-resolution vegetation phenology (NDVI) from Sentinel-2 and thermal data (LST) from Landsat 8/9 over the Elfeija watershed, Morocco (2020–2024). By isolating the thermal-phenological anomaly (ΔLST) between irrigated plots and natural reference areas, we translated the surface cooling effect into evapotranspiration fluxes and pumping volumes. To ensure physical plausibility and avoid circularity, the approach was constrained by independent bottom-up agronomic benchmarks alongside a multi-tier cross-comparison against global models (MOD16, ERA5-Land, GLEAM), and a rigorous Monte Carlo uncertainty propagation integrating both parametric and recharge uncertainties.

Results:

The analysis reveals a distinct seasonal signature of irrigation. Probabilistic modeling indicates that the inferred annual pumping for 2022 (central estimates: 5.85–6.79 Mm³) strongly suggests a state of severe overdraft, with a high probability (P > 80%) of exceeding the estimated renewable aquifer recharge (4.70 Mm³). Furthermore, the framework successfully disentangled climatic triggers from anthropogenic forcing, capturing an anomalous, water-intensive late-season agricultural cycle in 2024.

Discussion:

While the lack of in-situ piezometric data limits absolute ground-truthing, the multi-source convergence (global models and local infrastructure data) robustly brackets the overdraft scenario. The proposed semi-empirical inference workflow provides a scalable, cost-effective framework for data-scarce regions, offering river basin agencies a proactive tool to identify over-extraction hotspots and negotiate sustainable water quotas.

1 Introduction

Water security in arid and semi-arid regions represents one of the most pressing challenges of the 21st century, exacerbated by the compounding pressures of climate change and growing food demand. In these regions, groundwater has become a critical buffer resource, yet its intensive use is leading to widespread depletion, threatening both agricultural sustainability and ecosystem health globally (Bierkens and Wada, 2019; Scanlon et al., 2023). Agriculture is the dominant consumer, with its expansion into water-scarce areas driving an unprecedented and often unregulated reliance on aquifers (Du et al., 2024; Rosa, 2022). This challenge is compounded by weak governance and inadequate monitoring, particularly in regions such as South Asia Ricciardi et al. (2020), the Middle East, and North Africa (Hasan and Fatima, 2025; Ahmad and Al-Ghouti, 2020). The core of the management problem lies in the diffuse and “invisible” nature of this abstraction, where thousands of individual unmetered small-scale farmers contribute collectively to a significant, yet largely unquantified, regional water stress.

Morocco serves as a representative case study of this global challenge. As a key sector of the national economy, agriculture consumes about 88% of the country’s water, with a heavy and increasing reliance on groundwater (Baba et al., 2025). Government policies, such as the Green Morocco Plan, have often encouraged the cultivation of high-value, water-intensive crops, leading to severe aquifer over-extraction (Ouhbi and Boudhar, 2025). This issue is particularly acute for crops like watermelon in the Drâa Valley, where a shift from traditional surface irrigation to intensive private well pumping has exacerbated local water stress (Bossenbroek, et al., 2023). This agricultural pressure is compounded by a steady decline in annual precipitation, further diminishing natural aquifer recharge (Baba et al., 2025). Recent studies in key agricultural basins like the Middle Atlas have documented alarming declines in groundwater levels, with drops of up to 5 m per year, leading to the complete drying of surface water bodies such as Lake Aoua (El-Bouhali et al., 2023). The proliferation of private wells, especially for small-scale farming, has been identified as a major driver of this depletion (Massuel et al., 2017).

In response to this critical data gap, satellite remote sensing has emerged as a powerful tool for monitoring agricultural water use over large areas (Massari et al., 2021). Satellites provide a means to observe the biophysical consequences of water application. While sophisticated Surface Energy Balance (SEB) models (e.g., SEBAL, METRIC) provide highly accurate actual ET estimations (Bastiaanssen et al., 1998), their operational application for continuous, multi-annual monitoring of smallholder farming is often hindered by the extensive requirement for high-resolution, in-situ micrometeorological data. Therefore, there is a critical need for simplified, proxy-based operational workflows. The Normalized Difference Vegetation Index (NDVI) tracks vegetation health and phenology (Tucker, 1979), while Land Surface Temperature (LST) serves as a proxy for the surface energy balance. The process of evapotranspiration (ET) induces a detectable cooling effect on irrigated canopies compared to dry surroundings (Eamus et al., 2015). The synergistic analysis of NDVI and LST time-series offers a pragmatic physical basis for monitoring crop water status when ground data are scarce.

Recent advances in satellite technology, particularly the complementary missions of Sentinel-2 and Landsat, have opened new frontiers for this approach. The integration of high-resolution multispectral data from Sentinel-2 (10–20 m, 5-day revisit) and the rich thermal history of Landsat (30–100 m, 8–16 days revisit) provides an unprecedented combination of spatial detail and temporal frequency (Singh et al., 2020). This synergy allows for the reduction of data gaps caused by cloud cover and improves the accuracy of ET estimations, which is crucial for effective water resource management (Saadi et al., 2015). However, while many studies have successfully used this synergy to map irrigated areas or detect water stress (Conrad et al., 2020), a critical research gap remains in translating these observations into robust, multi-annual quantifications of abstracted groundwater volumes and assessing their direct contribution to local water stress UN-Water (2023).

2 Materials and methods

2.1 Study area: the elfeija watershed, Morocco

This study focuses on the Elfeija watershed, located within the larger Drâa River Basin in southeastern Morocco (approximate center: N 30°10′00″, W 6°10′00″) (Figure 1). The region is characterized by a hyper-arid climate, with annual precipitation averaging less than 100 mm and experiencing high inter-annual variability. Potential evapotranspiration rates are extremely high, exceeding 2,000 mm per year, resulting in a significant structural water deficit that makes rainfed agriculture unviable (Amiha et al., 2024). The landscape is a typical arid agro-ecosystem, comprising sparse natural vegetation (steppe) interspersed with concentrated pockets of irrigated agriculture, primarily located along the wadi systems and in plains where groundwater is accessible.

FIGURE 1

FIGURE 2

Table 1 summarizes the key characteristics of these polygon groups. The 16“Irrigated Plots” cover a total area of 8,780 ha and contain the vast majority of the agricultural activity, with an estimated 580 active wells identified through high-resolution imagery. The 16“Reference Areas,” covering a slightly larger area of 11,126 ha, were selected to represent the surrounding natural landscape. Both groups are situated at a similar mean elevation (≈775 m) and on very low slopes (≈1.5%), indicating that the primary discriminating factor between them is land use rather than topography.

To characterize the overall hydro-climatic conditions of the watershed, we analyzed the mean monthly dynamics of NDVI and LST across the entire study area (125,175 ha) for the 2020–2024 period, as derived from Sentinel-2 and Landsat data respectively. The results are presented in Figure 3. The watershed exhibits a typical arid climate signature, with LST peaking at over 60 °C during the hot, dry summer months (June-August). A clear inverse relationship is observed between temperature and vegetation, with the basin-wide mean NDVI reaching its lowest point during these summer peaks. A slight increase in NDVI is noticeable in the cooler, wetter winter and early spring months, reflecting the opportunistic response of the sparse natural vegetation to precipitation events. The large standard deviation observed for LST, particularly during transitional seasons, highlights the significant thermal heterogeneity of the landscape. This climatological baseline provides the essential context for interpreting the specific anomalies induced by irrigated agriculture, which will be detailed in the Results section.

FIGURE 3

2.2 Satellite and ancillary data

To conduct this multi-annual analysis (1 January 2020, to 31 December 2024), a suite of multi-source satellite and climate reanalysis datasets was acquired and processed within the Google Earth Engine (GEE) platform.

2.2.1 Sentinel-2 MSI

Sentinel-2 Level-2A surface reflectance data from the COPERNICUS/S2_SR_HARMONIZED collection were used to derive vegetation information. This dataset provides multispectral imagery at 10–20 m spatial resolution with a high revisit frequency (∼5 days). The red (Band 4, 10 m) and near-infrared (NIR, Band 8, 10 m) bands were used to compute the Normalized Difference Vegetation Index (NDVI) as NDVI = (NIR - Red)/(NIR + Red). The QA60 quality band was used to mask out pixels contaminated by clouds and cirrus.

2.2.2 Landsat 8/9 OLI/TIRS

Land Surface Temperature (LST) data were derived from the Landsat eight and 9 Collection2 Level-2 Science Products (LANDSAT/LC08/C02/T1_L2 and LANDSAT/LC09/C02/T1_L2). The combined constellation provides thermal imagery at an 8-day revisit frequency. We used the Surface Temperature band (ST_B10), to which the appropriate scale factor (multiply by 0.00341802 and add 149.0) was applied to convert digital numbers to Kelvin, and then to Celsius. The QA_PIXEL band was used to mask pixels affected by clouds and cloud shadows, ensuring data quality.

2.2.3 Climate reanalysis data (ERA5-Land)

To provide a climatic baseline for evapotranspiration, we used the ECMWF/ERA5_LAND/MONTHLY_AGGR dataset. This product provides monthly aggregated data at a ∼9 km resolution. We used the total_evaporation_sum variable, which represents the accumulated monthly ET. Following the product’s convention, a factor of −1000 was applied to convert the data from meters of water depth (with a negative sign for an outgoing flux) to a positive value in millimeters per month (mm/month).

2.2.4 Local hydrogeological data

To assess the sustainability of the estimated pumping volumes, we utilized ancillary data from local literature for the year 2022 (Amiha et al., 2025; Agence du Bassin Hydraulique de Guelmim-Oued Noun, 2014). These data include an estimated annual lateral groundwater recharge of 2.72 Mm³ and a direct infiltration recharge of 1.98 Mm³, for a total annual renewable resource of 4.70 Mm³.

2.3 Methodology

To address the need for continuous monitoring, we developed a replicable, three-step framework designed to translate raw satellite observations into quantified water volumes and sustainability metrics. All spatial processing and time-series extractions were performed within the Google Earth Engine (GEE) platform. The 5-year timeframe (2020–2024) was specifically selected to ensure dense data availability from the combined Sentinel-2 and Landsat 8/9 constellations, while capturing significant inter-annual climatic variability, including a major drought and an exceptional late-season anomaly. The workflow is summarized in Figure 4.

FIGURE 4

2.3.1 Time-series extraction

2.3.2 Estimation of irrigation water volume

The surplus evapotranspiration from irrigation (ET_irrigation) was estimated by translating the observed thermal anomaly (ΔLST) into a water flux. This approach is conceptually rooted in the surface energy balance, where the difference in available energy between a dry and a wet surface is primarily dissipated as latent heat flux (i.e., ET). Numerous studies have demonstrated a strong, often linear, inverse relationship between ET and LST under energy-limited conditions (Jackson et al., 1981; Carlson et al., 1994). Building on this principle, we modeled the relationship as Equations 1, 1a:

With

The coefficient c can be physically derived from the surface energy balance. Under the assumption that the available energy (R_n−G) is approximately equivalent between the irrigated and reference surfaces—warranted by their similar topography, albedo, and soil type (Table 2)—the surplus in latent heat flux (ΔLE) is balanced by the deficit in sensible heat flux (ΔH) (Equation 3):

By converting the latent heat flux into a monthly water depth equivalent, the thermo-evaporative conversion coefficient cc is expressed as Equation 4:where ρ_air​ is the air density (1.15 kg m^-3 at 775 m elevation and ∼30 °C), c_p​ is the specific heat of air at constant pressure (1005 J kg^-1·°C⁻¹), Δt_month​ is the duration of an average month (2.63 × 10⁶ s), r_a​ is the aerodynamic resistance (s·m^-1), λ is the latent heat of vaporization (2.43 × 10⁶ J kg^-1 at 30 °C), and ρ_w​ is the density of water (1000 kg m^-3).

For drip-irrigated sparse canopies typical of watermelon fields in arid environments, the aerodynamic resistance rara​ generally ranges between 150 and 400 s m^-1 (Moran et al., 1994; Allen et al., 1998). Propagating this range into Equation 4 yields the theoretically admissible bounds for c (Equation 5):

We acknowledge that r_a​ is not strictly constant and varies dynamically with crop phenology, wind speed, and canopy structure. This variability introduces a structural (model-form) uncertainty in the linear ΔLST-to-ET conversion. Rather than assuming a fixed r_a, we treat c as a probability distribution constrained within these physical bounds, which is subsequently propagated through our Monte Carlo uncertainty analysis (Figure 4).

To evaluate the sustainability of the groundwater system under uncertainty, the probability of aquifer overdraft was framed as a classic risk and reliability problem (Tung and Yen, 2005). The overdraft probability, P(V > R), occurs when the estimated annual pumping volume (V, the “load”) exceeds the annual aquifer recharge (R, the “capacity”). Assuming V and Rare independent random variables with probability density functions fV​(v) and fR​(r), this probability is theoretically defined as:

Because the analytical integration of these complex, mixed distributions (triangular, normal, and uniform) is intractable, we solved this equation numerically using the Monte Carlo method (Morgan and Henrion, 1990). For N = 100,000 iterations, random samples were simultaneously drawn from the uncertainty distributions of the pumping model parameters (to generate vi​) and the recharge distribution (to generate ri​). The probability of overdraft was approximated by the law of large numbers using the indicator function I:

To provide an independent, data-driven calibration of c, we derived monthly reference evapotranspiration (ET₀) for the 2020–2024 period using the FAO-56 Penman-Monteith equation (Allen et al., 1998), computed within GEE from ERA5-Land reanalysis variables (2-m temperature, dew-point temperature, 10-m wind speed, and surface solar radiation). The crop evapotranspiration for drip-irrigated watermelon was then estimated as:where K_c is the stage-specific crop coefficient (K_c,ini = 0.40; K_c,mid = 0.90; K_c,end = 0.75, following Allen et al., 1998), f_w = 0.70 is the fractional soil wetted area correction for drip irrigation (Allen et al., 1998), and f_irr = 0.10 is the mean fractional irrigated cover within the sampled polygons, estimated from high-resolution imagery. For April—the peak growing month—this yields:

Given the mean observed ΔLST for April months over the study period (mean ΔLST_April = 2.54 °C, derived from Table 2), this implies:

It is crucial to note that the fractional irrigated cover (f_irr) is NOT used to compute the total pumped volume, which is directly extrapolated over the entire polygon area to avoid introducing sub-pixel aggregation bias. The f_irr parameter is exclusively utilized in Section 4.2 for the cross-scale validation against field-level agronomic benchmarks.

This calibrated value is in close agreement with the adopted c = 5 mm/month/°C (relative difference of ∼14%), providing first-order validation of the empirical coefficient. The slight overestimate of c__cal relative to c reflects expected spatial heterogeneity in canopy coverage across large polygons and confirms that c = 5 represents a conservative lower bound on actual irrigation-induced evapotranspiration.

This agronomy-informed temporal integration yields the conservative lower-bound estimate (c = 5.0) of V_annual,2022​ = 5.85 Mm³, and the calibrated central estimate (c_cal​ = 5.8) of V_annual,2022​ = 6.79 Mm³. The temporal scaling introduces an additional layer of uncertainty (±10%), which is explicitly accounted for in the subsequent Monte Carlo error propagation (Section 3.5). The complete monthly breakdown used for this annual integration is provided inSupplementary Table S1).

2.3.3 Sustainability assessment

Finally, the estimated annual pumping volume (extrapolated from the monthly values based on the known crop cycle of watermelon) for the year 2022 was compared directly to the total annual aquifer recharge (4.70 Mm³) to assess the sustainability of the practice at the watershed scale. The Methodological workflow is presented below (Figure 5):

FIGURE 5

3 Results

3.1 Hydro-climatic context of the elfeija watershed

To provide a comprehensive overview of the hydro-climatic context of the Elfeija watershed, key water balance components were analyzed for the 2020–2024 period (Figure 6). The precipitation regime is typically Mediterranean-arid, with rainfall concentrated in the autumn and winter months, followed by a pronounced dry season from May to September (Figure 5c). This precipitation pattern directly governs the soil moisture dynamics (Figure 5a), which peaks in late winter before rapidly declining to its minimum during the hot summer months. The cumulative actual evapotranspiration (Figure 5b) illustrates the high atmospheric water demand, particularly during spring and summer. Notably, this analysis visually confirms the exceptional rainfall event in September 2024, which aligns with the anomalous agricultural activity detected in our remote sensing analysis. This strong seasonal disconnect between natural water supply and evaporative demand underscores the critical role of groundwater as the sole reliable water source for the region’s agricultural activities.

FIGURE 6

3.2 Seasonal and inter-annual dynamics of irrigated agro-ecosystem

The synergistic time-series analysis of Sentinel-2 (NDVI) and Landsat (LST) data from 2020 to 2024 reveals a distinct and highly repeatable seasonal signature for the irrigated plots, which stands in stark contrast to the behavior of the natural landscape. Figure 6 presents the monthly aggregated time-series for the mean NDVI and LST for both “Irrigated” and “Reference” polygon groups.

The NDVI profile (Figure 7a) clearly delineates the agricultural calendar. Each year, the irrigated plots exhibit a pronounced vegetation cycle, with NDVI values starting to rise in the winter months (December/January), reaching a peak in early spring (March/April), and rapidly declining to baseline levels by early summer (May/June). The peak NDVI values for irrigated plots consistently reached between 0.20 and 0.28, indicative of a vigorous crop canopy. In contrast, the Reference areas, representing the arid background, show a much flatter and lower NDVI profile, with baseline values consistently remaining below 0.10 throughout the year, showing only a slight response to seasonal precipitation events. This persistent and significant gap between the two curves during the growing season is the primary phenological signature of irrigation.

FIGURE 7

The LST time-series (Figure 7b) provides a complementary thermal perspective that corroborates the NDVI findings. While both irrigated and Reference plots follow the same broad seasonal temperature pattern, a critical divergence occurs during the peak growing season. From approximately February to May each year, the irrigated plots are systematically cooler than the adjacent Reference plots. This period of relative cooling perfectly coincides with the period of elevated NDVI. Following the harvest period in May, the LST of the irrigated plots rapidly converges with that of the Reference plots, and both surfaces exhibit nearly identical thermal behavior throughout the summer, autumn, and early winter (LST difference <1.0 °C), with peak summer temperatures consistently exceeding 60 °C. This recurring pattern of spring cooling followed by summer convergence is the distinct thermal signature of irrigation-induced evapotranspiration.

3.3 Quantifying the impact of irrigation on vegetation and surface temperature

To isolate the net impact of irrigation, we calculated the monthly difference between the irrigated and Reference plots for both NDVI (ΔNDVI) and LST (ΔLST). The results, presented in Figure 8, quantify the magnitude and timing of the agricultural anomaly.

FIGURE 8

The ΔNDVI time-series (Figure 8a) reveals sharp, positive peaks each spring, confirming that the “greenness” surplus is concentrated during the crop growth period. The magnitude of this peak varies from year to year, reaching a maximum of approximately 0.18 in 2020. The ΔLST, calculated as LST_Reference - LST_Irrigated, shows a remarkably similar pattern (Figure 8b). The relative cooling effect is also concentrated in the spring season, with positive peaks indicating that irrigated plots were up to 5 °C cooler than the surrounding natural landscape. A strong and statistically significant temporal correlation was found between the seasonal vegetation surplus (ΔNDVI) and the relative cooling effect (ΔLST) across the study period (Pearson’s r = 0.72, p = 0.019, N = 10). This correlation provides compelling physical evidence that the increased vegetation density is directly linked to an active cooling process, namely, evapotranspiration. Outside the growing season, both ΔNDVI and ΔLST values hover close to zero, confirming the absence of significant agricultural activity.

The physical basis for the observed correlation between ΔNDVI and ΔLST is further illustrated in the LST-NDVI feature space (Figure 9a, b). The scatter plot reveals a clear separation between the two land cover types. The Reference Areas (red dots) are confined to a low-NDVI, high-LST region, forming a well-defined “dry edge” characteristic of an arid landscape. In contrast, the Irrigated Plots (blue dots) expand into a high-NDVI domain that is inaccessible to the natural ecosystem. Significantly, for any given NDVI value, the irrigated points consistently exhibit lower LST values than the reference points, particularly during the irrigation season (blue circles).

FIGURE 9

This vertical displacement within the feature space is a direct visualization of the surface cooling driven by the evapotranspiration of irrigation water, and it forms the conceptual foundation of our ET estimation method.

3.4 Estimation of groundwater pumping volumes and inter-annual variability

The analysis reveals substantial groundwater abstraction during the peak growing season. For the month of April, estimated pumping volumes consistently ranged between 1.1 and 1.8 million m³ across the years. In stark contrast, the estimated volumes for October were minimal, typically less than 0.8 million m³, confirming that the majority of groundwater pumping is concentrated in the spring.

The aggregated monthly results, summarizing the net impact of irrigation and the total estimated pumped volumes for the peak (April) and off-season (October) months, are presented in Table 3 and Figure 10.

FIGURE 10

Significant inter-annual variability is evident. The highest peak-season pumping was estimated for April 2023 at 1.27 Mm³, while April 2024 showed a comparatively lower peak-season volume (0.85 Mm³). Most strikingly, the analysis reveals a major anomaly in late 2024. The estimated pumping volume for October 2024 surged to 1.20 Mm³, a value comparable to a typical peak season and an order of magnitude higher than any other off-season month in the study period. This quantitative result confirms the presence of an anomalous, water-intensive agricultural cycle in late 2024, corroborating the qualitative observations from the NDVI and LST time-series. Based on these seasonal estimates, the total extrapolated annual groundwater abstraction for the studied plots in 2022 is reported as follows: (i) conservative lower-bound estimate (c = 5.0): 5.85 Mm³; (ii) calibrated central estimate (c_cal = 5.8): 6.79 Mm³ [Monte Carlo IC 95%: 3.5–10.8 Mm³]. Both values are used in the sustainability assessment of Section 3.5 and discussed in Section 4.33.5. Water Stress Context and Sustainability Baseline.

To provide a broader context for our pumping estimates, we analyzed two independent, large-scale datasets. First, the long-term climatic water balance for the study area was calculated using the TerraClimate dataset. The results indicate a severe structural water deficit, with a 30-year average annual climatic water balance (Precipitation - Potential Evapotranspiration) of −2284 mm/year (Table 4).

Second, we assessed the baseline water stress using the WRI Aqueduct 3.0 global water risk atlas. The mean Baseline Water Stress score (bws_score) across all our study polygons is 3.76. According to the Aqueduct classification scheme, scores between 3.0 and 4.0 fall into the “High” risk category, bordering on “Extremely High” risk (score >4.0). Figure 11 presents a map of this water stress score over the study area.

FIGURE 11

3.5 Quantitative uncertainty propagation and overdraft risk assessment

While the volume estimates establish the magnitude of the agricultural water demand, providing a single deterministic value (e.g., 6.79 Mm³) obscures the structural and parametric uncertainties inherent in the remote sensing workflow. To rigorously characterize the confidence in these estimates, and to respond to the critical need for uncertainty quantification in data-scarce environments, we performed a Monte Carlo uncertainty propagation analysis (N = 100,000 iterations).The analysis integrated five independent sources of error pertaining to the model inputs (the coefficient c, instrumental ΔLST precision, inter-annual extrapolation, ERA5 ET₀​ precision, and polygon delineation) along with the uncertainty associated with the regional aquifer recharge estimate. The distributions, parameters, and physical justifications for each error source are detailed in Table 5.

The propagation analysis yields a combined relative uncertainty of approximately ±33% around the central estimate. Despite this substantial, though physically realistic, uncertainty bound for a satellite-derived water balance, the probabilistic risk assessment yields a critical finding. When confronted with the aquifer’s sustainable recharge limit (4.70Mm3 ± 15%), the simulation indicates an 82% probability that the actual annual groundwater abstraction exceeds the renewable recharge (P(V > Recharge)). Even at the lower limit of the 68% parametric uncertainty bound (4.55 Mm³), the estimated pumping consumes nearly 97% of the basin’s total annual recharge. This establishes that the diagnosis of severe groundwater overdraft is statistically robust to the dominant sources of parametric uncertainty in our inference model.

3.6 Independent water balance validation results

The three-tier validation framework produced mutually consistent results that corroborate the physical plausibility of our ΔLST-based ET estimates.

For Tier 1, the mean ΔET_MOD16 during the growing season peak (April) across the study period averaged 5.8 mm/month—systematically lower than our ΔLST-derived estimate of 11.9 mm/month (c = 5) but directionally consistent, capturing a positive ET surplus in irrigated plots during each growing season and near-zero values in off-season months. The systematic underestimation of ΔET_MOD16 relative to ΔET_ΔLST is consistent with the documented tendency of MOD16 to underrepresent ET in intensively irrigated arid environments, where the MODIS LAI product fails to fully capture the irrigated canopy density (Velpuri et al., 2013). The seasonal patterns remain in strong agreement (Pearson r = 0.78, p < 0.001, N = 30), validating the temporal structure of our irrigation signal.

For Tier 2, the ERA5-Land water balance residual (ΔGW) was consistently negative during the spring growing season across all five study years (mean ΔGW_April = −10.2 mm/month over irrigated polygons), confirming a net depletion of groundwater storage during the irrigation period. Integrated over the study area (8,780 ha), this corresponds to an estimated net monthly depletion of 0.89 Mm³ during peak season—consistent with our satellite-derived pumping estimates within the bounds of ERA5-Land uncertainty (±10%–15%). During off-season months, ΔGW returned toward neutral or slightly positive values, reflecting limited aquifer recovery during periods of reduced abstraction.

For Tier 3, GLEAM v3.8 ET for the Elfeija region exhibited a clear seasonal peak during the January–May window in all 5 years, with growing-season ET values of 8–13 mm/month—bracketing the lower end of our ΔLST-derived estimates. The exceptional September 2024 anomaly was independently captured by GLEAM as an abrupt ET spike exceeding all other months in the study period, confirming the exceptionality of this event and its attribution to enhanced water availability rather than a gradual trend.

Together (Table 6), these three independent lines of evidence establish a consistent picture: a reproducible and significant enhancement of actual ET over irrigated plots during the growing season, driving a net depletion of groundwater storage that exceeds natural recharge. This multi-source convergence substantially reduces the residual uncertainty attributable to the empirical ΔLST-to-ET conversion, and confirms that our central finding of groundwater overdraft is not an artifact of any single methodological choice.

To avoid any circularity in our validation framework, it is crucial to distinguish between ‘Top-Down’ global model comparisons and ‘Bottom-Up’ local validations. While Tier one and Tier two rely partially on global reanalysis products (ERA5-Land, MOD16), our primary assessment of physical plausibility relies on entirely independent local benchmarks: (1) the agronomic water requirements of watermelon (Table 6), and (2) the local infrastructure density (580 active wells). These bottom-up validations do not depend on ERA5-Land and provide an independent check on the magnitude of our satellite-derived estimates.

4 Discussion

The results of this study provide a multi-faceted quantification of small-scale groundwater pumping in a representative arid agro-ecosystem. By synergistically using Sentinel-2 and Landsat time-series, we have not only detected the clear biophysical signature of irrigated agriculture but also translated it into estimates of water volumes, revealing a critical state of resource unsustainability. This section interprets these findings, discusses their implications in the context of the broader scientific literature, and addresses the limitations of our approach.

4.1 The unambiguous signature of irrigation and the power of multi-sensor synergy

Our primary finding is the clear and repeatable thermal-phenological anomaly observed annually from 2020 to 2024 (Figure 6). The strong positive correlation between the vegetation surplus (ΔNDVI) and the surface cooling effect (ΔLST) (Figures 7, 8A) provides compelling biophysical evidence of irrigation-driven evapotranspiration. This result aligns perfectly with the foundational principles of remote sensing for water stress monitoring, where the LST-NDVI relationship is a cornerstone for assessing surface moisture status and ET fluxes (Sandholt et al., 2002). The convergence of LST values between irrigated and Reference plots immediately after the harvest season (Figure 6B) further strengthens this conclusion, effectively delineating the period of active water application. This capacity to precisely identify irrigated areas and their growing season duration is a significant advancement over methods that rely solely on vegetation indices, which can struggle to differentiate irrigated crops from opportunistic rain-fed vegetation (Conrad et al., 2020). The success of our approach confirms the findings of Massari et al. (2021) and Singh et al. (2020), who highlighted the critical importance of multi-sensor data fusion to overcome the limitations of individual satellites and improve the accuracy of ET estimations in agricultural settings.

4.2 From thermal anomaly to water volumes: inferring the “invisible pumping”

The central contribution of this work is the translation of satellite observations into quantitative estimates of groundwater abstraction. Our estimated peak-season (April) pumping volumes, ranging from 1.1 to 1.8 Mm³ for the studied plots, represent a tangible measure of the previously “invisible” water withdrawals. While the direct quantification of groundwater pumping via remote sensing remains a methodological challenge (Massari et al., 2021; Parra-López et al., 2025), the simplified approach, based on a linear relationship between ΔLST and ET, provides a conservative, first-order estimate of the magnitude of this abstraction. The estimated water consumption is consistent with the high water demand of watermelon cultivation, a practice known to be a major driver of groundwater depletion in other arid regions (Mokhtar and Aram, 2017).

The significant inter-annual variability observed in pumping volumes (Figure 9) reflects the dynamic nature of smallholder farming systems. These fluctuations can be attributed to a complex interplay of factors, including climate variability, farmer’s economic decisions, and energy costs, which are known to influence irrigation strategies in similar contexts (Sidhu et al., 2020; Chaudhuri et al., 2021). Furthermore, the spatial heterogeneity of pumping, as suggested by our per-polygon analysis (Figure 10), underscores the findings of Massuel et al. (2017), who noted significant differences in withdrawal rates among various farm types in North Africa. This highlights the necessity of high-resolution monitoring, as regional or aggregated approaches often fail to capture the localized ‘hotspots’ of over-extraction that are critical for targeted management interventions.

To assess the plausibility of our estimates, we confronted them with local agronomic benchmarks for watermelon cultivation. A direct comparison between scales is, however, only meaningful after accounting for the fractional irrigated cover (f_irr) within the sampled polygons. The satellite-derived ET_irrigation represents an area-averaged flux over large polygons (mean area: 548 ha) that encompass both actively cultivated plots and non-vegetated inter-field areas, whereas the agronomic requirements in Tables 4, 6 refer to the actively irrigated parcel area alone. For April—the peak demand month—our analysis yields a mean polygon-scale ET_irrigation of 11.9 mm/month (119 m³/ha). Scaling this value to the field level requires dividing by the fractional active cropping area within the polygons (f_irr). From high-resolution imagery, the 580 active wells identified within the 8,780 ha irrigated zone, each serving an estimated 1–2 ha of cultivated area (Bossenbroek et al., 2023; Massuel et al., 2017), yield a total active cultivated surface of approximately 580–1,160 ha, corresponding to f_irr ≈0.07–0.13. Applying this correction:

This field-scale estimate falls squarely within the agronomic requirement range of 1,200–1,800 m³/ha for the fruit maturation stage of watermelon (Table 7), providing robust cross-scale validation of our approach Food and Agriculture Organization of the United Nations (2021). Importantly, this reconciliation also clarifies why a conservative c value of 5 mm/month/°C—which operates at the polygon scale—produces aggregate volume estimates that are physically consistent with field-scale crop physiology: the two scales are connected through f_irr, and any systematic underestimation at the polygon scale is a direct reflection of sub-polygon spatial heterogeneity rather than a methodological bias in the thermal anomaly detection itself.

A second layer of validation for our estimation framework comes from in-situ infrastructure data and local hydro-agricultural surveys. Field surveys indicate a high pumping density of approximately one well per 15 ha within the irrigated plots. By allocating our conservative satellite-derived peak water application rate (131 m³/ha in April) across this density, we infer an average peak extraction volume of approximately 1,965 m³/well/month. To cross-check the plausibility of this ‘top-down’ inference, we collated available well pumping estimates from a 2022 survey conducted across six localities within the Elfeija watershed (Table 8; Supplementary Table S2 in Supplementary Material).

The survey of 15 active irrigation wells indicates a mean reported annual pumping volume of 17,601 m³ per well. Given that watermelon cultivation is concentrated over an active growing season of roughly 5–6 months, this annual volume suggests peak in-season monthly pumping rates that are in reasonable agreement with our satellite-derived estimate of 1,965 m³/well/month. This consistency between space-borne thermal anomalies and independent ground-truth data provides strong circumstantial evidence supporting the physical validity of our remote sensing workflow.

This cross-validation also supports our parameterization strategy. The FAO-56/ERA5 calibration procedure yields c_cal​=5.8 mm/month/∘C as the best available central estimate. However, to maintain a cautious approach, we adopt a dual-estimate reporting convention throughout this discussion: a conservative lower-bound estimate based on c = 5.0 (V = 5.85 Mm³) and the calibrated central estimate (c_cal​ = 5.8, V = 6.79 Mm³). The choice of c = 5 for the primary tabulated results reflects a deliberate decision to prioritize the minimum plausible abstraction. By ensuring that our core finding relies on a conservative parameterization, we strengthen the reliability of the overdraft conclusion. The finding that even this lower bound exceeds estimated aquifer recharge by 25% — with the calibrated central estimate exceeding it by 45% — provides a highly probable indication of unsustainable water use.

4.3 A critical state of unsustainability: Pumping vs. recharge

The most critical finding of this study is the direct confrontation of our pumping estimates with local hydrogeological data. Our extrapolated annual pumping of 5.85 Mm³ for 2022 surpasses the total renewable aquifer resource of 4.70 Mm³ by nearly 25% (Figure 5). This strongly suggests a state of groundwater overdraft, where withdrawals exceed the natural replenishment rate. This situation is emblematic of many arid regions worldwide, including the North China Plain (Du et al., 2024), the High Plains Aquifer in the USA (Ajaz et al., 2020), and Iran (Khorrami and Malekmohammadi, 2021), where intensive irrigation has led to severe groundwater depletion. Our results provide quantitative, field-scale evidence for the broader trend of non-renewable groundwater use described by Gleeson et al. (2012), Bierkens and Wada (2019) and Scanlon et al. (2023).

This local, field-scale conclusion is independently corroborated by global-scale water risk assessments. The WRI Aqueduct dataset assigns our study area a Baseline Water Stress score of 3.76, classifying it as a “High” to “Extremely High” risk hotspot. Our study provides the missing link by identifying the primary driver of this high-risk score: the intensive, unmetered groundwater pumping for small-scale agriculture. The “invisible pumping” is, in fact, the principal agent transforming a natural state of aridity into a critical state of human-induced water stress.

This state of unsustainability is further contextualized by the extreme climatic deficit of the region, evidenced by a climatic water balance (P-PET) of −2284 mm/year. Such hyper-arid conditions severely limit natural recharge, making the aquifer highly vulnerable to any level of abstraction. The consequences of such overdraft are well-documented and include declining water tables, increased pumping costs, land subsidence, and soil salinization (Pulido-Bosch et al., 2018; Abd-Elaty et al., 2023), all of which threaten long-term agricultural productivity. Our findings echo the concerns raised by Bossenbroek et al. (2023) regarding the unsustainability of current agricultural policies and practices in Morocco.

4.4 Disentangling climatic signals from anthropogenic forcing: the 2024 anomaly case study

A primary challenge in monitoring arid agro-ecosystems is the ability to distinguish the effects of sporadic, high-intensity rainfall events from those of deliberate irrigation, a difficulty noted in many remote sensing-based water accounting studies (Sachindra et al., 2018; Foster et al., 2020). The anomalous late season of 2024 provided a unique natural experiment to test the capacity of our comparative approach to resolve this ambiguity. Analysis of the ERA5-Land reanalysis data confirmed that this period was characterized by an exceptional rainfall event, which triggered a massive basin-wide evapotranspiration response (Figure 12). The critical question, however, was whether this climatic signal alone could explain the intense vegetation response observed in our time-series.

FIGURE 12

A methodological note is warranted before interpreting the spatially explicit ET estimates of Figure 12. The two quantities compared—ET for Reference Areas and total ET for Irrigated Plots—derive from two physically consistent but formally distinct estimation chains that must be made explicit. For the Reference Areas, the reported mean ET of 9.03 mm/month represents the ERA5-Land total actual evapotranspiration (total_evaporation_sum), providing the regional-scale, rain-fed energy-balance response to the September 2024 precipitation event. This value is extracted at the ERA5-Land native resolution (∼9 km) and applies identically to all polygons within each homogeneous land cover type.

For the Irrigated Plots, the total ET is computed as the additive combination of the same ERA5-Land baseline and the irrigation-induced surplus estimated from the ΔLST approach: ET_total,irr = ET_ERA5,ref + c × ΔLST_irr

For September 2024 specifically: ET_ERA5,ref = 9.03 mm/month (mean over reference polygons, identical spatial driver) and c × ΔLST_irr = 5.0 × 2.57 = 12.84 mm/month, yielding ET_total,irr = 21.87 mm/month—precisely the value reported. The surplus of 12.84 mm/month (∼143%) represents the irrigation-amplified component that cannot be attributed to precipitation alone. This additive decomposition is physically justified under the hypothesis of linear separability of the energy fluxes: if available energy (Rn − G) is approximately equal between the two land cover types—a reasonable assumption given their similar topography, elevation, and soil albedo (Table 1; Table 3) — then the additional latent heat flux in irrigated plots is entirely attributable to supplemental water application, and adds linearly to the rain-fed ET signal captured by ERA5-Land. Three potential residual biases are acknowledged: (1) ERA5-Land at 9 km resolution may slightly overestimate the rain-fed ET contribution over irrigated polygons, introducing a positive bias of ±5%–10% on the surplus estimate; (2) the assumption of equal (Rn − G) introduces an uncertainty of ±3%–7%, manageable given the similar albedo documented by Sentinel-2 imagery; and (3) the temporal mismatch between the monthly ERA5-Land product and the 8-day Landsat revisit cycle is negligible when both are aggregated to monthly composites. These combined biases are incorporated into the uncertainty envelope reported in Section 4.5.

To address this, a spatially explicit analysis of the estimated ET was conducted for all 32 polygons during September 2024 (Figure 13). The results reveal a clear and quantifiable divergence between the two land use types. The Reference Areas, which represent the baseline response of the natural landscape to the precipitation, exhibited a heightened but still moderate mean ET of 9.03 mm/month. In stark contrast, the Irrigated Plots displayed a significantly stronger response, with a mean ET of 21.87 mm/month.

FIGURE 13

This represents a surplus ET of nearly 12.9 mm/month (an increase of ∼143%) in the agricultural plots that cannot be attributed to the rainfall alone. This “anthropogenic amplification” of a natural climatic signal is the direct signature of supplemental water application via groundwater pumping. Figure 7 spatially demonstrates this divergence, showing that the Reference Areas fall almost exclusively into the lower ET classes, while the Irrigated Plots dominate the medium to extreme ET classes. This finding underscores the complex interplay between climate variability and human water use, which has been identified as a key amplifier of groundwater stress in other arid regions (Chen et al., 2023; Elsaidy et al., 2025).

This result is a powerful demonstration of the robustness of our comparative methodology. By systematically contrasting ‘managed’ plots with ‘reference’ areas, our workflow can effectively disentangle climatic drivers from anthropogenic forcing (Figure 14). It proves that while the 2024 rainfall acted as a trigger, it was the active human intervention of pumping that enabled and sustained the intense, water-consumptive agricultural cycle. The increasing severity of droughts in North Africa, as reviewed by Tanarhte et al. (2024), makes such opportunistic water use a critical factor to monitor. This capacity to accurately attribute water use is critical for avoiding the misinterpretation of natural climatic variability as changes in irrigation practices, a critical function for effective water governance.

FIGURE 14

4.5 Methodological limitations and future outlook

A primary limitation of this study is the reliance on indirect lines of evidence and the absence of direct, in-situ validation data. In regions dominated by informal, small-scale farming, records of actual pumping rates or continuous groundwater level time-series (piezometry) are virtually non-existent. Consequently, our satellite-derived volumes must be strictly interpreted as semi-empirical estimates of magnitude and physical plausibility, rather than absolute, exact measurements. To compensate for the lack of in-situ data and avoid circularity, we constrained these estimates using independent bottom-up validations (agronomic benchmarks and well-density infrastructure, Section 4.2), which confirm the physical plausibility of our lower-bound estimates. We acknowledge that extrapolating an instantaneous thermal snapshot (∼10:30 AM) to a monthly water balance introduces temporal scaling uncertainty. While the ‘self-preservation’ hypothesis suggests that the Bowen ratio remains relatively stable during the central hours of the day, diurnal variations in wind speed and cloud cover can decouple the instantaneous LST from the daily cumulative flux. The coefficient c absorbs this structural uncertainty, which is why we treat it as a probability distribution in our Monte Carlo analysis rather than a fixed constant.

To transparently address the remaining uncertainty, we formalized our error propagation in a rigorous Monte Carlo framework (Section 3.5; Table 6). While the resulting combined parametric uncertainty is substantial (±33%), it is a realistic reflection of the challenges inherent in scaling high-resolution satellite thermal anomalies to annual basin-level water budgets without ground-truth calibration. Crucially, this quantitative risk assessment demonstrates that the central conclusion of our study is decision-robust: despite the accumulated uncertainties in both the satellite model parameters (notably the empirical coefficient c) and the regional recharge estimates, the probability of groundwater overdraft remains alarmingly high (P ≈ 82%). Furthermore, our explicit sensitivity analysis (Table 2) demonstrates that the overdraft conclusion holds for any coefficient c > 4.02 mm/month/°C. Since the theoretically admissible range is [3.1; 8.3], the finding of unsustainability is robust across approximately 85% of the physically possible parameter space.

However, it is important to emphasize that our Monte Carlo simulation primarily propagates parametric uncertainty. The linear ΔLST-ET conversion introduces additional structural (model-form) uncertainty that is more difficult to quantify. For instance, extreme surface temperatures (>60 °C) may introduce LST retrieval biases, and dense crop canopies can lead to NDVI saturation effects. Furthermore, our assumption of a homogenous energy balance (Rn​−G) between irrigated and reference plots, while justified as a first-order approximation, may deviate under extreme advective conditions (the ‘oasis effect’) or highly heterogeneous canopy structures. Consequently, from an epistemological standpoint, this distinction is crucial: because our approach relies on a semi-empirical linearization, the resulting distributions do not satisfy the strict frequentist definition of ‘Confidence Intervals’, which implicitly assume a perfect, error-free model structure. Rather, the ranges reported in Table 5 must be interpreted as conditional parametric uncertainty bounds (or model-conditional envelopes). They represent the probability distribution of the pumping volume conditional upon the validity of our first-order structural assumptions. While they do not capture unquantified structural model-form errors, they provide a rigorous, transparent mapping of how parametric variability propagates through the proposed inference chain.

Despite these limitations, the proposed synergistic workflow represents a significant step towards operational monitoring of un-gauged groundwater use. Direct ground-truthing remains the most decisive path toward reducing residual uncertainty in the reported volumes. We recommend prioritizing instrumentation at the highest-stress polygon clusters identified in Figure 11 (bws_score >3.97, e.g., the Anagham and Rgabntel sectors—to be confirmed against polygon-level stress values), notably telemetric piezometers capable of sub-daily water table monitoring, combined with ultrasonic or electromagnetic smart flow meters installed at a representative subset of the 580 identified wells. Such infrastructure would enable a Bayesian recalibration scheme in which the empirical coefficient c is iteratively updated as metered volumes accumulate, progressively narrowing the uncertainty envelope reported in Table 7. The future of sustainable water management in arid regions will likely depend on such hybrid systems, integrating multi-sensor satellite data with ground measurements and hydrological models (Jarlan et al., 2015). Beyond refining the physical energy balance models discussed above, an important direction for future work is to formally assess the portability of the proposed workflow across different arid agro-ecosystems and cropping systems. While the underlying ET_irrigation ≈ c × ΔLST relationship is grounded in generic surface energy balance principles and should hold across contexts, the coefficient c is itself sensitive to canopy architecture and aerodynamic resistance (r_a), and is therefore expected to vary systematically between crop types—e.g., low sprawling cucurbits such as watermelon versus upright cereals or perennial tree crops such as date palms and olives, which dominate other Moroccan oases such as Tafilalet and the Ziz valley. The FAO-56/ERA5 calibration protocol introduced in this study (Section 2.3.2) offers a direct and low-cost pathway for such recalibration, requiring only locally available reanalysis data and crop coefficient tables rather than new field infrastructure. Testing this transferability across a gradient of arid basins with contrasting cropping systems—from other Drâa-Tafilalet oases to broader MENA drylands facing similar unmetered small-scale groundwater abstraction (Hasan and Fatima, 2025; Ahmad and Al-Ghouti, 2020) — would establish whether the workflow can serve as a general-purpose diagnostic tool for regional water agencies. The increasing adoption of machine learning and digital technologies offers promising avenues to improve water demand estimation and crop classification accuracy (Bwambale et al., 2022; Parra-López et al., 2025). Ultimately, our framework offers River Basin Agencies a scalable tool to transition from reactive crisis management to proactive governance—for example, by spatially identifying over-extraction ‘hotspots’ to target field inspections, and by providing data-driven baselines to negotiate sustainable water quotas with local farming communities.

4.6 Economic drivers of groundwater overexploitation

The remote sensing evidence above documents groundwater overdraft but does not explain why individually rational farm-level decisions aggregate into a regionally unsustainable outcome. This is a textbook common-pool resource dilemma (Ostrom, 1990): the private benefit of pumping accrues directly to the farmer, while the cost of aquifer depletion is diffused across all users and deferred in time. Bossenbroek et al. (2023), studying this same watermelon-groundwater system in the Drâa Valley, document that crop profitability, combined with limited well-metering, creates strong incentives for continued private well development. The ongoing shift toward solar-powered pumping under Morocco’s Generation Green strategy—shown to reduce marginal pumping costs by 33%–61% relative to LPG-based pumping (Telleria Juárez et al., 2020) — may further amplify this dynamic, as near-zero marginal costs remove the fuel-cost disincentive that previously constrained extraction.

To ground this narrative empirically, informal farmer consultations conducted during our field visits indicated standard costs for the dominant tenant-farming arrangement: land and water rental (20,000 MAD/ha) plus cultivation setup (30,000 MAD/ha) for a total production cost of ≈50,000 MAD/ha, against average harvest revenue of ≈120,000 MAD/ha—a net margin of ≈70,000 MAD/ha (58%). Normalized against the crop water requirement in Table 6 (3,500–5,150 m³/ha), this implies a private economic value of pumped water on the order of 14–20 MAD/m³ — likely exceeding any feasible volumetric tariff, indicating room for demand-management instruments without threatening profitability. Applied to the f_irr-corrected cultivated area (Section 4.2, 580–1,160 ha), this implies an aggregate annual net profit of ≈41–81 million MAD within the studied zone, illustrating the scale of the private incentive that any intervention must address. The full sourced breakdown is provided in Supplementary Table S3.

These farmer-reported figures are indicative order-of-magnitude estimates from informal consultations rather than a structured survey capturing the full distribution of costs and revenues across farm types (owner-operators versus tenant farmers, local smallholders versus outside investors; Bossenbroek et al., 2023); a dedicated hydro-economic assessment therefore remains a priority for future work, ideally alongside the smart-metering infrastructure discussed in Section 4.5. Beyond validating pumping volumes, such metering would enable the design of incentive-based interventions, tiered water pricing or targeted crop-substitution subsidies, following existing regional precedents, including subsidies for lower-water-demand crops (tomato, onion) and a seasonal water-allocation fee in neighboring Tata province, detailed in Supplementary Table S3.

5 Conclusion

This study presented a synergistic multi-sensor remote sensing workflow to quantify unsustainable groundwater pumping for small-scale agriculture, making the “invisible pumping” in an arid Moroccan watershed both visible and quantifiable. By integrating high-resolution vegetation indices from Sentinel-2 with thermal data from Landsat, a 5-year time-series analysis (2020–2024) successfully identified and monitored the distinct thermal-phenological signature of irrigation.

The main findings are threefold. First, the analysis revealed a clear and repeatable seasonal pattern of agricultural water use, characterized by a distinct spring greening (peak ΔNDVI ≈ +0.10) strongly correlated with a significant surface cooling effect (peak ΔLST ≈2.7 °C). Second, by translating this thermal anomaly into water volumes, a quantitative estimate of groundwater abstraction was produced. The results confirm a critical and unambiguous state of resource unsustainability. The calibrated central estimate of annual pumping for the studied plots in 2022–6.79 Mm³ [IC 95%: 3.5–10.8 Mm³], derived using the FAO-56/ERA5 calibrated coefficient c_cal = 5.8 mm/month/°C—exceeds the total renewable aquifer recharge (4.70 Mm³) by 45%. The conservative lower-bound estimate (5.85 Mm³, c = 5.0) independently confirms overdraft at +25%. The probability that pumping exceeds recharge across the full Monte Carlo uncertainty space is 82%, rising to >95% within the theoretically constrained parameter range. Together, these convergent lines of evidence establish that the observed groundwater overdraft is a robust, parameterization-independent finding, not an artifact of any single modeling choice.

The workflow presented here provides a vigorous, scalable, and cost-effective framework for monitoring un-gauged agricultural water withdrawals in data-scarce regions, offering a tangible tool for data-driven water governance. For policymakers and River Basin Agencies, this framework offers concrete operational applications to transition from reactive to proactive groundwater governance. Specifically, it can be implemented to: (1) spatially identify over-extraction ‘hotspots’ to target field inspections and well-metering prioritization; (2) evaluate compliance with spatial agricultural expansion limits; and (3) provide objective, data-driven baselines to negotiate water quotas and subsidy allocations with local farming communities. By providing the first spatially explicit, multi-annual quantification of groundwater pumping in the Elfeija basin, this study establishes a fundamental quantitative baseline where none previously existed. Future research should prioritize the installation of targeted in-situ piezometers and telemetric smart well-meters at the highest-risk hotspots identified in this study, alongside dedicated hydro-economic surveys of pumping costs and crop profitability. These ground-truth measurements will not only validate the long-term trends detected by this remote sensing framework but will also allow for the dynamic, Bayesian recalibration of the empirical coefficient c, transitioning the proposed workflow into an operational, real-time groundwater management system capable of supporting both regulatory enforcement and economic incentive design.

References

• 1

  Abd-ElatyI.FathyI.KuriqiA.JohnA. P.StrafaceS.RamadanE. M. (2023). Impact of modern irrigation methods on groundwater storage and land subsidence in high-water stress regions. Water Resour. Manag.37 (5), 1827–1840. 10.1007/s11269-023-03457-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  Agence du Bassin Hydraulique de Guelmim-Oued Noun (ABH) (2014). “Etude hydrogéologique de la nappe de la Feija-SMD. Hydroleader, Mission II and III,”. Maroc.

  • Google Scholar
  • View reference in article

• 3

  AhmadA.Al-GhoutiM. (2020). Approaches to achieve sustainable use and management of groundwater resources in Qatar: a review. Groundw. Sustain. Dev.11, 100367. 10.1016/j.gsd.2020.100367

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AjazA.DattaS.StoodleyS. (2020). High Plains aquifer–state of affairs of irrigated agriculture and role of irrigation in the sustainability paradigm. Sustainability12 (9), 3714. 10.3390/su12093714

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AllenR. G.PereiraL. S.RaesD.SmithM. (1998). Crop Evapotranspiration - Guidelines for Computing Crop Water Requirements - FAO Irrigation and Drainage Paper 56.

  • Google Scholar
  • View reference in article

• 6

  AmihaR.KabbachiB.Ait HaddouM.HammaK.BouchritiY.GougueniH. (2024). Climate vulnerability and adaptation strategies in the Zagora Oasis, southern Morocco: a time series analysis and projection to 2050. Theor. Appl. Climatol.155, 9767–9781. 10.1007/s00704-024-05213-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  AmihaR.KabbachiB.Ait HaddouM.BouchritiY. (2025). Spatiotemporal assessment of groundwater quality in the Feija aquifer, southern Morocco: impacts of agricultural practices and climate change. Mediterr. Geosci. Rev.7, 855–873. 10.1007/s42990-025-00172-8

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  BabaW. M.ChehbouniA.OuassanouanY.GascoinS.PaganiniM.OttavianelliG.et al (2025). Monitoring water crisis from space across a Mediterranean region. Sci. Rep.15, 23262. 10.1038/s41598-025-06240-1

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BastiaanssenW. G. M.MenentiM.FeddesR. A.HoltslagA. A. M. (1998). A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation. J. Hydrology212-213, 198–212. 10.1016/S0022-1694(98)00253-4

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BierkensM.WadaY. (2019). Non-renewable groundwater use and groundwater depletion: a review. Environ. Res. Lett.14 (4), 043002. 10.1088/1748-9326/ab1a5f

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  BossenbroekL.FtouhiH.KadiriZ.KuperM. (2023). Watermelons in the desert in Morocco: struggles around a groundwater commons-in-the-making. Water Altern. ISSN 1965-017516 (1), 87–107.

  • Google Scholar
  • View reference in article

• 12

  OuhbiH.BoudharA. (2025). Water sustainability and Policy in the Moroccan context: challenges, solutions, and future perspectives. Int. J. Account. Finance, Auditing, Manag. Econ6 (8), 601–618. Available online at: https://revue.ijafame.com/index.php/home/article/download/1974/1975.

  • Google Scholar
  • View reference in article

• 13

  BwambaleE.AbagaleF. K.AnornuG. (2022). Smart irrigation monitoring and control strategies for improving water use efficiency in precision agriculture: a review. Agric. Water Manag.269, 107324. 10.1016/j.agwat.2021.107324

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  CarlsonT. N.GilliesR. R.PerryE. M. (1994). A method to make use of thermal infrared temperature and NDVI measurements to infer surface soil water content and fractional vegetation cover. Remote Sens. Rev.9 (1-2), 161–173. 10.1080/02757259409532220

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  ChaudhuriS.ParakhD.RoyM.KaurH. (2021). Groundwater-sourced irrigation and agro-power subsidies: Boon or bane for Small/Marginal farmers in India?Groundw. Sustain. Dev.14, 100690. 10.1016/j.gsd.2021.100690

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  ChenH.WuM.DuanZ.ZhaY.WangS.YangL.et al (2023). Forecasting the human and climate impacts on groundwater resources in the irrigated agricultural region of North China Plain. Hydrol. Process.37 (4), e14853. 10.1002/hyp.14853

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  ConradC.UsmanM.Morper-BuschL.Schönbrodt-StittS. (2020). Remote sensing-based assessments of land use, soil and vegetation status, crop production and water use in irrigation systems of the Aral Sea Basin. A review. Water Secur.11, 100078. 10.1016/j.wasec.2020.100078

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  DuJ.LaghariY.WeiY.WuL.HeA.LiuG.et al (2024). Groundwater depletion and degradation in the North China Plain: challenges and mitigation options. Water16 (2), 354. 10.3390/w16020354

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  EamusD.ZolfagharS.Villalobos-VegaR.CleverlyJ.HueteA. (2015). Groundwater-dependent ecosystems: recent insights from satellite and field-based studies. Hydrology Earth Syst. Sci.19 (10), 4229–4256. 10.5194/hess-19-4229-2015

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  El-BouhaliA.AmyayM.El Ouzani Ech-ChahdiK. (2023). Recent variations of water area in the Middle Atlas lakes (Morocco): a response to drought severity and land use changes. Res. Square. 10.21203/rs.3.rs-3224966/v1

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  ElsaidyA.YimerE. A.MogheirY.HuysmansM.VillaniL.Van GriensvenA. (2025). Groundwater drought and anthropogenic amplifiers: a review of assessment and response strategies in arid and semi-arid areas. Sci. Total Environ.978, 179406. 10.1016/j.scitotenv.2025.179406

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  Food and Agriculture Organization of the United Nations (2021). “Watermelon,” in FAO Crop Information. Available online at: https://www.fao.org/land-water/databases-and-software/crop-information/watermelon/.

  • Google Scholar
  • View reference in article

• 23

  FosterT.MienoT.BrozovićN. (2020). Satellite-Based monitoring of irrigation water use: assessing measurement errors and their implications for agricultural water management Policy. Water Resour. Res.56 (10), e2020WR028378. 10.1029/2020WR028378

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  GleesonT.WadaY.BierkensM. F.van BeekL. P. (2012). Water balance of global aquifers revealed by groundwater footprint. Nature488 (7410), 197–200. 10.1038/nature11295

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  HasanF.FatimaB. (2025). A review of drivers contributing to unsustainable groundwater consumption in Pakistan. Groundw. Sustain. Dev.29, 101414. 10.1016/j.gsd.2025.101414

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  JacksonR. D.IdsoS. B.ReginatoR. J.Pinter JrP. J. (1981). Canopy temperature as a crop water stress indicator. Water Resour. Res.17 (4), 1133–1138. 10.1029/WR017i004p01133

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  JarlanL.KhabbaS.Er-rakiS.Le PageM.HanichL.FakirY.et al (2015). Remote Sensing of Water Resources in Semi-Arid Mediterranean Areas: the joint international laboratory TREMA. Int. J. Remote Sens.36 (22), 5679–5717. 10.1080/01431161.2015.1093198

  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  KhorramiM.MalekmohammadiB. (2021). Effects of excessive water extraction on groundwater ecosystem services: vulnerability assessments using biophysical approaches. Sci. Total Environ.799, 149304. 10.1016/j.scitotenv.2021.149304

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  MassariC.ModanesiS.DariJ.GruberA.De LannoyG. J. M.GirottoM.et al (2021). A review of irrigation information retrievals from space and their utility for users. Remote Sens.13 (20), 4112. 10.3390/rs13204112

  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  MassuelS.AmichiF.AmeurF.CalvezR.JenhaouiZ.BouarfaS.et al (2017). Considering groundwater use to improve the assessment of groundwater pumping for irrigation in North Africa. Hydrogeology J.25 (5), 1573–1588. 10.1007/s10040-017-1573-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  MeskourA.AhattabJ.AachibM.HasnaouiM. (2025). Assessing the impact of drought and upstream dam construction on agriculture in arid and semi-arid regions: a case study of the Middle Draa Valley, Morocco. Environ. Monit. Assess.197 (3), 236. 10.1007/s10661-025-13631-4

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  MokhtarA.AramS. (2017). Systemic insights into agricultural groundwater management: case of Firuzabad Plain, Iran. Water Policy19 (5), 867–885. 10.2166/wp.2017.159

  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  MoranM. S.ClarkeT. R.InoueY.VidalA. (1994). Estimating crop water deficit using the relation between vegetation indices and canopy temperature. Remote Sens. Environ.49 (3), 246–263. 10.1016/0034-4257(94)90025-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  MorganM. G.HenrionM. (1990). Uncertainty: A Guide to Dealing with Uncertainty in Quantitative Risk and Policy Analysis. Cambridge University Press. 10.1017/CBO9780511840609

  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  OstromE. (1990). Governing the Commons: The Evolution of Institutions for Collective Action. Cambridge University Press.

  • Google Scholar
  • View reference in article

• 36

  Parra-LópezC.AbdallahS. B.Garcia-GarciaG.HassounA.TrollmanH.JagtapS.et al (2025). Digital technologies for water use and management in agriculture: recent applications and future outlook. Agric. Water Manag.309, 109347. 10.1016/j.agwat.2025.109347

  • CrossRef
  • Google Scholar
  • View reference in article

• 37

  Pulido-BoschA.Rigol-SánchezJ. P.VallejosÁ.AndreuJ. M.CerónJ. C.Molina-SánchezL.et al (2018). Impacts of agricultural irrigation on groundwater salinity. Environ. Earth Sci.77 (10), 386. 10.1007/s12665-018-7386-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  RicciardiV.WaneA.SidhuB. S.GoddeC.SolomonD.McCulloughE.et al (2020). A scoping review of research funding for small-scale farmers in water scarce regions. Nat. Sustain.3 (10), 836–844. 10.1038/s41893-020-00623-0

  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  RosaL. (2022). Adapting agriculture to climate change via sustainable irrigation: biophysical potentials and feedbacks. Environ. Res. Lett.17 (7), 074039. 10.1088/1748-9326/ac7408

  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  SaadiS.SimonneauxV.BouletG.RaimbaultB.MougenotB.FaniseP.et al (2015). Monitoring irrigation consumption using high resolution NDVI image time series: calibration and validation in the Kairouan Plain (Tunisia). Remote Sens.7 (10), 13005–13028. 10.3390/rs71013005

  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  SachindraD. A.AhmedK.RashidM. M.ShahidS.PereraB. J. C. (2018). Statistical downscaling of precipitation using machine learning techniques. Atmos. Res.212, 240–258. 10.1016/j.atmosres.2018.05.022

  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  SandholtI.RasmussenK.AndersenJ. (2002). A simple interpretation of the surface temperature/vegetation index space for assessment of surface moisture status. Remote Sens. Environ.79 (2-3), 213–224. 10.1016/S0034-4257(01)00274-7

  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  ScanlonB. R.FakhreddineS.RatebA.de GraafI.FamigliettiJ.GleesonT.et al (2023). Global water resources and the role of groundwater in a resilient water future. Nat. Rev. Earth and Environ.4 (2), 87–101. 10.1038/s43017-022-00378-6

  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  SidhuB. S.KandlikarM.RamankuttyN. (2020). Power tariffs for groundwater irrigation in India: a comparative analysis of the environmental, equity, and economic tradeoffs. World Dev.128, 104836. 10.1016/j.worlddev.2019.104836

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  SinghR. K.KhandK.KagoneS.SchauerM.SenayG. B.WuZ. (2020). A novel approach for next generation water-use mapping using Landsat and Sentinel-2 satellite data. Hydrological Sci. J.65 (13), 2235–2248. 10.1080/02626667.2020.1817461

  • CrossRef
  • Google Scholar

• 46

  TanarhteM.De VriesA. J.ZittisG.ChfadiT. (2024). Severe droughts in North Africa: a review of drivers, impacts and management. Earth-Science Rev.253, 104701. 10.1016/j.earscirev.2024.104701

  • CrossRef
  • Google Scholar

• 47

  Telleria JuárezR. A.Serghini IdrissiH.SinghabhandhuA.Aw-HassanA.LaamariA. (2020). Economic and environmental benefits of solar-powered groundwater pumping over liquified petroleum gas in Morocco. Asian J. Agric. Food Sci.8 (3), 20–30. 10.24203/ajafs.v8i3.6232

  • CrossRef
  • Google Scholar

• 48

  TuckerC. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sens. Environ.8 (2), 127–150. 10.1016/0034-4257(79)90013-0

  • CrossRef
  • Google Scholar

• 49

  TungY. K.YenB. C. (2005). Reliability and Uncertainty Analyses in Hydraulic Design. New York: American Society of Civil Engineers ASCE, 305.

  • Google Scholar

• 50

  UN-Water (2023). The United Nations World Water Development Report 2023: Partnerships and Cooperation for Water. Paris: UNESCO.

  • Google Scholar

• 51

  VelpuriN. M.SenayG. B.SinghR. K.BohmsS.VerdinJ. P. (2013). A comprehensive evaluation of two MODIS evapotranspiration products over the conterminous United States: using point and gridded FLUXNET and water balance ET. Remote Sens. Environ.139, 35–49. 10.1016/j.rse.2013.07.013

  • CrossRef
  • Google Scholar