Dominant role of grazing and snow cover variability on vegetation shifts in the drylands of Kazakhstan

Drylands are critical ecological and biodiversity hotspots supporting pastoral livelihoods but are vulnerable to large-scale changes in vegetation structure and function. Socioeconomic shifts (e.g., cropland expansion, grassland degradation, herder sedentarization) and climate variability/change have altered dryland ecosystems, with global losses in ecosystem services estimated at ~US$10 trillion annually. Methodological variability has hindered robust detection and attribution of vegetation change drivers. In response, the IPCC detection-attribution frameworks emphasize separating detection (has change occurred?) from attribution (what are the causes?). Existing attribution studies often rely on single models or causal graphs that struggle with generalization, spatial heterogeneity, and multiple interdependent processes, and they typically underrepresent complex land-use practices due to limited spatial data. Kazakhstan (KZ) provides a pertinent test bed due to substantial anthropogenic impacts and policy shifts (e.g., Virgin Lands campaign, Aral Sea shrinkage, post-1991 cropland abandonment, institutional shifts, and ongoing climate stressors). This study asks: (1) What are the extent and spatial patterns of vegetation degradation hotspots across KZ? (2) What are the relative contributions of land use (LU), climate change (CC), climate variability (CV), and other factors (OF) to observed vegetation degradation? (3) Which SES drivers primarily contribute to vegetation degradation? The purpose is to produce a robust, scalable pixel-wise framework that disentangles natural versus anthropogenic drivers and quantifies their relative importance.

Prior studies document vegetation trends in Kazakhstan and Central Asia, but pixel/regional drivers remain uncertain due to differences in conceptualization, indicators, spatial/temporal scales, and data sources. Coarse-resolution products (e.g., GIMMS) have missed finer-scale degradation/greening patterns later observed by higher-resolution sensors. Many works emphasize climate change as the dominant driver, but critiques note that inadequate treatment of interannual climate variability (CV) biases attribution in short vegetation systems. Land-use practices (livestock grazing, cropland expansion/abandonment) are often rudimentarily represented due to sparse spatial data, downplaying direct human impacts. Formal detection-attribution standards from climate science have been suggested for land-use and biodiversity assessments. This study addresses these gaps using high-resolution datasets, explicit decomposition (LU vs CC vs CV vs OF), and causal/machine learning analyses that consider concurrent and lagged effects.

Study area: Kazakhstan (2.72 million km²), dominated by grasslands (~85.4%), with average annual precipitation 100–300 mm, subdivided into 17 regions and 215 districts; analyses focus on drylands (forests, water, wetlands, snow/ice, barren masked). Data: MODIS NDVI (MOD13Q1, 250 m, 16-day aggregated to monthly), LST (MOD11A2, 1 km), Snow cover fraction (MOD10A1, 500 m), Land cover (MCD12Q1, 500 m). Climate: CHELSA monthly precipitation, air temperature, VPD (1 km); TerraClimate shortwave radiation (SRad, 4 km resampled to 1 km). Socioeconomic: WorldPop population density (1 km), DMSP-OLS and VIIRS nighttime lights (1 km), gridded GDP (1 km), livestock densities (sheep & goats, horses; 1 km) downscaled from census data. CO2 concentration (CMIP5 RCP8.5). Monthly time step 2000–2019; data projected to Albers Equal Area; masking and resampling in GEE/ArcGIS. Framework (three stages): Stage 1 – Detection and decomposition (TSS-RESTREND). Per-pixel Theil–Sen slope and Spearman’s rho on NDVI to detect trends (2000–2019). Adjust NDVI to remove CO2 fertilization using a theoretical photosynthesis–CO2 relationship to generate NDVIadj. Compute vegetation–climate relationships (VCR) using accumulated precipitation and temperature (1980–2019), separating climate change (long-term trend; 20-year smoothing) from climate variability (detrended interannual variability). Derive contributions: LU (residuals from OLS of NDVIadj vs climate), CC (difference between climate-driven and CV-driven NDVI components), CV, and OF (remainder). Compute absolute slopes and identify per-pixel dominant factor; segregate LU-dominant degraded pixels from climate-driven (CC+CV) degraded pixels for subsequent analysis. Stage 2 – Causation and contribution (Granger causality, GC). For degraded pixels, compute anomalies and test stationarity (ADF). Fit bi-variate VAR models (optimal lag 1–6 months via AIC; 5-fold CV) for NDVI and each driver. Evaluate F-statistic and p-value; retain pixels with significant GC (p < 0.05). LU drivers: sheep & goats, horses, GDP, POPD, nighttime lights. Climate drivers: PCP, SRad, LST, VPD, snow cover. Stage 3 – Attribution (pixel-wise Random Forest with Shapley values). Build per-pixel RF models separately for LU-dominant degraded pixels (NDVI ~ five LU drivers) and climate-driven degraded pixels (NDVI ~ five climate drivers). Use 80/20 train-test split, 10-fold CV, hyperparameter tuning (trees 50–300; learning rate 0–1) to minimize RMSE. Compute R², RMSE, MSE; mask pixels with R² < 0.5. Use Shapley values to attribute dominant driver impacts and to identify hotspots. All GC/RF/Shapley in MATLAB; TSS-RESTREND in R; spatial operations in ArcGIS.

- Extent of change: 45.71% of Kazakhstan (out of 2.72 million km²) experienced vegetation degradation (browning). Considering the non-masked area (2.39 million km²), 52.02% degraded and 47.98% greening; 42.18% of total KZ showed greening. - Spatial patterns: Degradation hotspots in western and southern provinces (Caspian lowland desert, Central Asian northern desert parts of Aqtobe and Mangghystau, Kazakh steppe region of Aqtobe, Altai–Western Tian Shan steppe, Tian Shan foothill arid steppe, parts of Aqmola and Pavlodar). Greening primarily in semi-desert, steppe, upland steppe of Qaraghandy, North, and East KZ and parts of Almaty, Pavlodar, Zhambyl. - Decomposition of drivers (dominant factor by absolute slope): Overall vegetation change across KZ dominated by climate variability (CV) 41.3% (0.99 million km²), followed by land use (LU) 33.54% (0.81 million km²), other factors (OF) 15.45% (0.35 million km²), and climate change (CC) 9.71% (0.23 million km²). In degraded areas only: LU 22.54% (0.54 million km²) was dominant, followed by CV 21.93% (0.52 million km²), OF 4.38% (0.105 million km²), and CC 3.17% (0.076 million km²). - GC results (significant pixels within degraded areas): • Land-use dominated degraded pixels: Sheep & goats showed significant GC in 75.6% of pixels; horses 73.39%; nighttime lights 62.6%; population density 29.71%; GDP 21.94%. Small ruminants generally had negative effects on vegetation. • Climate-driven degraded pixels: Snow cover significant in 96.9%; LST 49.25%; PCP 40.35%; SRad 37.75%; VPD 18.81%. LST, SRad, and VPD generally had negative effects; precipitation positive overall; snow cover mixed spatially (positive in NW, negative in S). - Attribution via RF + Shapley (dominant driver shares among degraded pixels): • LU-dominant degraded pixels (~0.525 million km² modeled): Sheep & goat density dominated 11.78% of pixels, horses 8.06%, POPD 1.3%, GDP 0.97%, nighttime light 0.76%. Small ruminants dominated impacts in the Kazakh steppe, Pontic steppe (West), and upland steppe/semi-desert in Qaraghandy and East KZ. • Climate-driven degraded pixels (~0.588 million km² modeled): Snow cover dominated 23.3%, precipitation 14.87%, SRad 13.37%, LST 4.58%, VPD 4.02%. Snow cover dominance spanned 8 of 14 provinces; precipitation dominated parts of South KZ and western/northern steppe; SRad dominated parts of West KZ, Aqtobe, Qostanay, and South KZ. - Overall attribution coverage: Dominant SES driver hotspots identified for 41.57% of degraded pixels (~1.14 million km² degraded area basis), corresponding to 19.81% of KZ (out of 2.39 million km² non-masked area).

The framework successfully addresses the research questions by separating detection of vegetation trends from attribution to LU, CC, CV, and OF at pixel scale. Results highlight LU as the leading contributor to degradation in KZ’s drylands, particularly via grazing pressure from small ruminants, aligning with known post-Soviet socioecological transitions and livestock rebounds. Climate variability exerts the largest influence on overall vegetation dynamics, underscoring dryland sensitivity to interannual hydroclimatic fluctuations. Snow cover emerges as a pivotal climatic driver: increases or altered phenology can trigger early greening followed by late-season water stress, while declines in snow cover reduce cold-season moisture storage, both contributing to degradation depending on region. Precipitation deficits in western/northwestern KZ and warming amplify stress, while population/GDP and industrial activity have localized effects (e.g., Ustyurt Plateau). The findings delineate where land management (grazing systems, livestock numbers) versus climate adaptation (water and drought management) could be most effective, thus informing strategies toward SDG 15.3 (land degradation neutrality).

This study advances detection and attribution of dryland vegetation change by implementing a high-resolution, multi-stage framework that (1) detects pixel-wise NDVI trends, (2) decomposes observed change into LU, CC, CV, and OF, and (3) attributes dominant SES drivers using causal and machine learning approaches. In Kazakhstan, 45.71% of the landscape shows degradation, with LU dominating degraded areas and CV dominating overall variability. Small ruminant grazing pressure and snow cover variability are the primary LU and climate drivers, respectively, with additional roles for precipitation, solar radiation, and temperature. The framework identifies degradation hotspots and dominant drivers over 19.81% of KZ, providing actionable insights for targeted restoration and management (e.g., rotational/planned grazing, pastoral mobility, adjusting stocking densities, climate adaptation in water-limited regions). Future research could (i) apply the framework to other drylands, (ii) incorporate higher-resolution and additional vegetation indicators (e.g., soil-sensitive indices), (iii) expand multi-driver interactions and feedbacks, and (iv) refine uncertainty treatment (e.g., multiple testing adjustments, spatial autocorrelation controls).

- Data and proxies: NDVI used as a vegetation proxy may be confounded by shrub encroachment, agricultural intensification, non-native species, and biological soil crusts. Alternative indices (soil-sensitive) and higher-resolution sensors (Landsat/Sentinel) could improve detection. - Dataset uncertainties: Input climatic/demographic datasets carry uncertainties; SRad was resampled from 4 km; GDP held constant for 2016–2019 due to lack of gridded data. Although prior inter-dataset sensitivity suggests <4% average change for some combinations, broader sensitivity tests (e.g., multiple NDVI/LST products) are warranted. - Spatial resolution and scaling: Trends detected at 250 m were resampled to 1 km for GC/RF, introducing potential mixed-pixel effects and spatial autocorrelation. - Statistical considerations: Spatial autocorrelation and multiple hypothesis testing may increase false positives; future work could employ false discovery rate controls (Benjamini–Hochberg) or permutation tests. - Causality scope: GC identifies predictive causality but not physical causation; bidirectional effects (e.g., vegetation influencing land use/climate variables) were not tested. - Model filtering: Pixels with non-significant GC or RF R² < 0.5 were excluded, limiting attribution coverage but enhancing robustness of reported hotspots.