Causality-Constrained XGBoost–SHAP Reveals Nonlinear Drivers and Thresholds of kNDVI Greening on the Loess Plateau (2000–2019)

The Loess Plateau has experienced persistent vegetation greening over the past two decades, yet this recovery has occurred under a concurrent intensification of atmospheric evaporative demand and drying. This raises a key land–atmosphere question: which hydroclimatic processes most strongly constrain greening, and where do vegetation responses shift across environmental regimes? To address this issue, we integrated spatiotemporal trend analysis, Geographical Convergent Cross Mapping (GCCM)-based directional attribution, and an interpretable machine-learning framework combining Extreme Gradient Boosting (XGBoost) and SHapley Additive exPlanations (SHAP) to diagnose the dominant controls and threshold-like response patterns of vegetation activity. Using 1 km kernel Normalized Difference Vegetation Index (kNDVI) and eight hydroclimatic variables during 2000–2019, we found that regionally averaged kNDVI increased from 0.099 in 2000 to 0.164 in 2019, with a significant trend of 0.003 year−1, and greening trends covered 65.503% of the Loess Plateau. Over the same period, Vapor Pressure Deficit (VPD) increased from 0.142 to 0.275 kPa (+0.133 kPa), indicating that vegetation recovery did not occur under a more humid atmospheric background. GCCM results consistently showed stronger directional influence from hydroclimatic drivers to kNDVI than the reverse, with evaporation and thermal conditions, especially Tmin, emerging as the dominant constraints, followed by Tmax, VPD, and wind speed, whereas precipitation showed comparatively weaker recoverable influence. The tuned XGBoost model achieved strong out-of-sample performance (R2 = 0.9611, RMSE = 0.0188, MAE = 0.0131), and SHAP revealed clear nonlinear thresholds: evaporation and Tmin shifted into persistently positive contribution regimes beyond 302 mm and −17.6 °C, respectively; Tmax became predominantly inhibitory beyond −1.9 °C, and Palmer Drought Severity Index (PDSI) exhibited a multi-stage non-monotonic transition around −0.7. These results provide a coherent evidence chain linking directional influence, relative contribution, and threshold boundaries, offering quantitative support for identifying climate-sensitive zones and restoration risk regimes under continued warming and rising atmospheric dryness.