Hydrothermal and Vegetation-Mediated Controls on Soil Organic Carbon in an Alpine Headwater Region of the Tibetan Plateau: Implications for Sustainable Grassland Management

1. Introduction

Soil organic carbon (SOC) is the largest terrestrial carbon reservoir and plays a pivotal role in regulating ecosystem resilience, supporting carbon neutrality goals, and contributing to climate change mitigation efforts [1,2]. In alpine grasslands, SOC stability is intricately linked to hydrological regulation, vegetation productivity, and permafrost dynamics. Even slight variations in soil temperature and moisture under ongoing climate warming can significantly alter the carbon balance, rendering alpine SOC both a vulnerable component of the climate system and a critical target for sustainable land management [3].The Tibetan Plateau, often referred to as the “Third Pole,” is home to vast alpine grasslands and permafrost systems that are highly sensitive to warming [4]. These grasslands store approximately 7.4 Pg of SOC within the upper 100 cm of soil, with significant spatial variation largely driven by moisture availability [5,6]. However, warming-induced permafrost degradation may accelerate organic matter decomposition and carbon release, while increases in vegetation productivity may simultaneously enhance carbon inputs [7,8]. These opposing processes generate complex, spatially heterogeneous SOC responses, adding uncertainty to the assessment of regional carbon sink stability.

At regional scales, SOC in alpine environments is characterized by strong spatial heterogeneity, typically exhibiting right-skewed distributions with localized high-carbon patches embedded within broader low-carbon areas [4]. These spatial patterns suggest that regional carbon stability is strongly dependent on hydrothermally favorable zones, including wetlands and seasonally moist grasslands. Identifying and protecting these high-SOC areas is thus essential for sustainable grassland management and climate adaptation strategies[9]. The dynamics of SOC in these alpine systems are governed by interacting thermal, hydrological, biological, and edaphic processes[10]. Temperature influences microbial activity and freeze–thaw disturbances, moisture availability governs vegetation growth and soil redox conditions, and soil texture and mineral composition affect carbon stabilization [11]. These tightly coupled processes generate nonlinear and hierarchical responses, which are difficult to capture using conventional linear models[12,13,14]. Recent advancements in remote sensing have provided spatially explicit data on vegetation productivity, surface temperature, moisture status, and cryospheric processes, offering valuable tools for regional carbon assessments in remote alpine regions [15,16,17]. While machine learning methods, particularly Random Forest (RF) models, have significantly improved SOC prediction accuracy, many applications focus on predictive performance without addressing the underlying processes or providing actionable management insights[18,19]. Bridging predictive modeling with process-based interpretation is therefore critical for translating spatial carbon assessments into effective land management strategies[20].

This study focuses on the Yellow River Source Region (YRSR), a climate-sensitive alpine headwater area located on the northeastern Tibetan Plateau. By integrating field-measured SOC data with multi-source environmental predictors, the objectives of this study are to: (1) characterize the spatial heterogeneity of SOC, (2) identify the dominant hydrothermal and biotic drivers of SOC variation, and (3) uncover the direct and indirect regulatory pathways influencing SOC distribution. From a sustainability perspective, understanding these mechanisms is crucial for identifying vulnerable carbon pools, prioritizing the conservation of high-moisture areas, and developing adaptive grassland management strategies under climate change.

2. Materials and Methods

2.1. Study Area

The Yellow River Source Region (YRSR), located on the northeastern Tibetan Plateau, is the headwaters and a critical water supply area for the Yellow River. As an ecologically important region, it plays a central role in regulating regional hydrological processes, maintaining alpine ecosystem stability, and supporting socio-economic development downstream[21,22]. Consequently, environmental changes in the YRSR have implications not only for local ecosystem functioning but also for water security and sustainable development across the broader Yellow River Basin. Originating from the Bayan Har Mountains, the study area spans elevations of approximately 2680–6248 m, with topography generally decreasing from northwest to southeast (Figure 1) [23,24,25]. The landscape is characterized by extensive glaciers, seasonal snow cover, lakes, and permafrost, representing a typical alpine ecosystem of the Tibetan Plateau [26,27]. The regional climate is dominated by a cold continental regime, with low annual mean temperatures, a short growing season, and frequent freeze–thaw cycles. Precipitation occurs mainly between May and September and exhibits pronounced seasonal and spatial heterogeneity [28,29]. Vegetation is primarily composed of alpine meadows, alpine steppes, and wetlands, with localized floodplains, marshes, and seasonally inundated depressions. Wetland and marsh areas, characterized by low temperatures and high moisture conditions, often develop peat-rich or humus-enriched soils with relatively high SOC accumulation. In contrast, surrounding alpine grasslands and wind-eroded uplands typically exhibit lower SOC contents [23,30,31]. Overall, SOC dynamics in the YRSR are jointly regulated by climatic gradients, hydrological processes, vegetation productivity, and soil–topographic controls [32,33,34]. Low temperatures constrain decomposition rates, while spatial variability in moisture availability mediates vegetation growth and soil hydrothermal conditions. Given the region’s sensitivity to climate warming and hydrological shifts, understanding SOC spatial patterns in the YRSR is essential for safeguarding carbon stability and supporting sustainable alpine land management [35,36,37].

2.2. Data Sources

2.2.1. Field Sampling and Laboratory Analysis

Soil sampling was conducted during July–August 2023 across the YRSR on the Qinghai-Tibet Plateau. A total of 240 sampling sites were selected to capture the spatial variability of SOC. To ensure representative coverage of environmental gradients with a limited number of samples, conditional Latin hypercube sampling (cLHS) was used to design the sampling sites. Auxiliary variables representing topography, climate, vegetation, and hydrological conditions were resampled to a 500 m resolution. Road networks, settlements, elevation, and slope were incorporated as constraints to construct a sampling cost surface, assigning penalty weights to inaccessible areas (Figure 2). The sampling design was implemented in R (version 4.5.2) [38]. Samples were collected using a handheld soil auger (5 cm diameter) from the 0–20 cm soil layer to represent surface SOC characteristics. At each sampling location, three subsamples were collected using a randomized design and combined into a composite sample to reduce the influence of small-scale spatial heterogeneity. All samples were air-dried, cleared of plant residues and gravel, ground, and passed through a 0.25 mm sieve prior to analysis. SOC content was determined using the potassium dichromate–sulfuric acid external heating oxidation method following standard procedures [39,40].

2.2.2. Construction of Environmental Predictor Set

To ensure that variable grouping and pathway specification were based on ecological processes rather than purely statistical associations, we adopted the STEP–AWBH structural framework to organize the environmental predictors[41]. This framework conceptualizes SOC formation and spatial variability as emergent outcomes of interacting subsystems operating across hierarchical levels within the regional Earth system. Specifically, environmental drivers were partitioned into eight interacting components: soil properties (S), topography (T), ecological and vegetation processes (E), precipitation (P), atmospheric conditions (A), hydrological processes (W), biotic activity (B), and human influences (H) [41]. This process-oriented framework enabled the systematic representation of thermal regulation, moisture mediation, biological productivity, and anthropogenic disturbance within a unified conceptual model. Following this framework, 37 initial candidate environmental variables were compiled from multi-source remote sensing products and reanalysis datasets (Table 1), ensuring that variable selection reflected ecological mechanisms relevant to SOC regulation rather than purely statistical associations.

2.2.3. Sentinel-2 Remote Sensing Data and Cloud Masking

High-resolution environmental variables representing vegetation productivity, surface moisture status, and cryospheric dynamics were derived from Sentinel-2 Multispectral Instrument (MSI) Level-2A surface reflectance products (10–20 m spatial resolution). The red-edge and near-infrared bands were used to compute vegetation indices sensitive to canopy structure and physiological activity, which are closely linked to carbon input processes in alpine ecosystems [42,43]. To reduce the influence of frequent cloud cover and terrain-induced shadows, all imagery was subjected to strict cloud and cloud-shadow masking prior to index calculation. Variable-specific temporal compositing strategies were applied to reflect ecological seasonality. Leaf Area Index (LAI) was derived from Sentinel-2 Level-2A surface reflectance data using the ESA SNAP Biophysical Processor, which applies a neural network trained on radiative transfer simulations to retrieve biophysical variables from red-edge and shortwave infrared bands. Growing-season (May–September) LAI was composited using the median to represent peak vegetation productivity. Surface moisture-related indices were calculated as growing-season means to represent average hydrothermal conditions. Snow-related metrics, including the normalized difference snow index (NDSI), were computed for the cold season (October–April), and snow frequency was quantified as the proportion of snow-covered observations relative to the total number of valid observations during the cold season [44,45,46]. This seasonally stratified preprocessing framework minimized atmospheric noise and ensured that derived variables were ecologically consistent with SOC formation and stabilization processes.

2.3. Methods

2.3.1. Screening and Optimization of Environmental Predictors

To construct an independent and ecologically interpretable predictor set, a two-stage variable convergence strategy combining relevance screening and multicollinearity control was implemented. First, the Boruta algorithm was applied to evaluate the relevance of each environmental predictor to SOC. Boruta compares the importance of original predictors with that of randomly permuted shadow variables. For each predictor${\mathrm{X}}_{\mathrm{i}}$, its importance scoreis${\mathrm{I}(\mathrm{X}}_{\mathrm{i}})$ compared with the maximum importance of shadow features ${\mathrm{I}}_{\mathrm{s}\mathrm{h}\mathrm{a}\mathrm{d}\mathrm{o}\mathrm{w}}^{\mathrm{m}\mathrm{a}\mathrm{x}}$. A predictor was confirmed as relevant when:

$${{I(X}_i)>I }_{shadow}^{max}$$

(1)

consistently across repeated iterations. This procedure retained only predictors with statistically significant explanatory power.

Second, multicollinearity among retained predictors was assessed using the Variance Inflation Factor (VIF) [47]. For each predictor X_i, VIF is defined as:

$${VIF}_i ={\displaystyle \frac{1}{1-R_i^2}}$$

(2)

where ${R }_i^2$is the coefficient of determination obtained by regressing predictor Xᵢ against all remaining predictors. A high ${R }_i^{2 }$indicates strong linear dependence on other variables. Predictors with $VIF\ge 5$were iteratively removed, beginning with the variable exhibiting the highest VIF value, while maintaining ecological representativeness [48]. The final predictor set satisfied:

$${VIF}_i <5 \forall i$$

(3)

ensuring statistical independence and model stability for subsequent analyses.

2.3.2. Random Forest Modeling and SHAP Interpretation

Random forest (RF) models were constructed using the optimized predictor set to characterize nonlinear relationships between SOC and environmental drivers. RF is an ensemble learning algorithm that aggregates predictions from multiple regression trees built on bootstrap samples of the training data. The ensemble prediction for an observation x is expressed as:

$$\widehat{y_i}\left(x\right)={\displaystyle \frac{1}{T}} \sum _{t=1}^T{f}_t\left(x\right)$$

(4)

where T denotes the number of trees and $f_t\left(x\right)$ represents the prediction from the t-th tree.

Model performance was evaluated using the coefficient of determination ($R^2$) and root mean square error (RMSE) [49]. The coefficient of determination was calculated as:

$${\mathrm{R}}^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-\widehat{y_i})}^2}{{\sum }_{i=1}^n{(y_i-\overline{y_i})}^2}}$$

(5)

where $y_i$denotes observed SOC, ${\widehat{y}}_i$represents predicted SOC, and $\bar{y}$is the mean of observed SOC values. RMSE was defined as:

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n{\left(y_i-\widehat{y_i}\right)}^2}$$

(6)

To enhance interpretability, Shapley Additive exPlanations (SHAP) were applied to decompose RF predictions into additive contributions from individual predictors. For fea- ture $i$, the SHAP value ${\varphi }_i$quantifies its marginal contribution to the prediction and is de fined as:

$${\varphi }_i=\sum _{S\subseteq N\backslash \left\{i\right\}}{\displaystyle \frac{\left|s\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}}\left[f\left(S\cup \left\{i\right\}\right)-f\left(S\right)\right]$$

(7)

where $N$represents the set of all predictors and $S$denotes a subset not containing feature $i$.

Global variable importance was assessed using the mean absolute SHAP value:

$$\left|\overline{{\mathrm{\varnothing }}_i}\right|={\displaystyle \frac{1}{n}}\sum _{j=1}^n\left|{\mathrm{\varnothing }}_{ij}\right|$$

(8)

which reflects the average contribution magnitude of predictor $i$ across all observations. SHAP dependence patterns were examined to characterize nonlinear SOC responses to dominant environmental drivers.

SHAP attribution was used strictly for model interpretation and identification of influential predictors, rather than for causal inference.

2.3.3. Spatial Prediction and Distribution Analysis of SOC

The spatial distribution of SOC was derived from RF predictive model. The model-generated SOC raster was resampled to a 500 m spatial resolution and clipped to the boundary of the Yellow River Source Region. A multi-year mean SOC layer was subsequently calculated to represent the long-term spatial pattern. To evaluate differences among vegetation types, the predicted SOC raster was overlaid with land cover data. Zonal statistical analysis was performed to extract descriptive statistics [8,9]. This approach provided a spatially explicit characterization of SOC heterogeneity across grassland types and supported subsequent interpretation of environmental controls [50,51].

3. Results

3.1. Characteristics of Soil Organic Carbon Measurements

SOC at a depth of 0–20 cm in the Yellow River Source Region (YRSR) exhibited considerable variability, ranging from 2.96 to 169.32 g/kg (Table 2). The mean and median values were 46.61 g/kg and 40.75 g/kg, respectively, indicating moderate central tendency but substantial dispersion. The coefficient of variation (CV) was 64.97%, reflecting high relative variability across sampling sites. The SOC distribution was positively skewed (skewness = 1.22) and leptokurtic (kurtosis = 4.85), suggesting a right-skewed and heavy-tailed pattern. Most observations clustered in lower SOC ranges, while a few high-value samples extended the upper tail of the distribution. This distribution pattern indicates that localized high-SOC patches disproportionately contribute to regional variability. Spatially, elevated SOC concentrations were concentrated in discrete areas, whereas lower SOC values dominated the broader regional background. This patch-matrix configuration implies that SOC accumulation is strongly regulated by localized hydrothermal and ecological conditions, rather than by uniform regional controls. The observed non-normality and spatial concentration patterns underscore the importance of employing spatially explicit and nonlinear approaches to resolving SOC variability in alpine environments.

3.2. Identification of Key Environmental Predictors for SOC

The Boruta algorithm was applied to assess the relevance of 37 candidate environmental variables to soil organic carbon (SOC) in the Yellow River Source Region (YRSR). After iterative comparison with randomly permuted shadow variables, 27 predictors were identified as significant contributors to SOC variability, with their importance values significantly exceeding those of the random shadow variables (Figure 3). These retained variables encompassed key environmental factors, including thermal conditions (e.g., land surface temperature), moisture processes (e.g., precipitation and soil moisture), vegetation productivity (e.g., Leaf Area Index (LAI), Gross Primary Productivity (GPP)), and soil properties (e.g., pH, bulk density). To control for multicollinearity among the retained predictors, the Variance Inflation Factor (VIF) was applied. Predictors with VIF values greater than 5 were iteratively removed until all remaining variables exhibited acceptable VIF values (<5), ensuring statistical independence and reducing potential model instability (Table 3). The final predictor set included 16 variables, which adequately represent the dominant environmental gradients shaping SOC spatial variability in the region. These variables will serve as robust inputs for subsequent nonlinear modeling, offering a comprehensive, ecologically interpretable basis for assessing SOC in the YRSR.

3.3. Nonlinear Importance Patterns of SOC Revealed by RF–SHAP

The random forest (RF) model demonstrated robust predictive performance, explaining 68% of the variance in soil organic carbon (SOC), with an RMSE of 15.59 g/kg (Figure S1). These results indicate that the optimized environmental predictors effectively captured the primary spatial gradients in the Yellow River Source Region (YRSR). Among the predictors, vegetation productivity indicators such as LAI and GPP were the most significant contributors, with high SHAP values indicating their strong role in SOC accumulation. These variables demonstrated a nonlinear relationship with SOC, where SOC increased with higher productivity but plateaued at certain thresholds. Moisture-related variables (such as NDWI and SM) were also identified as crucial drivers, with SHAP highlighting their nonlinear contributions to SOC variation. Moderate to high moisture conditions were associated with higher SOC due to enhanced vegetation growth and reduced decomposition. In contrast, drier conditions corresponded to lower SOC values, underscoring the importance of maintaining a stable hydrothermal balance. Thermal variables, particularly daytime land surface temperature (Terra_LsTD) and the AI, showed significant contributions to SOC variability, with SHAP results indicating that lower temperatures were associated with higher SOC levels due to suppressed microbial activity. However, extreme thermal fluctuations, especially from freeze–thaw cycles, were found to reduce SOC stability, signaling the complex role of temperature variability in SOC dynamics (Figure 4). Overall, the results from RF–SHAP emphasize the nonlinear and context-dependent nature of SOC regulation in the YRSR, driven by the coupled effects of moisture, temperature, and vegetation productivity. These findings provide a strong empirical foundation for the subsequent modeling of SOC dynamics, which can inform the management of alpine ecosystems.

3.4. Spatial Heterogeneity and Vegetation-Type Differences in SOC in the YRSR

SOC in the Yellow River Source Region exhibited pronounced spatial heterogeneity in 2023(Figure 5). The multi-year mean SOC during the study period was 69.2 g/kg. Overall, SOC showed a clear spatial gradient, decreasing from southeast to northwest. High SOC values were primarily concentrated in the southeastern part of the region, whereas relatively low SOC levels were observed in the northwestern areas. Significant differences in SOC were identified among vegetation types. Mountain meadow exhibited the highest mean SOC (85.63 g/kg), representing the most carbon-rich ecosystem in the region. Forest land, improved grassland, lowland meadow, and wetlands also showed relatively high SOC levels, with mean values exceeding 76 g/kg. Alpine meadow and alpine steppe had intermediate SOC levels, with mean values of 70.63 g/kg and 53.12 g/kg, respectively. Temperate steppe showed a mean SOC of 58.75 g/kg, while temperate desert grassland exhibited the lowest SOC content (42.40 g/kg). In terms of variability, alpine meadow and wetland displayed relatively large SOC ranges, with standard deviations of 17.72 g/kg and 16.82 g/kg, respectively, indicating stronger spatial fluctuations. In contrast, improved grassland and lowland meadow showed smaller standard deviations, suggesting a more uniform spatial distribution of SOC. Overall, SOC in the Yellow River Source Region demonstrated clear vegetation-type differentiation and a distinct spatial gradient, reflecting substantial heterogeneity in carbon storage among different grassland ecosystems.

4. Discussion

4.1. Spatial Heterogeneity of Soil Organic Carbon in the Yellow River Source Region

The surface OC in the Yellow River Source Region (YRSR) exhibited significant spatial heterogeneity, which is consistent with the findings of previous studies in other alpine regions. The observed variation in SOC, ranging from 2.96 to 169.32 g/kg, indicates that SOC in the region is highly variable, with localized high-carbon patches embedded within broader low-carbon backgrounds. The mean value of 46.61 g/kg and the high coefficient of variation (CV = 64.97%) reflect a high degree of spatial variability in SOC concentrations. This variability is likely driven by the complex interactions of multiple environmental factors such as moisture, temperature, and vegetation. The positive skewness (skewness = 1.22) and leptokurtic distribution (kurtosis = 4.85) suggest a right-skewed and heavy-tailed distribution, where most of the SOC values are clustered in the lower range, while a few high-value samples extend the distribution’s upper tail. This distribution pattern implies that localized high-SOC patches play a disproportionate role in contributing to the regional SOC variability [52,53]. The spatial clustering of high SOC concentrations in specific areas indicates that localized hydrothermal and ecological conditions, rather than uniform regional controls, govern SOC accumulation in the YRSR [34,54]. These findings highlight the importance of using spatially explicit and nonlinear approaches to accurately capture SOC variability in heterogeneous alpine [55].

4.2. Key Environmental Drivers of SOC Variability Identified by RF–SHAP

The RF–SHAP analysis revealed that vegetation productivity, moisture availability, and thermal conditions are the primary contributors to SOC variability in the YRSR. These findings support existing literature, which highlights the strong regulatory effects of hydrothermal conditions and vegetation productivity on SOC dynamics in alpine ecosystems. Vegetation productivity indicators, such as the LAI and GPP, were found to be significant drivers of SOC accumulation. Higher vegetation productivity generally promotes SOC storage due to increased carbon inputs from aboveground biomass [56]. Moisture-related variables, including precipitation and soil moisture, were identified as critical drivers in SOC variability, displaying strong nonlinear contributions. SOC levels were higher under moderate to high moisture conditions, corresponding to enhanced vegetation growth and reduced decomposition in moist soils, while drier conditions led to lower SOC levels. These findings are consistent with previous studies that have shown moisture availability as a key regulator of SOC in cold alpine environments[57,58]. Thermal conditions, particularly freeze-thaw dynamics and surface temperature, were also significant contributors to SOC variability. These results align with studies that have shown that temperature strongly influences microbial decomposition and carbon turnover[59,60]. Lower temperatures were generally associated with higher SOC concentrations, reflecting suppressed microbial activity, while extreme thermal fluctuations, particularly freeze-thaw cycle s, reduced SOC stability [61,62,63].

4.3. The Role of Nonlinear Interactions in SOC Regulation

A key finding of this study is the importance of nonlinear interactions between environmental variables in shaping SOC patterns. The RF–SHAP analysis demonstrated that the relationship between SOC and its environmental drivers, such as moisture and temperature, is not linear [64,65,66,67]. For instance, higher vegetation productivity led to increased SOC accumulation up to a certain threshold, after which the relationship plateaued. This indicates that carbon stabilization processes may constrain SOC accumulation beyond certain productivity levels [68,69]. Moisture availability showed similar nonlinear behavior, with moderate to high moisture conditions promoting SOC accumulation, while drier conditions resulted in lower SOC levels. These findings suggest that SOC regulation is highly context-dependent and nonlinear, emphasizing the need for advanced models that can account for these complex interactions [70,71,72].

4.4. Context-Dependent Effects of Vegetation Productivity and Biological Inputs

The results of this study have important implications for the sustainable management of grasslands in the YRSR. Given the strong influence of moisture, vegetation productivity, and thermal conditions on SOC dynamics, it is crucial to prioritize the conservation of areas with high SOC potential, particularly those in hydrothermally favorable zones such as wetlands and seasonally moist grasslands[73,74]. These areas, which serve as significant carbon sinks, should be prioritized for conservation efforts to maintain regional SOC stocks. The findings also highlight the need for adaptive management strategies that consider the effects of climate warming on SOC dynamics [75]. As temperatures rise and precipitation patterns shift, the stability of SOC in the region may be threatened [76]. Therefore, management strategies should focus on enhancing vegetation productivity in moisture-limited areas and maintaining stable thermal conditions to protect SOC stocks [76,77]. These strategies are critical to ensuring the long-term sustainability of alpine grasslands under future climate scenarios [78].

4.5. Implications for Remote Sensing–Based SOC Modeling in Alpine Regions

Our study emphasizes the need to combine diverse remote sensing data sources to develop comprehensive SOC models that better capture interactions between these variables. diverse remote sensing data sources to build comprehensive SOC models that better reflect the interactions between these variables. Thermal and moisture-related remote sensing indicators, ncluding LST, freeze–thaw cycles, and NDWI, were identified as critical predictors of SOC variability in the YRSR. These variables provide a spatially explicit representation of hydrothermal conditions that influence microbial activity, carbon stabilization, and vegetation productivity. High-resolution moisture proxies, such as NDWI, are particularly effective in capturing local hydrological variability and SOC hotspots, complementing coarser-resolution data from reanalysis and microwave products[48]. The integration of such high-resolution moisture proxies with thermal and vegetation productivity data enhances the spatial reliability and ecological understanding of SOC distribution in alpine ecosystems [74]. One of the key advantages of using machine learning techniques, such as random forest (RF), in combination with remote sensing data is their ability to capture nonlinear relationships between SOC and environmental predictors. The RF–SHAP analysis demonstrated that moisture availability, temperature, and vegetation productivity interact in complex, nonlinear ways to regulate SOC [49]. These findings highlight the importance of using advanced modeling techniques that can account for the intricate dynamics of SOC accumulation, particularly in regions with strong environmental gradient Looking ahead, future research could focus on expanding the integration of remote sensing data with machine learning models to further improve SOC predictions in regions with limited field data. In particular, the use of high-resolution, multi-sensor data (e.g., optical, thermal, and radar remote sensing) could provide more comprehensive insights into SOC dynamics, especially in remote and data-scarce regions like the Tibetan Plateau [79,80].

5. Conclusions

This study quantified the spatial heterogeneity of surface soil organic carbon (SOC) in the Yellow River Source Region, identifying the dominant hydrothermal and biotic controls using an integrated RF–SHAP framework. The results revealed significant spatial variability in SOC, with moisture-related processes acting as a central mediator between thermal conditions, vegetation productivity, and soil properties. Thermal regimes and freeze-thaw dynamics directly influenced SOC, while vegetation productivity and soil properties primarily affected SOC through moisture-mediated pathways. These findings underscore the hierarchical and nonlinear nature of SOC regulation in alpine headwater ecosystems. From a sustainability perspective, maintaining favorable hydrothermal conditions is crucial for stabilizing soil carbon stocks. Prioritizing wetland conservation, adaptive grazing management, and reducing soil disturbance in high-SOC areas will help mitigate carbon loss, especially under climate warming. The integrated modeling approach developed in this study provides a transferable method for identifying vulnerable carbon pools and supporting sustainable land management in other alpine and climate-sensitive regions.