A deep learning approach for assessing environmental stress dynamics in the Makkah watersheds, Saudi Arabia

The integration of trend-based Environmental Stress Severity (ESS) and Artificial Neural Network (ANN) modeling provides a robust framework for quantifying natural environment dynamics and forecasting stress in arid ecosystems. This study focuses on the Makkah Watersheds, Saudi Arabia. ESS was derived from Kendall’s τ and significance thresholds applied to long-term Landsat indices (NDVI, NDWI, LST; 1991–2023) and hydroclimatic drivers from TerraClimate (temperature, precipitation, AET, CWD). Trend analysis revealed significant declines in NDVI (20.48 km²) and NDWI (99.57 km²) and an increase in LST (156.5 km²) at p < 0.05. Hydroclimatic pressures intensified as CWD expanded (1,936 km²) while AET (176 km²) and precipitation (128 km²) contracted. The optimized ANN models demonstrated strong predictive performance for Environmental Stress Severity, with the baseline ANN achieving R² = 0.881 (MAE = 0.200), and further improvements obtained using metaheuristic optimization (ANN-PSO: R² = 0.890, MAE = 0.186; ANN-GWO: R² = 0.883, MAE = 0.196). Zonal and landcover diagnostics showed the highest stress in sparsely vegetated zones (e.g., WS5: 3.49 ± 0.67, WS4: 3.31 ± 0.71), whereas built-up areas exhibited lower means (e.g., WS1: 2.88 ± 0.72). Right-skewed ESS distributions in natural covers highlight asymmetric dryland vulnerability, characterized by a larger tail of highly stressed pixels. The workflow is scalable, temporally sensitive, and directly actionable for restoration planning, land-use regulation, and climate-risk adaptation in drylands.

1 Introduction

Environmental degradation in arid and semi-arid regions threatens ecological sustainability, water security, and climate resilience. Covering >40% of Earth’s land surface, these regions are especially exposed to climate variability and change (UNCCD, 2022). Robust, spatially resolved understanding of vegetation dynamics under changing hydroclimatic conditions is therefore essential for effective restoration and adaptive management in drylands.

Saudi Arabia occupying ∼80% of the Arabian Peninsula is predominantly arid, with high temperatures and low, erratic rainfall (Al-Rowaily et al., 2016). In the Makkah Province, projected increases in temperature and evapotranspiration indicate escalating environmental stress under SSP2-4.5 and SSP5-8.5 scenarios (El-Rawy et al., 2023), while rapid urbanization elevates flood risk and impact (Al-Areeq et al., 2021). Satellite remote sensing enables consistent observation of these dynamics across scales; the near-infrared (NIR, ∼0.76–0.90 µm) and thermal infrared (TIR, ∼10.4–12.5 µm) bands are central for monitoring vegetation condition, surface water, and land-surface temperature (Lambin, 2001). The Landsat program’s multi-decadal record provides an empirical basis for trend detection (Wulder et al., 2022).

Recent advances in cloud platforms and open libraries as Google Earth Engine (Gorelick et al., 2017) and Google Colab (Bisong, 2019), TensorFlow and scikit-learn (Abadi et al., 2016; Pedregosa et al., 2011), and geospatial tooling in R/Python (R Core Team, 2024; Hijmans, 2024a; 2024b; Kuhn, 2008) enable scalable, reproducible analyses of large Earth-observation archives. These tools support machine-learning and deep-learning models capable of capturing non-linear responses and emerging degradation hotspots (Schmidhuber, 2015; Hasan et al., 2024; Mas et al., 2004). Recent studies have demonstrated the growing utility of deep learning in vegetation stress and desertification monitoring. Sun et al. (2023) applied a bidirectional long short-term memory (BiLSTM) model to predict NDVI-based vegetation stress across northern China using meteorological data, while Paul et al. (2025) provided a comprehensive review highlighting how deep learning architectures can effectively detect plant stress from multi-sensor imagery. Additionally, threshold-based vegetation index approaches have been successfully used for desertification assessment, reinforcing the rationale for trend- and index-driven stress diagnostics applied here Berdyyev et al. (2024).

This study develops a spatially explicit Artificial Neural Network (ANN) framework to predict an Environmental Stress Severity (ESS) index for watersheds in the Makkah Province. We integrate long-term Landsat-derived indices—NDVI, NDWI, LST (Tucker, 1979; Gao, 1996; Jiménez-Muñoz and Sobrino, 2003) with TerraClimate variables (temperature, precipitation, AET, CWD; Abatzoglou et al., 2018) and ESA land cover (ESA, 2022). Trend features are derived using Kendall’s τ and significance thresholds and supplied to the ANN to generate continuous ESS predictions. In the regional literature, most applications emphasize trend statistics or GAMs; while deep learning has been explored for desertification and urban ecosystem mapping (Said et al., 2021; Alqadhi et al., 2024), a multidecadal, spatially explicit ANN for environmental stress prediction across the Makkah watersheds remains under-developed. We complement the predictive layer with a multi-criteria evaluation (MCE) to identify restoration suitability.

The study addresses the following research objectives: (i) what are the long-term trends in NDVI, NDWI, and LST (1991–2023) in the study area; (ii) whether a shallow Artificial Neural Network, using normalized NDVI, NDWI, LST, temperature, precipitation, AET, and CWD, can accurately predict an Environmental Stress Severity (ESS) index; and (iii) how the predicted ESS varies spatially across watersheds and land-cover classes to reveal stress hotspots relevant to restoration.

Accordingly, we (1) quantify multi-decadal trends in NDVI, NDWI, LST and climatic conditions; (2) develop and validate an ANN to predict ESS from normalized spectral and climatic predictors; and (3) generate spatial ESS maps and conduct zonal analyses by watershed and land cover to identify stress hotspots for planning.

2 Materials and methods

2.1 Study area

The study was conducted across eight watersheds in the Makkah Province of western Saudi Arabia. These watersheds lie between 21.1458°N and 21.6047°N latitude, and 39.5708°E to 40.2881°E longitude, with elevations ranging from 124 m to 2,479 m above sea level (Figure 1). The water flow generally follows an east-to-west direction, shaped by the escarpment topography of the western Arabian Shield. The region spans a variety of landforms such as mountains, hills, plains, and tablelands, offering a diverse landscape for analyzing environmental stress across terrain types. The dominant climate zone is Tropical Desert, with additional presence of Subtropical Dry, Tropical Dry, and Warm Temperate Dry types (Karger et al., 2017). These climatic gradients influence vegetation dynamics, soil moisture, and surface energy fluxes, contributing to heterogeneous environmental stress patterns.

Figure 1

Figure 1. (a) Location of Makkah Province in Saudi Arabia (b) Watersheds in Makkah Province (c) The area of interest, highlighting Landuse/landcover (LULC) within the eight watersheds.

Land use and land cover (ESA, 2022) analysis highlights a predominance of bare areas (1,303.94 km²) and sparse vegetation (973.20 km²). Urban or artificial surfaces cover 185.37 km², while cropland (67.95 km²) and grassland, scrub, or shrubland (31.06 km²) represent smaller portions of the landscape. Only 0.18 km² is covered by deciduous forest. This LULC composition reflects the arid to semi-arid environment, where vegetation is limited and unevenly distributed.

The physical and climatic complexity of the study area makes it a suitable testbed for deep learning models aimed at identifying degradation hotspots and evaluating landscape vulnerability across environmental gradients.

2.2 Schematic workflow

The schematic workflow (Figure 2) illustrates the integrated framework applied to assess environmental stress and identify priority sites for vegetation restoration. Multi-source datasets, including Landsat imagery, TerraClimate variables, SoilGrids, and ESA CCI land cover, were systematically processed using Google Earth Engine to ensure consistency. Temporal trends were derived with Kendall’s tau and Mann–Kendall tests. These indicators, combined with climate predictors, were modeled through an Artificial Neural Network (ANN) to generate the Environmental Stress Severity (ESS) index.

Figure 2

Figure 2. Workflow showing integration of geospatial datasets and processing steps. Landsat, climate and land cover were pre-processed, analyzed with Kendall’s trend, and modeled using an ANN to generate the Environmental Stress Severity (ESS) index.

2.3 Data sources and processing

To assess environmental stress using deep learning techniques, multi-source geospatial datasets were employed and systematically processed, as shown in Table 1. The data sources include optical and thermal satellite imagery from the Landsat archive, climate variables from the TerraClimate dataset, and land cover information from ESA’s Climate Change Initiative (CCI). Data processing was performed using Google Earth Engine (GEE) to ensure scalability, reproducibility, and temporal consistency across datasets.

Table 1

Table 1. Data sources, variables, native resolution, and references used in this study.

2.3.1 Landsat data processing

Annual time series from Landsat 5 TM, 7 ETM+, and 8 OLI/TIRS (1991–2023) were obtained in Google Earth Engine for the study area (WRS-2 path 166, rows 44–45) (Chander et al., 2009; Wulder et al., 2012; Roy et al., 2014; Gorelick et al., 2017). Processing used Collection-2 Level-2 surface reflectance; pixels flagged as cloud, cloud shadow, snow, or saturation were removed with QA_PIXEL/BQA bitmasks. For each year, median surface-reflectance composites were generated and NDVI/NDWI computed. Land surface temperature (LST) was retrieved from the thermal bands, Band 6 (TM/ETM+) and Band 10 (TIRS) via a radiative-transfer approach with emissivity correction, converted to °C, and resampled to a 30-m analysis grid (Jiménez-Muñoz and Sobrino, 2003). This workflow limits phenological and sensor inconsistencies while preserving long-term trends.

2.3.2 Climate data processing

Monthly TerraClimate fields (tmax, precipitation, actual evapotranspiration, climatic water deficit) at ∼4 km resolution were retrieved for 1991–2023 (Abatzoglou et al., 2018). For each pixel and year, tmax was averaged, while precipitation, AET, and CWD were summed to produce annual series aligned with the Landsat period. Pixel-wise monotonic trends were then quantified using Kendall’s tau-b (Equation 6), with significance assessed at α = 0.05. All rasters were reprojected, resampled (bilinear), and snapped to the Landsat analysis grid for co-registration; resampling does not increase spatial detail beyond TerraClimate’s native resolution.

2.3.3 Landuse/landcover (LULC)

Annual ESA CCI Land Cover was retrieved and the latest year within the study window was selected (ESA, 2022; Harper et al., 2023). The native legend was collapsed to five study-relevant classes, bare ground, shrubland, cropland, urban/built-up, and sparse vegetation, using a documented crosswalk. The map was reprojected to the analysis CRS and snapped to the Landsat grid via majority resampling, then clipped to the AOI. The resulting LULC layer was used to support interpretation of modeled stress patterns.

2.3.4 Hydrological analysis

A 30 m global DEM was processed in Google Earth Engine. The DEM was hydrologically conditioned (depression filling), then used to derive elevation, slope, plan/profile curvature, flow direction, and flow accumulation. A channel network was extracted from flow-accumulation using an empirically tuned threshold and pruned to remove spurious headwater segments. Pour points were placed at the main drainage outlets within the AOI, and eight sub-basins (WS1–WS8) were delineated by watershed tracing. All rasters were reprojected and snapped to the Landsat analysis grid to ensure co-registration.

2.4 Trend analysis

Annual median spectral indices (1991–2023) were evaluated using Kendall’s tau-b, a nonparametric, tie-corrected measure appropriate for monotonic trends in non-normal environmental series (Kendall, 1975). Tau-b was computed using Equation 1:

${\tau }_b=\frac{C-D}{\sqrt{\left(C+D+T_x\right)\left(C+D+T_y\right)}}(1)$

where C and D are the counts of concordant and discordant pairs, and Tx and Ty denote ties in the time (x) and index (y) variables, respectively. Statistical significance of trends was assessed with the two-sided Mann–Kendall test (Mann, 1945; α = 0.05). To compare mapped extents of positive versus negative trends across predefined groups (watersheds or land-cover classes), we applied a paired t-test to the corresponding area measures, thereby testing for systematic differences in trend areas.

2.5 Model development, evaluation, and validation

An Artificial Neural Network (ANN) model was developed to assess spatial patterns of environmental stress by predicting the Environmental Stress Severity (ESS) index based on normalized environmental predictors (Figure 3).

Figure 3

Figure 3. Schematic of the feed-forward Artificial Neural Network (ANN) architecture used for ecological stress prediction. Input variables were normalized (AET, NDVI, LST, NDWI, temperature, rainfall, and climatic water deficit), processed through a single hidden layer of 10 neurons, and mapped to a continuous ecological stress output.

2.5.1 Environmental stress severity (ESS)

For each indicator (NDVI, NDWI, LST), a severity score (S_v) ranging from one to five was assigned according to the criteria in Table 2. The final Environmental Stress Severity (ESS) was defined as the dominant stress signal, expressed as Equation 2

$\text{ESS}=\max \left\{S_{\text{NDVI}},S_{\text{NDWI}},S_{\text{LST}}\right\}(2)$

Table 2

Table 2. Decision rules for per-indicator severity $S_v$.

The selected cut-offs also follow common practice in time-series and hydroclimatic trend applications, where τ magnitudes around 0.10, 0.30, and above 0.50 are treated as weak, moderate, and strong monotonic associations, respectively, and are used as interpretable bins for mapping spatial gradients (Yue et al., 2002; Wang et al., 2020; de Beurs and Henebry, 2005). A formal sensitivity analysis was not conducted because the classification was anchored to (i) these widely used τ-strength conventions and (ii) the observed τ distribution in the study area; additionally, the stability of the severity mapping was checked through independent cross-sensor agreement with MODIS NDVI trend behavior.

2.5.2 Input variables

Predictors comprised normalized trends of NDVI, NDWI, LST and annual climate drivers (tmax, precipitation, CWD, AET). Each predictor was scaled to [0,1] via min–max normalization to stabilize training using Equation 3:

$x_{normalised}=\frac{x-\mathit{min}\left(x\right)}{\mathit{max}\left(x\right)-\mathit{min}\left(x\right)}(3)$

2.5.3 ANN model and training

We used a feed-forward Artificial Neural Network (ANN) with a single hidden layer, an architecture widely applied for ecological prediction tasks where nonlinear but relatively smooth input–response relationships are expected. To avoid overfitting, we restricted model complexity and selected the final architecture (one hidden layer with 10 neurons) after exploratory hyperparameter tuning in which we evaluated combinations of one to two hidden layers and 5–20 neurons using fivefold internal validation. The selected configuration provided stable performance gains without unnecessary model depth. The hidden layer used a logistic (sigmoid) activation, which performed comparably to ReLU and Tanh during tuning and offered greater numerical stability for our normalized [0–1] predictors, while the output layer used a linear activation appropriate for continuous regression. Training was performed using backpropagation with weight decay (0.01) for regularization and a maximum of 200 iterations. Data were randomly split 70/30 into training and testing sets; only training data informed weight updates.

To further assess whether metaheuristic optimization improves predictive performance beyond a classical shallow ANN, we additionally implemented Particle Swarm Optimization (ANN-PSO) and Grey Wolf Optimization (ANN-GWO) to tune key ANN hyperparameters (number of hidden neurons and weight decay). Both algorithms searched predefined parameter ranges and optimized model performance using internal validation. The optimized ANN-PSO and ANN-GWO models were then evaluated on the same held-out test set and compared directly with the baseline ANN using identical predictors and evaluation metrics (MAE, RMSE, R²), enabling an explicit assessment of performance gains attributable to post-2020 ANN optimization strategies.

2.5.4 Model evaluation

Performance was quantified on the held-out test set using mean squared error (MSE), mean absolute error (MAE), and coefficient of determination (R²) (Chai and Draxler, 2014; Willmott and Matsuura, 2005; Kvalseth, 1985). These complementary metrics assess average error magnitude, variance explained, and sensitivity to outliers.

To interpret the ANN outputs at scale, we trained a compact gradient-boosting surrogate on the same predictors and computed SHAP values on a large subset.

2.5.5 Feature importance assessment

Global feature importance was assessed using a permutation-based, model-agnostic approach implemented in the iml framework. Each predictor was randomly permuted while holding others constant, and the resulting increase in mean absolute error (MAE) was recorded as its importance score. Larger values indicate greater contribution to predictive performance.

2.5.6 External cross-sensor validation with MODIS Terra/Aqua NDVI

To test whether ANN-based stress (ESS) captures vegetation dynamics independently of the training sensor, we performed an external, cross-sensor validation using MODIS Terra/Aqua NDVI (500 m) from 2000 to 2021 at five sites (points a–e) (Figure 6). Daily NDVI was aggregated to monthly medians, and we also formed deseasonalized anomalies (monthly value minus the site-specific monthly climatology). For each site we estimated trend magnitude with Sen’s slope (Theil–Sen; NDVI units · yr^-1) and monotonicity with Kendall’s τ (two-sided p). Site stress was fixed to the ANN based stress values. Agreement with the a priori hypothesis (higher stress ⇒ more negative NDVI trend) was evaluated two ways: (i) within-site trend signs and (ii) across-sites rank associations (Spearman’s ρ, Kendall’s τ) between stress and Sen’s slope. We also summarize decadal slopes (slope × 10) and total change over 2000–2021 (slope × 21).

3 Results

3.1 Spatiotemporal trends in environmental indicators (1991–2023)

Pixel-wise Kendall’s τ revealed heterogeneous but coherent trends across the eight watersheds (Table 3). Areas meeting the significance threshold (α = 0.05) and, separately, marginal significance (0.05–0.10) were extracted to quantify change.

Table 3

Table 3. Areal extent (km²) of significant monotonic trends (Mann–Kendall) by variable and sign for 1991–2023. Values are reported for p < 0.05 and for marginal significance (0.05–0.10).

Greenness and moisture increased over substantially larger areas than they declined: NDVI showed 282.46 km² of positive trends versus 20.48 km2 negative (p < 0.05), and NDWI 248.44 km2 positive versus 99.57 km2 negative. In contrast, surface warming and atmospheric demand intensified: LST increased over 156.45 km2 (versus 2.30 km2 decreasing; p < 0.05) and maximum temperature was positive over 2,016 km². Hydroclimate indicators pointed to drying pressure, AET and rainfall declined over 176 km² and 128 km², respectively, while climatic water deficit increased over 1,936 km² (all p < 0.05). Together, these results indicate a net rise in thermal and moisture stress, co-occurring with localized greening (Figure 4).

Figure 4

Figure 4. Spatial distribution of Kendall’s τ (1991–2023) for major environmental variables: (A) NDVI, (B) NDWI, (C) LST, (D) maximum temperature (Tmax), (E) rainfall, (F) actual evapotranspiration (AET), (G) climatic water deficit (CWD), and (H) elevation. Color bars represent Kendall’s τ values for panels A–G and elevation in meters (m (A) s.l.) for panel (H).

3.2 Model-based prediction of environmental stress severity

3.2.1 Model performance

The baseline simple ANN was compared with two metaheuristic-optimized variants, ANN–PSO and ANN–GWO. The baseline ANN achieved a mean absolute error (MAE) of 0.200 and an R² of 0.881, indicating strong baseline predictive skill. Optimization via Particle Swarm Optimization further improved performance (ANN–PSO: MAE = 0.186, R² = 0.890), yielding the lowest prediction error and highest explanatory power. The ANN–GWO model also outperformed the baseline ANN (MAE = 0.196, R² = 0.883), though with slightly lower gains than ANN–PSO (Table 4). The results demonstrate that metaheuristic optimization enhances ANN performance, with ANN–PSO providing the most robust balance between accuracy and generalization for environmental stress prediction. Model choices follow established practice for ecological prediction (Lek and Guégan, 1999; Ripley, 1996; Kuhn and Johnson, 2013).

Table 4

Table 4. Performance comparison of the baseline ANN and metaheuristic-optimized ANN models (PSO and GWO) for Environmental Stress Severity prediction.

3.2.2 Distribution of Predicted Environmental Stress Severity

The distribution of baseline ANN-predicted Environmental Stress Severity (ESS) values exhibits a pronounced right-skewed structure (Figure 5), with most pixels clustered within the low-to-moderate stress range (approximately ESS = 1.5–2.5) and a steadily declining frequency toward higher values. This asymmetric pattern reflects the heterogeneous spatial expression of environmental stress, whereby extensive portions of the landscape experience relatively moderate conditions, while smaller, spatially concentrated areas form distinct high-stress hotspots. The extended right tail indicates localized zones of intensified stress, consistent with persistent vegetation decline, moisture deficit, and elevated surface temperatures.

Figure 5

Figure 5. Histogram of baseline ANN-predicted Environmental Stress Severity (ESS) values across the study area.

These distributional characteristics of ESS across the study area given in Figure 6. Predicted stress forms a coherent belt across the central and eastern watersheds, particularly along steeper, high-elevation ridges, whereas flatter lowland areas exhibit substantially lower stress levels. Watershed-level statistics reinforce this gradient: WS4 and WS5 show the highest mean ESS (2.62 and 2.40, respectively), with WS4 also displaying the greatest variability and pronounced hotspots. WS1 and WS7 exhibit moderately elevated stress, especially across mid-to high-elevation terrain, while WS6 and WS8 consistently register the lowest stress levels; WS2 and WS3 fall within an intermediate range. Collectively, these patterns indicate that ESS intensifies with increasing elevation and terrain complexity, highlighting topographically driven vulnerability within the watershed system.

Figure 6

Figure 6. Predicted Environmental Stress Severity (continuous ESS) from the baseline simple ANN across eight watersheds (1–8) and points a–e for cross-sensor validation.

3.2.3 Feature importance of predictors

Permutation-based importance shows that LST trend is the dominant predictor (importance = 2.87), followed by NDWI trend (2.06) and NDVI trend (1.40), indicating that thermal and moisture stresses exert stronger control on environmental stress severity than rainfall or atmospheric water demand (Figure 7).

Figure 7

Figure 7. Global permutation-based feature importance of predictors used in the ANN model. The bars show the increase in mean absolute error (MAE) after permuting each predictor, computed using a model-agnostic approach (imlFeatureImp).

3.2.4 Cross-sensor validation with MODIS Terra/Aqua

Trends from monthly NDVI (anomalies) matched the hypothesized pattern at all sites ‘a’ to ‘e’ as shown in Figure 8.

• High-stress: a (stress = 4.138): τ = −0.430 (p < 0.001), Sen = −0.00289 years^-1 (≈−0.0289 per decade; ≈ −0.0607 total). c (3.603): τ = −0.251 (p = 1.7 × 10⁻⁹), Sen = −0.00102 years^-1 (≈−0.0102 per decade; ≈ −0.0214 total).

• Low-stress: b (0.673): τ = 0.672 (p < 0.001), Sen = +0.00976 years^-1 (≈+0.0976 per decade; ≈ +0.2050 total). d (1.041): τ = 0.379 (p < 0.001), Sen = +0.00160 years^-1 (≈+0.0160 per decade; ≈ +0.0336 total).

• Mid-stress: e (2.468): τ = 0.320 (p < 0.001), Sen = +0.00053 years^-1 (≈+0.0053 per decade; ≈ +0.0111 total).

Figure 8

Figure 8. External cross-sensor validation of ANN-derived environmental stress using MODIS NDVI (Terra/Aqua, 500 m; 2000–2021).

Across sites, stress ranked perfectly opposite to NDVI trend in both the raw monthly series and the anomalies (identical outcomes): Spearman ρ = −1.000 (p = 1.12 × 10⁻²³); Kendall τ = −1.000 (p = 0.0275). The mean NDVI trend in high-stress sites was −0.00196 years^-1, versus +0.00568 years^-1 in low-stress sites (difference = −0.00764 years^-1), indicating materially divergent trajectories consistent with the ESS signal.

These results demonstrate strong cross-sensor consistency: locations the ANN classified as highly stressed exhibit significant NDVI decline in an independent satellite record, while low-stress locations greened over the same period. Using deseasonalized anomalies confirms that the relationships are not artifacts of seasonality.

3.2.5 Zonal assessment by watershed and land cover

Baseline simple ANN–derived ESS values varied systematically by land-cover type and watershed (Figure 9). Sparse vegetation exhibited the highest means (e.g., WS5: 3.49 ± 0.67; WS4: 3.31 ± 0.71), indicating elevated vulnerability of open rangelands. Grassland/scrub also showed elevated stress in WS1 (3.04 ± 0.46). Cropland was heterogeneous, from high stress in WS3 (3.85 ± 0.77) to low values in WS7 (1.47 ± 1.05), reflecting management contrasts. Artificial surfaces had lower means overall (e.g., WS1: 2.88 ± 0.72; WS6: 2.15 ± 1.06) but included local extremes (maximum ESS = 6.92 in WS5), consistent with urban heat-island hotspots. Bare areas were comparatively stable (e.g., WS5: 3.04 ± 0.80; WS2: 2.75 ± 0.67).

Figure 9

Figure 9. Mean (±) Baseline simple ANN–derived ESS by land-cover class within each watershed.

To characterize distributional differences beyond means, kernel densities of ESS were estimated for land-cover classes within each watershed. Densities highlight class-specific stress regimes, e.g., higher stress concentration in sparsely vegetated zones (Figure 10).

Figure 10

Figure 10. Density distributions of Baseline simple ANN–derived ESS by land-cover class for WS1–WS8; peaks indicate dominant stress levels within each class.

4 Discussion

This study integrates multi-decadal satellite and climate records with an ANN to produce spatially explicit estimates of environmental stress. The approach responds to calls for scalable, data-driven diagnostics that can guide land restoration and policy targeting (UNCCD, 2022; Poldveer et al., 2019).

4.1 Data integration and modelling

The framework couples long-term spectral and hydro-climatic trends with a shallow ANN, capturing non-linear responses without strong distributional assumptions (LeCun et al., 2015). Trend-based features improve temporal sensitivity and reduce noise from interannual variability, enabling the model to detect legacy effects of stress (Ma et al., 2019; Bürgi et al., 2017). Normalization and regularization supported stable training and generalization, consistent with recent dynamic environmental modelling (Banerjee et al., 2025). Previous studies have reported ANN superiority over SVM and RF models in similar ecological forecasting contexts (Etensa et al., 2025; Kladny et al., 2024). Three ANN variants were evaluated: a baseline simple ANN and two metaheuristically optimized models, ANN-PSO and ANN-GWO. Performance improved with optimization, as ANN-PSO achieved the lowest error (MAE = 0.186) and highest explanatory power (R² = 0.890), followed by ANN-GWO (MAE = 0.196, R² = 0.883), compared to the simple ANN (MAE = 0.200, R² = 0.881), highlighting the benefit of post-2020 optimization strategies. The ANN therefore provided the best balance between accuracy, generalization, and smoothness of the predicted surface, reflecting its ability to capture nonlinear interactions among NDWI, NDVI, and LST interactions that are especially important in transitional ecotones where vegetation–climate responses diverge from linear assumptions. The harmonized, co-registered inputs further address spatial comparability challenges common in degradation mapping (Yang et al., 2025). Although the modelling components used here are well established, the integrated use of long-term trend features (Kendall’s τ) with ANN-based prediction provides a novel application for continuous environmental stress mapping in the Makkah watersheds. This trend-driven feature engineering and watershed-level zonal assessment extend existing approaches and offer a replicable framework for dryland restoration planning.

4.2 Patterns and environmental implications

Predicted stress exhibits strong spatial structure aligned with key biophysical gradients: elevated values dominate open, transitional ecotones. This configuration reflects established dryland vulnerability theory, where sparse canopies amplify thermal loading and water imbalance (D’Odorico et al., 2013; Huang et al., 2016). The SHAP analysis reinforces this interpretation by identifying NDWI trend as the dominant predictor, exceeding both NDVI and LST in marginal contribution to stress. In semi-arid systems, vegetation water content responds more quickly and sensitively to moisture deficits than greenness, and NDWI captures liquid-water dynamics in foliage and upper canopy layers, often signaling hydraulic stress before visible declines emerge in NDVI. The sharp separation between high- and low-stress watersheds thus suggests that sustained moisture loss, rather than greenness alone, is the primary driver of ecosystem vulnerability in this landscape.

Class-wise diagnostics show that shrublands and bare areas disproportionately absorb climatic variability, echoing patterns documented in North African and Central Asian arid belts (El Najjar et al., 2023; Ensor et al., 2019). Cropland responses are heterogeneous consistent with management-dependent buffering noted in precision agriculture (Kamilaris and Prenafeta-Boldu, 2018), while some built-up zones exhibit localized greening, likely tied to urban vegetation maintenance. Overall, the spatial and predictor-based evidence converges on moisture limitation as the central mechanism shaping stress distribution across the region.

4.3 Limitations and areas for advancement

Empirical reliability depends on input representativeness, and several constraints should be acknowledged. Medium-resolution LULC and 30 m Landsat inputs may overlook fine-scale stress processes within heterogeneous mosaics of shrubs, bare soil, croplands, and urban–rangeland edges. Likewise, the ESS index is derived from multi-decadal monotonic trends, providing robust long-term signals but potentially missing short-lived or acute disturbances such as anomalously hot years, pest outbreaks, or sudden drought pulses. The ANN is inherently correlative, identifying statistical associations rather than mechanistic drivers of stress. Hybrid models that couple data-driven learning with process constraints could therefore improve realism (Olawade et al., 2024), while temporal ensembles and attention mechanisms offer promise for capturing extremes and emerging regime shifts (Vaswani et al., 2017). Integrating higher-resolution imagery, seasonal indicators, and mechanistic components represents a logical next step to enhance sensitivity and predictive robustness. A key limitation of this study is that model evaluation was based on a random 70/30 train–test split rather than spatial cross-validation. Because the predictor variables and target response exhibit spatial autocorrelation, randomly partitioned samples may remain spatially proximate, potentially leading to mildly optimistic estimates of predictive accuracy due to spatial leakage.

4.4 Management implications and future outlook

By mapping stress and suitability jointly, the workflow supports pre-emptive restoration and adaptive land management. Concentrations of high stress near urban–marginal interfaces highlight priorities consistent with regional sustainability targets. The approach is implementable within environmental early-warning and planning systems and is transferable to other drylands confronting similar pressures (IPCC, 2023). As aridity intensifies under climate change, combining trend-aware learning with MCE offers a pragmatic, policy-ready pathway to allocate restoration where feasibility and impact intersect.

5 Conclusion

This study demonstrates that Artificial Neural Networks (ANNs) provide a robust, systems-based means to model environmental stress in arid landscapes. By integrating multi-decadal remote-sensing indices (NDVI, NDWI, LST) with climate variables (temperature, precipitation, AET, CWD), the framework captured non-linear controls on vegetation stress across the Makkah watersheds, achieving high predictive skill (R2 = 0.88; MAE = 0.20). The resulting Environmental Stress Severity (ESS) surfaces identified coherent hotspots, with sparsely vegetated zones and marginal croplands showing systematically higher stress under rising temperature, increasing climatic water deficit, and declining moisture availability. The right-skewed stress distributions in natural cover types underscore asymmetric vulnerability to cumulative pressures in drylands.

Operationally, the workflow delivers policy-relevant, spatially explicit guidance for restoration planning, land-use regulation, and climate adaptation. It is scalable in Google Earth Engine, transferable to other drylands, and can be expanded through field validation and sensitivity analyses.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

AA-H: Writing – original draft, Conceptualization, Funding acquisition, Writing – review and editing, Supervision. ZI: Data curation, Methodology, Writing – review and editing, Software, Formal Analysis, Writing – original draft.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This research was funded by the Deanship of Scientific Research, King Saud University, through the Vice Deanship of Scientific Research Chairs: Chair of Climate Change, Environmental Development and Vegetation Cover.

Acknowledgements

The authors are grateful to the Deanship of Scientific Research, King Saud University, for funding through the Vice Deanship of Scientific Research Chairs: Chair of Climate Change, Environmental Development and Vegetation Cover.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Abbreviations

ANN, Artificial Neural Networks; GEE, Google Earth Engine; LST, Land Surface Temperature; NDVI, Normalized Difference Vegetation Index; NDWI, Normalized Difference Water Index (vegetation water content); NIR, Near-Infrared; SWIR – Shortwave Infrared.

References

Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., and Dean, J. (2016). TensorFlow: a system for large-scale machine learning. arXiv preprint arXiv:1605.08695. Available online at: https://arxiv.org/abs/1605.08695.

Google Scholar

Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., and Hegewisch, K. C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. Sci. Data 5, 170191. doi:10.1038/sdata.2017.191

PubMed Abstract | CrossRef Full Text | Google Scholar

Al-Areeq, A. M., Al-Zahrani, M. A., and Sharif, H. O. (2021). Physically-based, distributed hydrologic model for Makkah watershed using GPM satellite rainfall and ground rainfall stations. Geomatics, Nat. Hazards Risk 12, 1234–1257. doi:10.1080/19475705.2021.1924873

CrossRef Full Text | Google Scholar

Al-Rowaily, A. M., Assaeed, S. A., Al-Khateeb, A. A., Al-Qarawi, A. A., and Arifi, F. S. A. (2016). Vegetation and condition of arid rangeland ecosystem in central Saudi Arabia. Saudi J. Biol. Sci. 25 (6), 1022–1026. doi:10.1016/j.sjbs.2016.12.015

PubMed Abstract | CrossRef Full Text | Google Scholar

Alqadhi, S., Bindajam, A. A., Mallick, J., Talukdar, S., and Rahman, A. (2024). Applying deep learning to manage urban ecosystems in arid Abha, Saudi Arabia: remote sensing-based modelling for ecological condition assessment and decision-making. Heliyon 10, e25731. doi:10.1016/j.heliyon.2024.e25731

PubMed Abstract | CrossRef Full Text | Google Scholar

Banerjee, S., Nandi, T., Sati, V. P., Mezlini, W., Alkhuraiji, W. S., Al-Halbouni, D., et al. (2025). Integrating remote sensing, landscape metrics, and random forest algorithm to analyze crop patterns, factors, diversity, and fragmentation in a Kharif agricultural landscape. Land 14, 1203. doi:10.3390/land14061203

CrossRef Full Text | Google Scholar

Berdyyev, A., Al-Masnay, Y. A., Juliev, M., and Abuduwaili, J. (2024). Desertification monitoring using machine learning techniques with multiple indicators derived from Sentinel-2 in Turkmenistan. Remote Sens. 16 (23), 4525. doi:10.3390/rs16234525

CrossRef Full Text | Google Scholar

Bisong, E. (2019). “Google Colaboratory,” in Building machine learning and deep learning models on google cloud platform (Berkeley: Apress), 59–64. doi:10.1007/978-1-4842-4470-8_7

CrossRef Full Text | Google Scholar

Bürgi, M., Östlund, L., and Mladenoff, D. J. (2017). Legacy effects of human land use: ecosystems as time-lagged systems. Ecosystems 20, 94–103. doi:10.1007/s10021-016-0051-6

CrossRef Full Text | Google Scholar

Chai, T., and Draxler, R. R. (2014). Root mean square error (RMSE) or mean absolute error (MAE)? Arguments against avoiding RMSE in the literature. Geosci. Model Dev. 7, 1247–1250. doi:10.5194/gmd-7-1247-2014

CrossRef Full Text | Google Scholar

Chander, G., Markham, B. L., and Helder, D. L. (2009). Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sens. Environ. 113, 893–903. doi:10.1016/j.rse.2009.01.007

CrossRef Full Text | Google Scholar

de Beurs, K. M., and Henebry, G. M. (2005). A statistical framework for the analysis of long image time series. Int. J. Remote Sens. 26 (8), 1551–1573. doi:10.1080/01431160512331326657

CrossRef Full Text | Google Scholar

D’Odorico, P., Bhattachan, A., Davis, K. F., Ravi, S., and Runyan, C. W. (2013). Global desertification: drivers and feedbacks. Adv. Water Resour. 51, 326–344. doi:10.1016/j.advwatres.2012.01.013

CrossRef Full Text | Google Scholar

El Najjar, P., Chidiac, S., Probst, J.-L., El Omari, K., Ouaini, N., and El Azzi, D. (2023). Geochemical signature of the bed sediments at the outlet of the Ibrahim river (Lebanon): temporal variation. Environ. Monit. Assess. 195, 11103. doi:10.1007/s10661-023-11103-1

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Rawy, M., Batelaan, O., Al-Arifi, N., Alotaibi, A., Abdalla, F., and Gabr, M. (2023). Climate change impacts on water resources in arid and semi-arid regions: a case study in Saudi Arabia. Water 15 (3), 606. doi:10.3390/w15030606

CrossRef Full Text | Google Scholar

Ensor, J. E., Wennström, P., Bhattarai, A., Nightingale, A. J., Eriksen, S., and Sillmann, J. (2019). Asking the right questions in adaptation research and practice: seeing beyond climate impacts in rural Nepal. Environ. Sci. & Policy 94, 227–236. doi:10.1016/j.envsci.2019.01.013

CrossRef Full Text | Google Scholar

ESA (2022). Land cover CCI product user guide v2. Tech. Report. Available online at: https://maps.elie.ucl.ac.be/CCI/viewer/download/ESACCI-LC-Ph2-PUGv2_2.0.pdf (Accessed May 16, 2025).

Google Scholar

Etensa, T., Alemu, T., and Yayo, M. (2025). Rethinking the measurements and predictors of environmental degradation in Ethiopia: predicting long-term impacts using a kernel-based machine learning approach. Environ. Sustain. Indic. 25, 100554. doi:10.1016/j.indic.2024.100554

CrossRef Full Text | Google Scholar

Gao, B. C. (1996). NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 58, 257–266. doi:10.1016/S0034-4257(96)00067-3

CrossRef Full Text | Google Scholar

Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D., and Moore, R. (2017). Google Earth engine: planetary-scale geospatial analysis for everyone. Remote Sens. Environ. 202, 18–27. doi:10.1016/j.rse.2017.06.031

CrossRef Full Text | Google Scholar

Harper, K. L., Lamarche, C., Hartley, A., Peylin, P., Ottlé, C., Bastrikov, V., et al. (2023). A 29-year time series of annual 300 m plant-functional-type maps for climate models. Earth Syst. Sci. Data 15, 1465–1499. doi:10.5194/essd-15-1465-2023

CrossRef Full Text | Google Scholar

Hasan, S. S., Alharbi, O. A., Alqurashi, A. F., and Fahil, A. S. (2024). Assessment of desertification dynamics in arid coastal areas by integrating remote sensing data and statistical techniques. Sustainability 16, 4527. doi:10.3390/su16114527

CrossRef Full Text | Google Scholar

Hijmans, R. J. (2024a). Raster: geographic data analysis and modeling. R. Package Version 3, 6–26. Available online at: https://cran.r-project.org/package=raster.

Google Scholar

Hijmans, R. J. (2024b). Terra: spatial data analysis. R. Package Version 1, 7–65. Available online at: https://cran.r-project.org/package=terra.

Google Scholar

Huang, J., Yu, H., Guan, X., Wang, G., and Guo, R. (2016). Accelerated dryland expansion under climate change. Nat. Clim. Change 6, 166–171. doi:10.1038/nclimate2837

CrossRef Full Text | Google Scholar

IPCC (2023). Climate change 2023: impacts, adaptation and Vulnerability. Contribution of working group II to the sixth assessment report. Cambridge: Cambridge University Press. doi:10.1017/9781009325844

CrossRef Full Text | Google Scholar

Jiménez-Muñoz, J. C., and Sobrino, J. A. (2003). A generalized single-channel method for retrieving land surface temperature from remote sensing data. J. Geophys. Res. 108, 4688. doi:10.1029/2003JD003480

CrossRef Full Text | Google Scholar

Kamilaris, A., and Prenafeta-Boldú, F. X. (2018). Deep learning in agriculture: a survey. Comput. Electron. Agric. 147, 70–90. doi:10.1016/j.compag.2018.02.016

CrossRef Full Text | Google Scholar

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., et al. (2017). Climatologies at high resolution for the Earth’s land surface areas. Sci. Data 4, 170122. doi:10.1038/sdata.2017.122

PubMed Abstract | CrossRef Full Text | Google Scholar

Kendall, M. G. (1975). Rank correlation methods. 4th edn. London: Griffin.

Google Scholar

Kladny, K. R., Milanta, M., Mraz, O., Hufkens, K., and Stocker, B. D. (2024). Enhanced prediction of vegetation responses to extreme drought using deep learning and Earth observation data. Ecol. Inf. 80, 102474. doi:10.1016/j.ecoinf.2024.102474

CrossRef Full Text | Google Scholar

Kuhn, M. (2008). Building predictive models in R using the caret package. J. Stat. Softw. 28, 1–26. doi:10.18637/jss.v028.i05

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuhn, M., and Johnson, K. (2013). Applied predictive modeling. New York: Springer. doi:10.1007/978-1-4614-6849-3

CrossRef Full Text | Google Scholar

Kvalseth, T. O. (1985). Cautionary note about R². Am. Statistician 39, 279–285. doi:10.2307/2683704

CrossRef Full Text | Google Scholar

Lambin, E. F. (2001). Remote sensing and geographic information systems analysis. Amsterdam: Elsevier eBooks, 13150–13155. doi:10.1016/B0-08-043076-7/04200-5

CrossRef Full Text | Google Scholar

LeCun, Y., Bengio, Y., and Hinton, G. (2015). Deep learning. Nature 521, 436–444. doi:10.1038/nature14539

PubMed Abstract | CrossRef Full Text | Google Scholar

Lek, S., and Guégan, J. F. (1999). Artificial neural networks as a tool in ecological modelling: an introduction. Ecol. Model. 120, 65–73. doi:10.1016/S0304-3800(99)00092-7

CrossRef Full Text | Google Scholar

Ma, X., Geng, Y., Li, Y., Wu, Y., and Lu, Y. (2019). Comparative analysis of machine learning models for vegetation health assessment. ISPRS J. Photogrammetry Remote Sens. 149, 14–28. doi:10.1016/j.isprsjprs.2019.01.014

CrossRef Full Text | Google Scholar

Mann, H. B. (1945). Nonparametric tests against trend. Econometrica 13, 245–259. doi:10.2307/1907187

CrossRef Full Text | Google Scholar

Mas, J.-F., Velázquez, A., Díaz-Gallegos, J. R., Mayorga-Saucedo, R., Alcántara, C., Bocco, G., et al. (2004). Assessing land use/cover changes: a nationwide multidate spatial database for Mexico. Int. J. Appl. Earth Observation Geoinformation 5 (4), 249–261. doi:10.1016/j.jag.2004.06.002

CrossRef Full Text | Google Scholar

Olawade, D. B., Wada, O. Z., Ige, A. O., Egbewole, B. I., Olojo, A., and Oladapo, B. I. (2024). Artificial intelligence in environmental monitoring: advancements, challenges, and future directions. Hyg. Environ. Health Adv. 12, 100114. doi:10.1016/j.heha.2024.100114

CrossRef Full Text | Google Scholar

Paul, N., Sunil, G. C., Horvath, D., and Sun, X. (2025). Deep learning for plant stress detection: a comprehensive review of technologies, challenges, and future directions. Comput. Electron. Agric. 229, 109734. doi:10.1016/j.compag.2024.109734

CrossRef Full Text | Google Scholar

Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., and Grisel, O. (2011). Scikit-learn: machine learning in python. J. Mach. Learn. Res. 12, 2825–2830.

Google Scholar

Põldveer, E., Korjus, H., Kiviste, A., Kangur, A., Paluots, T., and Laarmann, D. (2019). Assessment of spatial stand structure of hemiboreal conifer-dominated forests according to different levels of naturalness. Ecol. Indic. 110, 105944. doi:10.1016/j.ecolind.2019.105944

CrossRef Full Text | Google Scholar

R Core Team (2024). R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. Available online at: https://www.R-project.org/.

CrossRef Full Text | Google Scholar

Ripley, B. D. (1996). Pattern recognition and neural networks. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511812651

CrossRef Full Text | Google Scholar

Roy, D. P., Wulder, M. A., Loveland, T. R., Woodcock, C. E., Allen, R. G., Anderson, M. C., et al. (2014). Landsat-8: science and product vision for terrestrial global change research. Remote Sens. Environ. 145, 154–172. doi:10.1016/j.rse.2014.02.001

CrossRef Full Text | Google Scholar

Said, Y., Barr, M., Saidani, T., and Atri, M. (2021). Desertification detection in Makkah region based on aerial images classification. Comput. Syst. Sci. Eng. 40, 607–618. doi:10.32604/csse.2022.018479

CrossRef Full Text | Google Scholar

Schmidhuber, J. (2015). Deep learning in neural networks: an overview. Neural Netw. 61, 85–117. doi:10.1016/j.neunet.2014.09.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, Y., Lao, D., Ruan, Y., Huang, C., and Xin, Q. (2023). A deep learning-based approach to predict large-scale dynamics of normalized difference vegetation index for the monitoring of vegetation activities and stresses using meteorological data. Sustainability. 15 (8), 6632–6632. doi:10.3390/su15086632

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

UNCCD (2022). Global land outlook: Second Edition (GLO2). United nations convention to combat desertification, Bonn. Available online at: https://www.unccd.int/resources/global-land-outlook/glo2

Google Scholar

U.S. Geological Survey (USGS) and National Aeronautics and Space Administration (NASA) (2022). Landsat collection 2 Level-2 surface reflectance products (TM, ETM+, OLI/TIRS). Available online at: https://www.usgs.gov/landsat-missions/landsat-collection-2.

Google Scholar

Vaswani, A., Shazeer, N., Parmar, J., Uszkoreit, J., Jones, L., and Gomez, A. N. (2017). Attention is all you need. Adv. Neural Inf. Process. Syst. 30, 5998–6008. doi:10.48550/arXiv.1706.03762

CrossRef Full Text | Google Scholar

Wang, F., Shao, W., Yu, H., Kan, G., He, X., Zhang, D., et al. (2020). Re-evaluation of the power of the mann-kendall test for detecting monotonic trends in hydrometeorological time series. Front. Earth Sci. 8, 14. doi:10.3389/feart.2020.00014

CrossRef Full Text | Google Scholar

Willmott, C. J., and Matsuura, K. (2005). Advantages of the mean absolute error (MAE) over the root mean square error (RMSE) in assessing average model performance. Clim. Res. 30, 79–82. doi:10.3354/cr030079

CrossRef Full Text | Google Scholar

Wulder, M. A., Masek, J. G., Cohen, W. B., Loveland, T. R., and Woodcock, C. E. (2012). Opening the archive: how free data has enabled the science and monitoring promise of landsat. Remote Sens. Environ. 122, 2–10. doi:10.1016/j.rse.2012.01.010

CrossRef Full Text | Google Scholar

Wulder, M. A., Roy, D. P., Radeloff, V. C., Loveland, T. R., Anderson, M. C., Johnson, D. M., et al. (2022). Fifty years of landsat science and impacts. Remote Sens. Environ. 280, 113195. doi:10.1016/j.rse.2022.113195

CrossRef Full Text | Google Scholar

Yang, J., Jiang, Y., Song, Q., Wang, Z., Hu, Y., Li, K., et al. (2025). An approach for multi-source land use and land cover data fusion considering spatial correlations. Remote Sens. 17, 1131. doi:10.3390/rs17071131

CrossRef Full Text | Google Scholar

Yue, S., Pilon, P., and Cavadias, G. (2002). Power of the Mann–Kendall and Spearman’s rho tests for detecting monotonic trends in hydrological series. J. Hydrology 259 (1-4), 254–271. doi:10.1016/s0022-1694(01)00594-7

CrossRef Full Text | Google Scholar