Spatial and temporal patterns of vegetation carbon sources/sinks and driving factors in southeastern Xizang from 2000 to 2020

Data quality control and SLC-OFF gap handling

To address the SLC-OFF striping gaps and cloud/shadow interference in Landsat-7 imagery after 2003, a rigorous data quality control procedure was applied. Pixels affected by clouds, cloud shadows, snow/ice, and saturation were first removed using the QA flags from the Collection 2 Surface Reflectance (SR) dataset, with a cloud probability threshold of 20%. Within the growing season (June–September), multiple images acquired in the same month were screened for quality to derive monthly NDVI composites, and an annual medoid composite was then generated to represent each year, allowing the temporal synthesis process to naturally fill the striping gaps. Pixels with fewer than three valid observations were flagged for subsequent filtering. To reduce systematic biases caused by observation geometry and terrain effects, Bidirectional Reflectance Distribution Function (BRDF) normalization was performed using MCD43 parameters, combined with terrain correction based on the SCS+C model and SRTM DEM.

Spatial and temporal alignment of multi-source data

A unified spatial resolution of 30 m × 30 m was adopted as the target grid. The spatial resolutions of Landsat-7 NDVI and GLC_FCS30 (30 m) were already consistent with the target grid, and land-cover data were mapped by category without interpolation. TerraClimate variables, including monthly mean temperature, total precipitation, and total solar radiation (4 km), were resampled to 30 m using bilinear interpolation and co-registered with NDVI and land-cover grids to drive the CASA model on a monthly scale. MOD17A3H annual NPP (500 m) was used as an external reference and resampled to 30 m for pixel-to-pixel comparison. Temporally, all datasets were aligned on a monthly basis: monthly NDVI was used to estimate the fraction of photosynthetically active radiation (FPAR) and ɛ, which, together with monthly meteorological data, were used to calculate monthly NPP; soil heterotrophic respiration was also simulated monthly and subsequently accumulated to obtain annual NEP.

Discretization of continuous factors and robustness testing

For the analysis of driving factors, continuous variables such as precipitation, temperature, elevation, and slope were discretized into five classes using the natural breaks (Jenks) method to balance within-group variance minimization and sample size uniformity, while aspect, as a categorical variable, was used directly. To assess the robustness of results, both equal-interval and equal-frequency classifications were applied for comparison. Although the absolute q-statistics differed slightly across classification schemes, the relative ranking and interaction types of dominant factors remained consistent, indicating robust conclusions.

NPP estimation

The CASA model, proposed by Potter et al. (1993), is a typical large-area-scale light energy use model driven by a combination of remotely sensed data, meteorological data, and land cover data. The NPP depends on how much photosynthetically active radiation (APAR) it absorbs relative to its light energy conversion efficiency (ɛ). Zhu, Pan & Zhang (2007) enhanced the existing CASA model by incorporating vegetation cover classification into the model. The NPP formula for the improved CASA model is as follows: (1)${\displaystyle {\mathrm{NPP}}_{\mathrm{(x,t)}}={\mathrm{APAR}}_{\mathrm{(x,t)}}\times {\varepsilon }_{\mathrm{(x,t)}}}$ where NPP is the net primary productivity of vegetation, $\mathrm{APA}{\mathrm{R}}_{\left(\mathrm{x},\mathrm{t}\right)}$ is the photosynthetically active radiation absorbed by image x in month t (MJ m⁻²), and ${ϵ}_{\left(\mathrm{x},\mathrm{t}\right)}$ is the efficiency of light energy utilization by image x in time period t.

The absorbed photosynthetically active radiation (APAR) in a specific area is determined by the proportion of total incoming radiation that the vegetation absorbs, calculated using the following formula: (2)${\displaystyle {\mathrm{APAR}}_{\mathrm{(x,t)}}={\mathrm{SOL}}_{\mathrm{(x,t)}}\times {\mathrm{FPAR}}_{\mathrm{(x,t)}}\times 0.5}$ where SOL_(x,t) is the total solar radiation (MJ m²) received within pixel x in the t-th time period, FPAR_(x,t) is the fraction of photosynthetically active radiation absorbed by the vegetation within pixel x in the t-th time period, and 0.5 denotes the proportion of total solar radiation used for photosynthesis.

The total solar radiation was computed using an established formula developed by He and Xie (He and Xie., 2010), which is specific to the development of western China: (3)${\displaystyle \mathrm{Q}={\mathrm{Q}}_0\times \left(\mathrm{a}+\mathrm{b}\times {\mathrm{T}}_{\mathrm{S}}/{\mathrm{T}}_{\mathrm{A}}\right)}$ where Q₀ is astronomical radiation, which is calculated on the basis of geographic latitude, solar declination and other parameters; T_S is observed sunshine hours; T_A is theoretical sunshine hours, which is calculated on the basis of latitude and solar declination; and a and b are empirical coefficients (a = 0.185; b = 0.595).

The formula for calculating the fraction of photosynthetically active radiation absorbed (FPAR) by vegetation is as follows: (4)${\displaystyle {\mathrm{FPAR}}_{\mathrm{(x,t)}}=\alpha {\mathrm{FPAR}}_{\mathrm{NDVI}}+\left(1-\alpha \right){\mathrm{FPAR}}_{\mathrm{SR}}}$ (5)${\displaystyle FPA{R}_{NDVI}=\frac{NDV{I}_{\left(x,t\right)}-NDV{I}_{\left(i,min\right)}}{NDV{I}_{\left(i,max\right)}-NDV{I}_{\left(i,min\right)}}\times \left(FPA{R}_{max}-FPA{R}_{min}\right)+FPA{R}_{min}}$ (6)${\displaystyle FPA{R}_{SR}=\frac{S{R}_{\left(x,t\right)}-S{R}_{\left(i,min\right)}}{S{R}_{\left(i,max\right)}-S{R}_{\left(i,min\right)}}\times \left(FPA{R}_{max}-FPA{R}_{min}\right)+FPA{R}_{min}}$ (7)${\displaystyle S{R}_{\left(x,t\right)}=\frac{1+NDV{I}_{\left(x,t\right)}}{1-NDV{I}_{\left(x,t\right)}}}$

where t is the month and the rest of the parameters are as above. FPAR_max and FPAR_min are the maximum and minimum photosynthetically active radiation ratios taken as 0.95 and 0.001, respectively. SR_(i,max) and SR_(i,min) are the maximum and minimum values of the simple ratio of the normalized index for the ith vegetation type, obtained in a probability distribution of 95% and 5%, respectively. NDVI_(i,max) and NDVI_(i,min) are the maximum and minimum values of the NDVI for the ith vegetation type, taking the DN values corresponding to a probability distribution of 95% and 5%. SR_(x,t) is the simple ratio of the normalized vegetation index within pixel x in time period t.

The actual light use efficiency (ɛ) is influenced primarily by temperature, water availability, and vegetation type. The effective light use efficiency is determined by multiplying the vegetation’s maximum light use efficiency by factors related to water and temperature stress. The calculation formula is as follows: (8)${\displaystyle {\varepsilon }_{\mathrm{(x,t)}}=T_{\varepsilon 1\left(\mathrm{x,t}\right)}\times T_{\varepsilon 2\left(\mathrm{x,t}\right)}\times W_{\varepsilon \left(\mathrm{x,t)}\right)}\times {\varepsilon }_{\max }}$ where T_ɛ1 and T_ɛ2 are the stress effect coefficients of the maximum and minimum temperatures on the actual light use efficiency ɛ_(x,t), W_ɛ(x,t) is the water stress effect coefficient, and ɛ_max is the maximum light energy utilization under ideal conditions (gC MJ⁻¹).

Table 1 lists the NDVI_(i,min), NDVI_(i,max), NDVI_(i,min), NDVI_(i,max) and ɛ_max values of the different vegetation types (Zhu, Pan & Zhang, 2007).

NEP estimation

NEP serves as a critical indicator of regional vegetation carbon dynamics, representing the difference between NPP and carbon lost through soil microbial respiration, excluding other natural and human-induced influences, which is calculated as follows: (9)${\displaystyle \mathrm{NEP(x,t)}=\mathrm{NPP(x,t)}-\mathrm{RH(x,t)}}$ where NEP(x, t) represents the net ecosystem productivity of vegetation for pixel x during time period t (gC m⁻²), NPP(x, t) represents the net primary productivity of vegetation for pixel x during time period t (gC m⁻²), and RH(x, t) represents the soil microbial respiration of pixel x during time period t (gC m⁻²).

The calculation formula for soil heterotrophic respiration is as follows (Zhiyong et al., 2003): (10)${\displaystyle {\mathrm{R}}_{\mathrm{H}}=3.069\times \left({\mathrm{e}}^{0.0912\mathrm{T}}+\mathrm{In}\left(0.3145\times \mathrm{R}+1\right)\right)}$ where R_H represents soil heterotrophic respiration (gC m⁻²), T represents the temperature (°C), and R represents precipitation (mm).

Theil–Sen trend analysis and testing methods

The Theil–Sen median (Sen) method is less affected by outliers, making it more robust. (Lu et al., 2023). It is a widely used metric for analyzing trends in long-term data series. It measures the direction and magnitude of trends in NEP: a positive β signifies an upward trend, while a negative β signifies a downward trend. The formula for calculating β is as follows: (11)${\displaystyle \mathrm{\beta }=\mathrm{Median}\left({\mathrm{NEP}}_{\mathrm{j}}-{\mathrm{NEP}}_{\mathrm{i}}\right)/\left(\mathrm{j}-\mathrm{i}\right)2000\le i<j\le 2020}$ where b is Sen’s slope value, and NEP_i andNEP_j are the NEPs in years i and j, respectively. The Mann−Kendall test was used to test the significance of the results of the β-trend analyses.

GeoDetector

The GeoDetector model is a tool for detecting spatial heterogeneity and is capable of revealing the drivers behind spatial differentiation (Zhang et al., 2023b); it includes four modules: the factor detector, interaction detector, risk detector, and ecological detector. This study employed the GeoDetector model to investigate the driving factors of spatial variation in the NEP of vegetation in southeastern Xizang.

Factor detection: this tool assesses the extent to which the driving factor X (Average multiyear total precipitation (2000–2020, mm); Average multiyear average temperature (2000–2020, ^∘C); Digital Elevation Model (DEM, m); Slope (^∘); Aspect (^∘) (°C)) influences Y (NEP), employing the q value metric to quantify this impact. which is calculated as follows: (12)${\displaystyle q=1-\frac{1}{f{\sigma }^2}\sum _{h=1}^{L}fk{\sigma }^{2}k}$

where L is the classification of factor X or variable Y; f and fk represent the sample sizes in the entire study area and within type h, respectively; and σ² and σ²k represent the discrete variances of Y values across the entire study area and within type h, respectively. The value range of q is [0, 1]. A higher q value signifies that the driver factor X has a stronger explanatory power on the NEP, while a lower q value indicates a weaker explanatory power.

Interaction detector (GeoDetector software 1.6, developed by Professor Jiumiao Cheng’s team at Renmin University of China; freely available at http://geodetector.cn/index.html): a statistical tool used to evaluate the impact of multiple factors on a specific variable. By detecting the interactions between different factors, the interaction detector can reveal how these factors work together and whether their combined effects are synergistic, antagonistic, or independent. Details of the interaction types are given in Table 2.