Nonlinear Associations and Candidate Statistical Thresholds Between Ecosystem Services and Groundwater Storage Anomalies: An Analysis Using CatBoost–SHAP

Ecological restoration in arid regions can enhance ecosystem functions but may coincide with groundwater stress. Using nine headwater basins of the Tarim River Basin, we analyzed a GRACE/GLDAS-derived groundwater storage anomaly (GWSA) proxy for 2003–2025 and four ecosystem services and the integrated index (ESI) for 2000–2025; joint analyses covered 2003–2025. Correlation analysis, joint-trend zoning, and CatBoost–SHAP characterized spatial mismatches and nonlinear associations. The GWSA proxy declined by −0.31 mm month−1, with negative anomalies intensifying after 2020 and a northwest-decrease/southeast-increase pattern. Comparisons with groundwater-level anomalies showed moderate agreement at Aksu (r = 0.52) but weak, nonsignificant agreement at Cele (r = −0.08), indicating spatially heterogeneous reliability. Carbon storage, habitat quality, net primary productivity, and ESI increased, whereas sand retention declined. Areas with increasing ESI but decreasing GWSA accounted for 42.5% of the study area. Five CatBoost models achieved R2 values of 0.868–0.977. Regarding ESI, SHAP identified a candidate transition near −18.6 mm, below which GWSA was associated with negative model contributions; human-footprint contributions became increasingly negative with disturbance. These statistical associations, rather than causal effects, demonstrate that ecosystem-service improvement does not necessarily indicate groundwater recovery.