The Tianshan Mountains form the central mountain system of Xinjiang Uygur Autonomous Region, China (39°24′40″-45°23′08″N, 73°50′28″-95°33′56″E). This study focused on the Chinese sector, spanning 11 prefecture-level administrative units and covering about 37.02×10⁴ km² (Fig. 1). The range has an average ridgeline elevation of about 4000 m, and Tomur Peak reaches 7443 m a.s.l. (Lu et al., 2022). The Tianshan Mountains form the natural and climatic divide between the Junggar Basin to the north and the Tarim Basin to the south, bordering the Taklimakan Desert and the Gurbantunggut Desert, and representing a typical mountain-oasis-desert complex system (Wang et al., 2016; Zhou and Lei, 2018). Climate conditions exhibit pronounced north-south contrasts: the annual mean temperature ranges from 2.50°C-5.00°C on the northern slope to 7.50°C-10.00°C on the southern slope (Li et al., 2019). Water scarcity is a defining feature, with strong spatial precipitation gradients: on the northern slope, annual precipitation decreases from approximately 250.00-350.00 mm in the western sector to 150.00-200.00 mm in the central sector, while the southern slope generally receives less precipitation than the northern slope (Wang et al., 2016). In the adjacent Tarim Basin, annual total precipitation is typically less than 50.00 mm, highlighting the hyper-arid background of the region (Zhou and Lei, 2018). As a "water tower of Central Asia", the Tianshan Mountains supply the headwaters of major regional rivers through combined contributions of glacier/snow meltwater and mountain precipitation (Chen et al., 2016). Their strong elevation gradients produce clear vertical zonation that supports diverse vegetation belts and ecosystem functions relevant to NPP dynamics in arid mountains (Li et al., 2019).

This study assembled and preprocessed multi-source datasets to support the integrated analytical framework (Table 1). Core environmental variables included land use, the normalized difference vegetation index (NDVI), climatic and topographic factors, which served as the primary inputs to CASA model for estimating NPP. Monthly remote-sensing and meteorological data were used to drive CASA and to generate annual composites and cumulative values for the period 2000-2020. Anthropogenic pressure indicators included population density, gross domestic product (GDP) density, nighttime light index, human footprint, and grazing intensity, and accessibility metrics included distances to roads, highways, water bodies, and government centers; these variables were compiled to quantify anthropogenic pressure and infrastructure accessibility influences on NPP. Spatial consistency across datasets was ensured through reprojection, resampling to a 1-km resolution, and grid alignment. To refine the hydrothermal fields used by CASA, we complemented TerraClimate data with higher-resolution temperature and precipitation from the National Tibetan Plateau Data Center (NTPDC; https://data.tpdc.ac.cn/home). To represent land-cover classification uncertainty and recent changes, we jointly used two land-cover products: the Land Use and Cover Change (LUCC) dataset provided by Resource and Environmental Science Data Platform (http://www.resdc.cn) and the China Land Cover Dataset (CLCD) (https://zenodo.org/records/15853565). In addition, a 2024 data snapshot was compiled for supplementary near-term assessments and ongoing monitoring only and was not included in the core 2000-2020 analyses. All supplementary datasets used for robustness checks were harmonized to match the projection, spatial resolution, and grid of the core dataset.

This study estimated annual NPP during 2000-2020 at the pixel scale using CASA model, a remote-sensing light use efficiency model widely applied for regional and global NPP estimation (Zhou et al., 2025). Evaluations across diverse ecosystems indicated that CASA-based estimates showed good generality and stability (Zeng et al., 2024; Dong et al., 2025). In CASA, NPP was calculated as the product of absorbed photosynthetically active radiation (APAR) and actual light use efficiency (ε):

where NPP_(x,t) is the NPP at pixel x in month t (g C/(m²•month)); APAR_(x,t) is canopy absorbed photosynthetically active radiation (MJ/(m²•month)); ε_(x,t) is the actual light use efficiency (g C/MJ); SOL_(x,t) is total incoming shortwave solar radiation derived from the TerraClimate shortwave radiation product (MJ/(m²•month)); FPAR_(x,t) is the fraction of absorbed photosynthetically active radiation; 0.5 represents the proportion of photosynthetically active radiation in total shortwave radiation; ε_max is the maximum light use efficiency (g C/MJ) (Zhu et al., 2006); and T_ε1(x,t), T_ε2(x,t), and W_ε(x,t) are the temperature- and moisture-stress scalars used to down-regulate ε_max, with T_ε1(x,t) and T_ε2(x,t) accounting for temperature constraints via complementary response functions and W_ε(x,t) representing water limitation (Zhu et al., 2007). The formulations of T_ε1(x,t), T_ε2(x,t), and W_ε(x,t) have been applied and validated in arid and mountainous environments similar to the Tianshan Mountains (Liu et al., 2022; Heng and Wang, 2023; Zhang et al., 2024). Monthly NPP was summed to obtain annual NPP for each pixel for 2000-2020.

Field measurements were used to derive plot-based NPP estimates as independent references for evaluating CASA-derived NPP. A total of 265 vegetation sampling plots were established across the Tianshan Mountains in July 2018 and July 2024, stratified by elevation and ecosystem type (grassland and Picea schrenkiana Fisch. & C. A. Mey. forest). For grasslands, aboveground vegetation was harvested from 1 m by 1 m quadrats, oven-dried and plot-based grassland NPP was calculated as:

where NPP_g is the plot-based annual grassland NPP (g C/(m²•a)); AGB is the aboveground dry biomass (g dry matter/m²); and k is the biomass-to-carbon conversion coefficient and was set to 0.45 in this study (Abdi et al., 2016; Wu et al., 2022).

For P. schrenkiana dominated forests, 10 m×10 m plots were established to quantify stand structure, including tree density, diameter at breast height, and tree height, and understory vegetation was surveyed in nested 1 m×1 m quadrats. Forest NPP was estimated as the sum of annual community growth and annual litter input (Wang et al., 2001).

where G is the annual community growth (kg C/(m²•a)); L is annual litter input (kg C/(m²•a)); NPP_f is the annual forest NPP (kg C/(m²•a)); B is the stand biomass per unit area (kg dry matter/m²); T is the stand age (a); and c, d, e, and f are the forest-type constants. For spruce forests, c=0.4598, d=0.0069, e=10.1320, and f=0.0874 (Wang et al., 2009).

In MCD12Q1, "urban and built-up" was defined as pixels with at least 30.00% impervious surface and was termed "construction land" in this study. CASA-derived NPP over construction land therefore represented the productivity of vegetation within mixed pixels rather than impervious surfaces, which is essential for interpreting NPP patterns in urbanized areas.

To document uncertainty in key model inputs, we conducted consistency diagnostics for the meteorological forcing and land-cover information. Temperature and precipitation fields from TerraClimate and the NTPDC were compared in representative years (2000, 2005, 2010, 2015, and 2020) using correlation- and error-based metrics, and pixel-wise differences for 2020 were mapped. Land-cover products were also cross-compared for Built% and Crop%, and grids with absolute inter-product differences of at least 5.00% were identified as hotspots. To quantify the agreement between the TerraClimate and NTPDC datasets, we employed a suite of statistical metrics: the Pearson correlation coefficient and the Spearman's rank correlation coefficient to assess linear and monotonic relationships, respectively; the root mean square error (RMSE) and the median absolute error (MdAE) to measure the magnitude of average and typical deviations; and the mean bias error (MBE) to evaluate systematic over- or under-estimation. Detailed diagnostics are provided in the Supplementary Material (Tables S1 and S2; Figs. S1 and S2).

Based on the pixel-level annual NPP product for 2000-2020, this study quantified spatiotemporal variability of NPP by integrating interannual dynamics, spatial patterns, stability, and trend significance. Interannual dynamics were summarized at the regional scale using total NPP and area-weighted mean NPP. Spatial patterns were characterized using multi-year mean NPP maps and zonal statistics across major land-cover types and environmental strata.

Interannual stability at the pixel level was quantified using the CV (Hou et al., 2015; Wei et al., 2022):

where CV_x is the coefficient of variation of annual NPP at pixel x; σ_x is the standard deviation of annual NPP at pixel x during 2000–2020 (g C/(m²•a)); and $ {\overline{\mathrm{NPP}_{x}}} $ is the corresponding multi-year mean annual NPP at pixel x (g C/(m²•a)). Long-term trends were estimated using the Theil–Sen median slope (Eastman et al., 2013; Tian et al., 2015):

where Slope_x is the Theil-Sen slope of annual NPP at pixel x (g C/(m²•a)); NPP_(x,i) and NPP_(x,j) are the annual NPP values at pixel x in years i and j, respectively (g C/(m²•a)); and n is the number of years from 2000 to 2020, so here n=21.

Trend significance was evaluated using the Mann-Kendall test. The Mann-Kendall statistic S_x for pixel x was computed as (Hamed, 2008):

When tied values occurred (i.e., two or more years had identical NPP values within the time series at a same pixel), the variance of S_x (Var(S_x)) was calculated with tie correction (Blain, 2013):

where g is the number of tied groups; and t_k is the size of the k^th tied group.

The standardized test statistic Z_x was then obtained as (Gocic and Trajkovic, 2013):

The two-sided P-value (P_x) was derived from Z_x under the standard normal distribution (Sonali and Nagesh Kumar, 2013):

where Ф is the cumulative distribution function of the standard normal distribution. Following the criteria summarized in Table 2, P_x was used together with slope_x to assign trend categories.

To assess the long-term persistence and future tendencies of NPP trends, we employed the Hurst exponent (H), a robust non-parametric method derived from rescaled range (R/S) analysis. The H value quantifies the autocorrelation or memory effect within a time series, indicating whether an observed trend is likely to persist, reverse, or behave randomly in the future (Bui and Ślepaczuk, 2022).

For a given annual NPP time series $\left\{N P P_{u}\right\}_{u=1}^{n}$ at a pixel (here n=21 for 2000–2020), the Hurst exponent (H) was estimated using the rescaled range (R/S) approach. For each segment length m (2≤m≤n; i.e., the number of annual observations in a subseries), the cumulative deviate series X_t,m was first computed as follows (Hamed, 2007):

where $\overline{\mathrm{NPP}}_{m}$ is the mean NPP of the same m-year subseries; u is the index in the summation; and t denotes the cumulative time step within the subseries. The range R_m and standard deviation S_m for each segment length m were then defined as (Mielniczuk and Wojdyłło, 2007):

The rescaled range (R/S)_m was calculated as (R/S)_m=R_m/S_m. According to R/S theory, log(R/S)_m exhibits a linear relationship with log(m) (Hamed, 2007):

where a is the intercept. The Hurst exponent (H) was estimated as the slope of the ordinary least squares regression fitted to the log(m)-log(R/S)_m relationship.

The Hurst exponent (H) indicates the persistence of a time series. Within an ecological context, this study established a comprehensive analytical framework: (1) determining whether the historical NPP trend represents degradation or improvement based on the direction of its slope; and (2) assessing the likely future trajectory by incorporating the persistence characteristics suggested by the H value. By integrating these two dimensions, 9 distinct categories describing future NPP dynamics are derived, as detailed in Table 3.

Fifteen variables capturing natural geography, climate, and anthropogenic pressure were used as drivers of NPP in this study. These include X1 (land use), X2 (temperature), X3 (precipitation), X4 (digital elevation model, DEM), X5 (slope), X6 (aspect), X7 (population density), X8 (GDP density),X9 (nighttime light index), X10 (human footprint), X11 (grazing intensity), X12 (distance to road), X13 (distance to highway), X14 (distance to water body), and X15 (distance to government center). These variables encompassed factors such as climate (temperature and precipitation), topography (elevation, slope, and aspect), human activities (population density, economic activity, grazing intensity, and accessibility), and land use. Using these variables, XGBoost modeled the nonlinear responses of NPP to multiple factors, and SHAP improved interpretability. XGBoost models complex nonlinear relationships by iteratively building an ensemble of decision trees, and its regularized objective reduces overfitting (Ma et al., 2025b). SHAP, grounded in Shapley values, attributes each variable's marginal contribution across feature combinations and quantifies both effect size and direction (Wu et al., 2024b):

where y is the model-predicted NPP for a given sample; ϕ₀ is the expected (baseline) model output, i.e., the mean prediction of NPP; ϕ_k is the SHAP value representing the contribution of the k^th feature to y; M is the number of input features (here M=15, corresponding to X1-X15).

Locally estimated scatterplot smoothing (LOESS) regression smoothed the SHAP dependency plots of major variables to explore their nonlinear responses and threshold effects on NPP. This method fits low-order polynomials to locally weighted samples without requiring predefined functions, effectively capturing variable trends and identifying critical intervals (Austin and Steyerberg, 2013). As a local fitting method, LOESS relies on data's local structure, necessitating a large and densely sampled dataset for stability and reliability. A 2 km×2 km grid with 90,348 cells was created in ArcGIS Pro to match the characteristics of LOESS fitting. NPP estimates and values of 15 drivers were assigned to each grid-cell center, enabling fine-scale analysis of driving mechanisms.

Land-use change in the Tianshan Mountains from 2000 to 2020 was assessed with a seven-category transfer matrix, and land flows and transfer magnitudes were visualized (Fig. 3). From 2000 to 2010, grassland increased to 43.63% (+3.02%), mainly from unutilized land (14,295.99 km²), cropland (1871.76 km²), and wetland (1789.30 km²). Unutilized land decreased to 48.06%, mainly converting to grassland, cropland, and other ecological land. Cropland increased to 7.02%, with a net gain of 24.5% mainly from grassland. Wetland area shrank dramatically by 75.12% to 0.16% of the study area, with 1789.30 km² converted to grassland, while forest area decreased by 45.65% to 0.28% of the study area. Meanwhile, water bodies decreased slightly from 0.48% to 0.43%, and construction land remained approximately 0.42%. From 2010 to 2020, ecological restoration and structural adjustment became evident. Forest rebounded by 20.20% and cropland expanded by 9.08%, largely at the expense of grassland (6739.33 km²). Grassland expansion (+2.34%) was driven primarily by conversion from unutilized land, which declined by 3.80% to 46.23%. Ecological recovery included wetland rebounding by 50.45% to 0.24% and water bodies increasing to 0.46%, while construction land remained near 0.42% with an increase of 0.84%. Overall, this period combined ecological recovery with cropland intensification.

From 2000 to 2020, cropland increased to 7.65% (+35.80%), forest decreased to 0.34% (-34.67%), grassland reached 44.65%, and unutilized land declined to 46.23%. These changes informed the subsequent NPP analysis, as gains in grassland and cropland corresponded to higher NPP areas, whereas unutilized and degraded grasslands aligned with low-NPP zones. Percentages of construction land (Built%) and cropland (Crop%) were highly consistent among the land use products (MCD12Q1, LUCC, and CLCD), with points clustering along the 1:1 line and most deviations within ±5.00% (Fig. S1). Cross-product comparisons showed that Built% and Crop% were broadly consistent among land-cover datasets at the 6 km grid scale (Table S1), and most grid cells clustered around the 1:1 line within the ±5.00% band (Fig. S1). Nevertheless, 3.03%-18.41% of grids exceeded the ±5.00% threshold across product pairs (Table S1), indicating localized discrepancies in derived land-use proportions and thus spatially heterogeneous uncertainty.

Combining Hurst exponent (Fig. 7a) and slope (Fig. 6c), NPP trends in the Tianshan Mountains were classified into nine types, revealing marked spatial heterogeneity and distinct transition patterns (Table 4; Fig. 7b). NPP was projected to decline across much of the Tianshan Mountains, with SPD (36.41%) and WPD (17.52%) collectively covering 53.93% of the study area, primarily in central and northern areas, indicating a significant and spatially continuous degradation trend. WAPI at 13.63% and WAPD at 12.81% occurred mainly in central and western Tianshan Mountains, indicating reversals relative to historical trends with moderate dynamics. Improving trends were limited: SPI covered 4.97% and SAPI covered 0.36%, mainly in southern Tianshan Mountains, indicating low improvement potential. Regions classified as NSC covered 12.83%, were spatially scattered, and showed limited future variation. SAPD at 1.39% and MAPI at 0.07% comprised a very small share and were highly localized. Overall, persistent degradation was likely to dominate future NPP in the Tianshan Mountains, with localized and weak improvement and a strongly heterogeneous spatial pattern.

Integrated diagnostics confirmed that persistent degradation would affect over half of the region, concentrated in specific transition zones. This trend was primarily driven by the dominant role of land use, coupled with critical hydroclimatic thresholds and amplifying anthropogenic pressures. Complex nonlinear interactions among these drivers further intensified the degradation risks in identified hotspots. At the 6 km scale, we cross-validated two consistent four-class disturbance schemes with an identical type definition (I-IV; Table S3): the human disturbance four-type hydrothermal-regime clustering (HD4) and an independent rule-based disturbance classification. HD4 delineates the four types via clustering based on hydrothermal conditions and human-pressure proxies (Built%, Crop%, grazing intensity, population density, and GDP), whereas the rule-based scheme assigns the same four types using predefined thresholds of these indicators. The two schemes showed strong agreement (overall accuracy=83.79% and Kappa coefficient=0.71; Table S3). Most mismatches were mainly concentrated in mixed transition cells between types II and III and between types I and IV (Fig. S3).

The 2024 diagnostics revealed synergy hotspots along oasis-desert ecotones and urban fringes. Threshold screening identified high-risk grids covering 3.88% of the study area (Fig. S4), with widespread low NPP values along arid margins. Across all 2 km grid cells in the study area, the 2024 mean values were 202.31 g C/(m²•a) for NPP, 269.05 mm for precipitation, 3.53°C for temperature, 0.51% for Built%, and 6.51% for Crop%. The spatial pattern aligned with SHAP-based thresholds, indicating that hydrothermal stress combined with anthropogenic pressure governs near-term NPP risk.

NPP, a key indicator of terrestrial carbon cycling, is highly sensitive to land use change, reflecting the combined effects of human and natural processes on ecosystem structure and function (Zeng et al., 2025). In the Tianshan Mountains, the overall NPP increase from 2000 to 2020 was closely linked to land use changes, primarily driven by grassland and cropland expansion under favorable mid- to low-elevation hydrothermal conditions (Liu et al., 2021). This expansion enhanced regional carbon fixation through higher photosynthetic efficiency and biomass accumulation. This pattern aligned with Hou (2020), who reported that cropland and grassland expansion dominated NPP increase in the Weigan-Kuqa River Basin, Xinjiang Uygur Autonomous Region, China. Converting unutilized land to ecological types like grassland optimizes ecological space, enhancing soil-water conservation and carbon sequestration (Hu et al., 2025). This positive link between intensive land management and enhanced ecosystem services is common in arid western China (Guo and Zhang, 2021).

Despite localized recovery under return-to-forest programs, forest did not offset prior degradation. Due to limited forest area and slow growth rates, forest's short-term contribution remained constrained (Pang et al., 2024). Similarly, wetland exhibited high per-area NPP due to moisture availability and national greening efforts, despite their limited spatial extent. However, both forests and wetlands experienced net declines over the study period, weakening their overall regional NPP support (Xiang et al., 2023; Ren et al., 2025). Construction land expanded slowly but reduced NPP through surface sealing and heat-island effects, especially along urban fringes (Niu et al., 2024). Note that NPP for construction land does not originate from impervious surfaces, but reflects the productivity of vegetation components (e.g., urban parks, lawns, and street trees) within these mixed pixels, as defined by our land cover classification (Liu et al., 2024). Consequently, land use management should prioritize optimizing ecological space and protecting high-value ecosystems to sustain productivity and carbon sequestration.

NPP in the Tianshan Mountains resulted from a complex coupling of human activity, land use, and climate. Among these factors, land use type explained the largest share of NPP variability. Ecological lands enhanced NPP by increasing vegetation cover and improving local conditions, while construction land reduces it through surface sealing and habitat fragmentation. The SHAP analysis correctly captured this suppressive effect. However, the non-zero NPP value assigned to construction land underscored that our model estimates reflect the productivity of the integrated pixel (Luo and Wang, 2009). This interpretation aligns with the land cover definition and highlights that the CASA model, while powerful, attributes productivity to the entire mixed pixel rather than isolating the null contribution of impervious surfaces. Unutilized land currently had low productivity due to poor soils but held significant potential for ecological restoration, aligning with evidence that ecological lands promote productivity whereas non-ecological lands constrain it (Li et al., 2023). The pattern of wetland NPP exceeding that of grassland in the Tianshan Mountains, driven by stable water supply alleviating moisture stress, underscored the primary control of precipitation (Chen et al., 2023). This key climatic driver exhibited an optimal window of 300.00-400.00 mm, alongside a temperature threshold near 0.00°C, which collectively confirmed the strong hydrothermal regulation in arid areas (Huxman et al., 2004). Anthropogenic pressures were nonlinear. Moderate population and GDP density promoted NPP, whereas high-density development imposed stress and reduced productivity. Grazing showed low-intensity promotion and high-intensity suppression, indicating grassland sensitivity and the benefits of moderate grazing (Suo et al., 2023). These controls aligned with national and regional evidence on climate and human effects (Wei et al., 2025). At the 6 km scale, robustness tests showed strong agreement between the human disturbance hydrothermal-regime clustering and the rule-based scheme. Consistent with the SHAP-based attribution, land use and hydroclimatic conditions dominated NPP variability, while anthropogenic pressure acted as an important modifier in hotspot areas.

Near-term diagnostics for 2024 linked mechanisms to management (Fig. S4). High-risk grids clustered along these belts, accounting for 3.88% of the study area effectively evaluated region. Priority actions should include water-saving irrigation and allocation, grazing control or rotation, and peri-urban greening and ecological corridors. The spatial concordance between high-risk areas and projected persistent-degradation zones supported using the identified thresholds and interaction mechanisms for early warning and targeted intervention. Unlike plain-type arid areas (e.g., Ningxia Hui Autonomous Region), degradation in the Tianshan Mountains clusters along ecotones and peri-urban fringes under orographic amplification (Cong et al., 2025). This reflected the co-occurrence of water or heat stress with anthropogenic pressure, consistent with known vegetation-specific drought responses and the dominant role of topography and climate in mountain basins (Wu et al., 2024a; Zhang et al., 2025).

CASA performed well in the Tianshan Mountains, consistent with Pan and Li (2015), who reported a high correlation (R²=0.870) for an improved version in the Northwest China. However, model-level uncertainties remain. In extremely arid zones, CASA can be biased due to the simplified moisture stress term (Bao et al., 2016). Transferring the ε_max parameter across regions may introduce magnitude uncertainty, necessitating future calibration for the Tianshan Mountains using site-level flux and plot data (Guan et al., 2022). CASA also tends to underestimate NPP compared with MODIS products and field measurements. Potential reasons include uncertainties in validation data (e.g., sampling bias toward high-productivity sites), the limited ability of MOD17A3 to capture vegetation-soil-moisture heterogeneity, and spatiotemporal scale mismatches (Wen et al., 2022; Wang et al., 2023).

Since MODIS and CASA both rely on remote sensing and meteorological drivers, they are not fully independent. Our dual validation provides convergent evidence but does not eliminate dependence. At the data level, the main analyses employ TerraClimate data resampled from 4 km to 1 km, a process that may attenuate fine-scale climatic signals. A scale check against higher-resolution NTPDC data at 1 km and 6 km scales showed close agreement with small deviations, indicating limited scale-related uncertainty for analysis. Cross-product validation at the 6 km scale confirmed high overall consistency. In contrast, a 1 km analysis revealed localized discrepancies, particularly within ecological transition zones. The anthropogenic footprint data lacked explicit disturbance types, which constrained the diagnosis of underlying mechanisms (Ellis, 2011). The 500 m resolution of MCD12Q1 limited the detection of fine-scale land use changes and exacerbated class mixing in urban-oasis transition zones. This led to overestimated NPP for construction land, causing local overestimation or underestimation and regional aggregation errors (Huang et al., 2021; Wu et al., 2023). Because wetland-water boundaries and seasonal inundation can blur at 500 m, local overestimation of wetland NPP is possible.

Future research should incorporate advanced moisture stress algorithms to improve accuracy in arid areas and calibrate the ε_max parameter locally using tower flux and plot-level data. To achieve a more accurate attribution of NPP dynamics, future work must incorporate high-resolution data on climate, land use, and explicit anthropogenic disturbance types, including grazing intensity and infrastructure development, to overcome scale-related biases. We can further integrate multi-source data, such as unmanned aerial vehicle (UAV) remote sensing and satellite-derived soil moisture, across scales to reconcile model discrepancies and reduce uncertainties, enabling more reliable assessments of ecosystem productivity in vulnerable transition zones.