Machine learning-based modeling of soil respiration components across four forest types in a subtropical mountainous forest with implications for modeling strategies and hydrothermal controls

Soil respiration (Rs) is a key component of the terrestrial carbon cycle, yet its accurate prediction remains challenging due to complex interactions among environmental drivers and the distinct mechanisms governing its components. In this study, we conducted monthly 24-h measurement campaigns at hourly intervals of total soil respiration (Rs) and its components, autotrophic respiration (Ra) and heterotrophic respiration (Rh), across four typical forest types in the Wuling Mountains, a representative subtropical mountainous region in central–southern China. Multiple machine learning (ML) models were developed under different modeling strategies, feature engineering schemes, and validation protocols. Component-based modeling markedly outperformed unified and forest-type-based approaches, suggesting that under the current experimental design, differences among respiration components exerted a stronger influence on predictive performance than differences among forest types. Under random validation, ensemble tree models achieved high accuracy for Rs and Rh (R 2 > 0.90), but their advantage narrowed substantially under leave-one-month-out blocked validation. Notably, CatBoost trained without temporal variables showed the most consistent temporal generalization and was therefore selected for ecological interpretation. SHAP analysis revealed that deepest monitored soil hydrothermal coupling and vertical thermal structure were the dominant controls on Rs and Rh, particularly the interaction between the deepest monitored soil temperature and moisture; depth-ablation analysis further confirmed the independent predictive contribution of deep-layer measurements. In contrast, Ra was primarily related to deep-layer daytime temperature and moisture balance but remained the least predictable component, reflecting unmeasured plant physiological controls. Our results highlight the importance of integrating appropriate modeling strategies, validation protocols, and machine learning techniques to improve prediction of soil respiration components. The findings also emphasize the role of hydrothermal coupling and vertical heterogeneity in regulating soil carbon fluxes in subtropical forest ecosystems.

Soil respiration (Rs) is one of the largest carbon fluxes from terrestrial ecosystems to the atmosphere and a fundamental component of the global carbon cycle (Raich et al., 2002). As a major pathway of CO 2 release from soils, its annual flux is far greater than that of anthropogenic fossil-fuel emissions, meaning that even small changes in Rs may have important consequences for climate-carbon feedbacks (Raich et al., 2002; Moyes and Bowling, 2016). In forest ecosystems, which represent the largest terrestrial carbon pool and an important long-term carbon sink, soil respiration forms a critical link between aboveground productivity and belowground carbon turnover (Smith, 2008; Bond-Lamberty et al., 2018). Improving the prediction of soil respiration is therefore essential for understanding ecosystem carbon dynamics and for reducing uncertainty in regional carbon budget assessments under climate change.

The dynamics of soil respiration are jointly regulated by multiple abiotic and biotic factors, including soil temperature, soil moisture, substrate supply, vegetation characteristics, and microbial activity (Klimek et al., 2020). Among these, soil temperature and moisture are widely regarded as the dominant environmental controls and can explain a substantial proportion of the temporal variation in Rs (Moyano et al., 2013; Ma et al., 2020; Fekadu et al., 2023). However, soil respiration responses are rarely governed by single factors alone. Instead, they are shaped by nonlinear interactions among hydrothermal conditions, vegetation, and soil properties across temporal and spatial scales (Bond-Lamberty and Thomson, 2010; Hereș et al., 2021). This complexity is particularly pronounced in subtropical mountainous forests, where strong environmental heterogeneity associated with topography and seasonal climate may generate highly variable soil carbon fluxes. As a result, accurately predicting soil respiration in such ecosystems remains a major challenge. Total soil respiration consists of two functionally distinct components: autotrophic respiration (Ra), derived mainly from roots and rhizosphere processes, and heterotrophic respiration (Rh), produced primarily by microbial decomposition of soil organic matter (Rey et al., 2002). These components differ substantially in their regulatory mechanisms. Rh is generally more directly constrained by environmental conditions and substrate decomposition, whereas Ra is additionally influenced by plant physiological activity, recent photosynthate supply, and carbon allocation dynamics (Kuzyakov and Gavrichkova, 2010; Tang et al., 2019). Consequently, Ra is often more difficult to predict and remains an important source of uncertainty in terrestrial carbon-cycle studies. Distinguishing Rs into Ra and Rh is therefore necessary not only for improving predictive accuracy, but also for clarifying the mechanistic basis of ecosystem carbon efflux.

Existing methods for estimating soil respiration can generally be divided into process-based models and empirical statistical models. Process-based models aim to represent ecosystem mechanisms explicitly and have been widely used to simulate soil carbon fluxes at multiple scales (Lawrence et al., 2019; Lienert and Joos, 2018). However, their application is often limited by structural assumptions, parameter uncertainty, and trade-offs between mechanistic detail and computational efficiency (Yuan et al., 2011; Jägermeyr et al., 2014). In contrast, empirical statistical models are easier to implement, but they typically assume relatively simple functional relationships and therefore often fail to capture the strong nonlinearity and interaction effects that characterize soil respiration processes. These limitations have motivated increasing interest in machine learning (ML) approaches. ML provides a flexible framework for modeling complex ecological processes because it can learn nonlinear relationships directly from data and integrate high-dimensional predictors without requiring explicit assumptions about functional form (Breiman, 2001; Friedman, 2001; Chen and Guestrin, 2016). In recent years, ML methods have been increasingly applied to soil carbon flux estimation, especially at global and regional scales and in agricultural ecosystems (Liu et al., 2023; Liu et al., 2024; Chen et al., 2026). In comparison, their application to forest soil respiration remains relatively limited, particularly for partitioned respiration components (Canteral et al., 2023; Vicentini et al., 2023; Ferdous et al., 2025). Moreover, most existing studies have focused primarily on total soil respiration (Rs) (Lu et al., 2021), while relatively few have explicitly evaluated whether different modeling strategies-such as unified modeling, vegetation-based grouping, or component-based modeling—affect predictive performance and ecological inference. For heterogeneous ecological datasets, model structure and grouping strategy may be as important as algorithm choice, yet this issue remains insufficiently addressed in current soil respiration research. Despite their strong predictive performance, ML models are often criticized for limited interpretability, which constrains their ecological usefulness beyond prediction (Lipton, 2018). In ecological studies, identifying dominant environmental drivers and understanding their interactions are as important as achieving high predictive accuracy. Recent advances in explainable artificial intelligence, particularly SHAP (SHapley Additive exPlanations), provide an effective way to interpret complex models by quantifying the contribution of individual predictors and revealing nonlinear effects and interactions (Lundberg and Lee, 2017). These methods make it possible to combine the predictive power of ML with ecologically meaningful interpretation, thereby offering new opportunities to investigate the mechanisms underlying soil respiration dynamics.

In this study, we investigated soil respiration in the Wuling Mountains of central-southern China, a representative subtropical mountainous forest region. Four typical forest types were selected: coniferous forest, coniferous–broadleaf mixed forest, broadleaf forest, and economic forest. Using the trenching method to partition soil respiration into autotrophic and heterotrophic components, we conducted continuous observations of total soil respiration and its components together with multilayer monitoring of soil environmental variables. Based on these data, we developed multiple machine learning models under different modeling strategies and feature schemes. The objectives of this study were to: (1) compare the performance of different modeling strategies and machine learning algorithms for predicting soil respiration and its components, and identify an optimal modeling framework; and (2) clarify the nonlinear effects of soil hydrothermal conditions and their vertical structure on different respiration components, and identify their dominant environmental controls using interpretable machine-learning methods.

Study area and experimental design

This study was conducted at the Hunan Cili Forest Ecosystem Research Station, Hunan Province, China (29°25′27″ N, 111°12′42″ E), a typical subtropical mountainous region with a monsoon climate characterized by hot, humid summers and mild winters. Four representative forest stands, all approximately 30 years old, were selected: (1) a pure coniferous forest dominated by Pinus massoniana; (2) a coniferous–broadleaf mixed forest dominated by Pinus massoniana, Liquidambar formosana, and Cinnamomum camphora; (3) an economic forest represented by a Citrus spp. orchard; and (4) a secondary broadleaf forest dominated by Quercus spp. In China, economic forests are defined as woodlands cultivated primarily for non-timber forest products including fruits, nuts, plant oils and medicinal herbs, instead of timber. A Citrus spp. orchard was used to represent this forest type in our study, because citrus orchards are the most widely distributed economic forest type in the study region. Within each forest stand, two adjacent subplots (3 m × 2 m), separated by approximately 2 m, were established to quantify total soil respiration and heterotrophic respiration, respectively. One subplot was left undisturbed for the measurement of total soil respiration (Rs), whereas the other was subjected to root exclusion using the trenching method following Hanson et al. (2000) to estimate heterotrophic respiration (Rh). The root-exclusion subplot was isolated on all four sides by trenches extending to the bedrock (approximately 1.5 m deep, depending on local soil depth), thereby severing both lateral surface roots and deeper roots. The trenches were lined with corrosion-resistant metal mesh (1.9 mm aperture) to prevent root re-ingrowth.

To maintain the root-exclusion treatment, newly emerged vegetation within the trenched subplot was removed weekly throughout the study period. The litter layer was retained to preserve natural surface substrate inputs, and root residues were not intentionally removed from the soil profile in order to minimize additional disturbance. Autotrophic respiration (Ra) was calculated as the difference between total and heterotrophic respiration (Ra = Rs − Rh). To reduce short-term disturbance effects caused by trenching, a two-month stabilization period was allowed before formal monitoring began. In each subplot, three PVC collars were permanently installed in a triangular arrangement, approximately 1 m apart from each other and at least 0.5 m from the subplot boundary. The three collars served as technical sampling points to characterize within-subplot variability, and measurements from the three collars were averaged to produce one subplot-level value per sampling time. The subplot was treated as the experimental unit in all analyses. Soil respiration monitoring was conducted from June 2024 to July 2025.

Soil respiration measurements and environmental monitoring

Soil CO 2 efflux was measured using a LI-870 portable soil CO 2 flux system (LI-COR Biosciences, Lincoln, NE, USA). Measurements were conducted approximately once per month under rain-free conditions. During each campaign, one forest stand was measured per day, and 24-h observations were performed at 1-h intervals for both subplots within that stand. Measurements were conducted approximately once per month under rain-free conditions. Each campaign lasted four consecutive days, with one forest stand measured per day, and 24-h observations were performed at 1-h intervals for both subplots within the measured stand. In total, 12 measurement campaigns were completed during the study period. At each hourly sampling time point, soil respiration was measured sequentially at the three permanently installed PVC collars within each subplot. This yielded a total of 12 campaigns × 4 forest stands × 24 h × 2 subplots × 3 collars = 6,912 raw collar-level respiration observations for Rs and Rh. For model development, the three collar-level measurements within each subplot at each hour were averaged to generate subplot-level hourly respiration values, resulting in 12 × 4 × 24 × 2 = 2,304 observations for subsequent analysis. These averaged subplot-level values were used to represent total soil respiration (Rs) in the undisturbed subplot and heterotrophic respiration (Rh) in the trenched subplot, respectively. These averaged subplot-level values were used to represent total soil respiration (Rs) in the undisturbed subplot and heterotrophic respiration (Rh) in the trenched subplot, respectively. Autotrophic respiration (Ra) was then calculated for each corresponding time point (Ra = Rs − Rh). Thus, machine learning models were trained using subplot-level averaged data rather than collar-level replicate measurements.

Soil environmental variables were continuously monitored in each forest stand using a ZL6 data logger (METER Group, Pullman, WA, USA). Soil temperature (°C), volumetric soil water content (SWC, %), and soil electrical conductivity (EC, dS m −1) were recorded at three depths: 0–10 cm (F1), 10–20 cm (F2), and 20–30 cm (F3). Environmental data corresponding to each soil respiration campaign were extracted and matched with the respiration observations for subsequent analyses. Because the two paired subplots within each forest stand were adjacent and shared the same canopy and stand background, stand structural variables were treated as common stand-level attributes for both subplots, whereas soil respiration measurements differed according to the experimental treatment.

Feature engineering and feature selection

To characterize the spatiotemporal variability of soil respiration and its relationships with environmental drivers, a multidimensional predictor set containing 68 candidate features was constructed. These features were grouped into seven categories: (1) basic features, including soil temperature (°C), soil water content (%), and electrical conductivity (dS/m) at three depths, together with stand structural variables; (2) temporal features, including sine-cosine encodings of month and hour, a day/night indicator, and Julian day; (3) lag features, including 1-, 2-, 3-, and 6-h lagged variables; (4) interaction features, including temperature-moisture coupling terms (°C×%), temperature-day/night interactions (°C), and quadratic terms (°C 2); (5) vertical-profile features, including depth gradients (°C or %), profile variability, and depth ratios; (6) rolling statistical features, including 3-h moving means and standard deviations; and (7) depth-averaged features (°C or %). Detailed definitions and formulas for all features are provided in Supplementary Table S1.

To reduce redundancy and improve model robustness, feature selection was performed in three steps. First, low-variance features (variance < 0.01) were removed. Second, features showing weak marginal relationships with the target variable were excluded using a Pearson correlation threshold of | r | < 0.05. Third, multicollinearity among the remaining predictors was assessed using the variance inflation factor (VIF), and features were iteratively removed until all retained variables satisfied VIF ≤ 10. After feature selection, 47 features were retained for model development (Supplementary Table S1). Feature selection was performed using the training set only, and the selected feature set was then applied to the test set. All numerical predictors were standardized using Z-score normalization. To avoid data leakage, scaling parameters were estimated using only the training data and then applied to the corresponding test data.

Comparison of modeling strategies

Because the dataset included four forest types and three respiration components (Rs, Rh, and Ra), substantial heterogeneity was expected in both predictor distributions and response mechanisms. To evaluate how data grouping influences predictive performance, four modeling strategies were compared. Strategy A was a unified modeling approach in which all forest types and all respiration components were pooled to train a single model. Strategy B was a respiration component-based approach in which three separate models were developed for Rs, Rh, and Ra, each using data from all forest types. Strategy C was a forest type-based approach in which four separate models were developed for the four forest types, each including all respiration components. Strategy D was a fully grouped approach in which forest type and respiration component were crossed, resulting in 12 independent sub-models (4 forest types × 3 respiration components).

To ensure that both the training and test datasets represented the full seasonal range of observations and to reduce sampling bias caused by uneven temporal distribution, the dataset was divided into training (80%) and test (20%) subsets using stratified shuffle splitting by month. During this stage, XGBoost was used as a baseline algorithm to compare the relative performance of different modeling strategies. Final algorithm selection was conducted separately in Section 2.5 based on the optimal strategy identified in this step. Model performance was evaluated using the coefficient of determination (R 2), root mean square error (RMSE), and mean absolute error (MAE). Model overfitting was quantified as Δ R 2 = R 2 _train − R 2 _test. The comparison of modeling strategies was based on random stratified splits with stratification by month. After identifying the optimal strategy, all subsequent algorithm comparisons, benchmark analyses, and interpretability analyses were conducted under this strategy.

Comparison of machine learning algorithms and feature schemes

Based on the optimal modeling strategy identified above, seven machine learning algorithms were evaluated, including both linear and ensemble approaches. The linear models were Linear Regression (Seber and Lee, 2003), Ridge Regression (Hoerl and Kennard, 1970), and Lasso Regression (Tibshirani, 1996), which were used as baseline methods. The ensemble tree-based models were Random Forest (Breiman, 2001), XGBoost (Chen and Guestrin, 2016), LightGBM (Ke et al., 2017), and CatBoost (Prokhorenkova et al., 2018)—together with ecologically interpretable benchmark models. The benchmark models included a simple exponential temperature-response model (Q 10) and a generalized additive model (GAM). The Q 10 model represented soil respiration as: R S = R 10 × Q 10 (T − 10) / 10.

where T was the topsoil temperature at 0–10 cm and the parameters R 10 and Q 10 were estimated by nonlinear least squares. The second benchmark was a generalized additive model (GAM) constructed with penalized smoothing splines for the original soil temperature and moisture predictors, including a tensor-product interaction term between deep-layer temperature and moisture. Because the benchmark models were defined on the original environmental variables only, their input features were identical regardless of feature scheme. To assess the influence of temporal information on model behavior, two feature schemes were designed for each algorithm. Scheme I included all predictors, including temporal encodings, whereas Scheme II excluded temporal variables and retained only environmental and derived environmental predictors. For the ecological benchmark models, the two schemes correspond to the same predictor set because no temporal encodings were included. Hyperparameters of the machine-learning models were optimized using grid search with five-fold cross-validation on the training set. Model performance was evaluated using the coefficient of determination (R 2), root mean square error (RMSE), and mean absolute error (MAE). The within-fold overfitting under random validation was quantified as Δ R 2 = R 2 _train − R 2 _test.

Two validation protocols were applied. The first was a repeated stratified random split, in which the dataset was randomly divided into a training set (80%) and a test set (20%) with stratification by month, and the procedure was repeated with different random seeds to obtain stable estimates. The second was a leave-one-month-out blocked cross-validation, which was designed to assess temporal generalizability and to avoid optimistic estimates caused by temporal autocorrelation. All observations were grouped into monthly blocks. In each fold, one complete monthly block was held out as the test set while the remaining blocks were used for training, and this procedure was repeated for every month. Predictions from all test folds were pooled to compute a single overall blocked R 2, RMSE, and MAE; in addition, per-fold RMSE and MAE values were averaged over folds for reporting. The final model used for ecological interpretation was selected based on the best balance of predictive performance, temporal generalization, and ecological interpretability under Scheme II. CatBoost achieved the most consistent performance under temporally blocked validation when trained without temporal variables, and was therefore selected as the final model for subsequent SHAP and partial-dependence analysis. The use of Scheme II ensured that the mechanistic interpretation reflects environmental controls rather than recurrent seasonal or diel temporal patterns.

Model interpretability analysis

To identify the dominant environmental drivers of soil respiration and to interpret the predictive behavior of the final model, SHAP (SHapley Additive exPlanations) was applied to the CatBoost model trained under Scheme II. Because Scheme II contains no temporal encoding variables, the resulting feature attributions reflect environmentally driven responses rather than seasonal or diel time cues. SHAP quantifies the contribution of each feature to an individual prediction based on Shapley value theory, enabling interpretation at both global and local scales. Global feature importance was calculated as the mean absolute SHAP value across all samples. For each respiration component (Rs, Rh, and Ra), the 10 most important predictors were identified as the dominant environmental drivers. SHAP dependence plots were then used to examine the nonlinear relationships between key predictors and model outputs. To facilitate visualization of interaction patterns in the dependence plots, one interacting feature was selected for each target variable based on its association with the residual variation in the SHAP dependence pattern of the focal feature, and was then used for color mapping. Local SHAP analysis was further performed on three representative samples corresponding to low, median, and high predicted respiration levels. For each sample, the major positive and negative feature contributions were examined to reveal how specific environmental conditions shaped the model prediction. To further quantify joint effects between major predictors, two-dimensional partial dependence plots (2D PDPs) (Friedman, 2001) were generated for the top-ranked features. These plots were used to evaluate the combined influence of environmental variables on predicted soil respiration while averaging over the remaining predictors.

To evaluate whether predictors from different soil depths provided independent predictive information, a controlled depth-wise feature-ablation analysis was additionally conducted under the same leave-one-month-out blocked validation framework. Starting from the full Scheme II predictor set, all variables derived wholly or partly from the 0–10 cm (F1), 10–20 cm (F2), or 20–30 cm (F3) soil layers were removed separately, including corresponding depth-specific and cross-depth derived features. The ablated models were refitted using the same CatBoost configuration, and changes in pooled R 2, RMSE, and MAE relative to the full model were used to quantify the incremental predictive contribution of each depth (Figure 1).

Seasonal dynamics of soil respiration components

Soil respiration exhibited clear seasonal dynamics across the four forest types (Figure 2). Overall, total soil respiration (Rs) and heterotrophic respiration (Rh) followed a pronounced unimodal seasonal pattern, with maxima occurring during the warm growing season (May–August) and minima during the cold season (December–February). In contrast, autotrophic respiration (Ra) showed weaker and more variable seasonal dynamics among forest types.

In the pure coniferous forest (Figure 2A), Rs displayed the largest seasonal amplitude among the four forest types, reaching a maximum in July–August at approximately 4.9 μmol CO 2 m −2 s −1. Rh closely tracked the seasonal pattern of Rs, whereas Ra remained comparatively low throughout the year, varying within a narrower range (approximately 0.2–0.9 μmol CO 2 m −2 s −1) with a moderate summer peak. In the coniferous-broadleaf mixed forest (Figure 2B), Rs peaked earlier, in May, at approximately 2.6 μmol CO 2 m −2 s −1, followed by a relatively stable summer period and a marked decline after September. Rh showed a similar temporal pattern, whereas Ra exhibited greater seasonal complexity, with relatively higher values in May and October and lower values during midsummer. In the economic forest (Figure 2C), Rs increased from winter to a peak in May and then gradually declined toward the end of the growing season. Rh again showed dynamics closely aligned with those of Rs. By contrast, Ra displayed weak seasonality and remained within a relatively narrow range (approximately 0.2–0.7 μmol CO 2 m −2 s −1) throughout the study period. In the secondary broadleaf forest (Figure 2D), Rs showed pronounced growing-season activity, with high values from May to August and very low rates during the dormant period from November to April. Rh generally mirrored the temporal pattern of Rs, while Ra was mainly active during the growing season, peaking in June and approaching zero during the dormant season.

Taken together, these results indicate that Rs and Rh shared consistent seasonal trajectories across forest types, whereas Ra exhibited weaker magnitude and greater temporal variability. The timing and amplitude of seasonal peaks differed among forest types, but the overall seasonal signal was robust across the study region.

Model performance under different modeling strategies and algorithms

The four modeling strategies showed substantial differences in predictive performance (Table 1). Strategy A, in which all forest types and respiration components were pooled into a single model, yielded the lowest test R 2 (0.508 ± 0.034) and relatively high overfitting (Δ R 2 = 0.042 ± 0.044). Strategy C, which grouped data only by forest type, performed similarly poorly (test R 2 = 0.478 ± 0.089), indicating that separating observations by vegetation type alone did not effectively reduce dataset heterogeneity. In contrast, Strategy B, which modeled Rs, Rh, and Ra separately, markedly improved predictive performance (test R 2 = 0.795 ± 0.198) while reducing overfitting (Δ R 2 = 0.068 ± 0.078). Under this strategy, Rs and Rh both achieved high predictive accuracy, whereas Ra remained less predictable. Strategy D, which fully crossed forest type and respiration component, achieved the highest mean test performance (test R 2 = 0.849 ± 0.178), but required 12 independent sub-models and showed a higher degree of complexity and overfitting than Strategy B. Considering predictive accuracy, model generalization, and structural simplicity together, Strategy B was selected as the optimal modeling framework.

Performance differences were also evident among machine learning algorithms (Supplementary Table S2). The simple exponential Q 10 model performed poorly across all respiration components, with test R 2 values below 0.18 even under random validation and near-zero or negative R 2 under blocked validation, indicating that a temperature-only response is insufficient for this subtropical mountainous forest. The GAM, as an ecologically interpretable benchmark, performed substantially better than Q 10, reaching random-split R 2 values of 0.58 for Rs and 0.56 for Rh, but its blocked-validation R 2 was limited to 0.33 and 0.29, respectively. Linear models showed intermediate predictive ability under random splits (R 2 = 0.59–0.72 for Rs and Rh, 0.10–0.29 for Ra) but their temporal generalization was poor, with some blocked R 2 values falling below zero. In comparison, ensemble tree-based models consistently outperformed both ecological benchmarks and linear models, with random-split R 2 generally exceeding 0.90 for Rs and Rh and reaching 0.52–0.57 for Ra. However, their advantage narrowed substantially under leave-one-month-out blocked validation, confirming that part of the machine learning predictive power originates from recurrent seasonal and diel temporal patterns rather than from environmental responses alone.

Comparing the two feature schemes, excluding temporal variables produced little change in random-split accuracy (R 2 differences ≤ 0.02 in most cases), but generally improved temporal generalization under blocked validation. Notably, CatBoost trained without temporal variables (Scheme II) achieved the most consistent performance blocked-validation R 2 for all three respiration components (Rs: 0.475; Rh: 0.544; Ra: 0.081), outperforming random forest (Rs: 0.290; Rh: 0.371; Ra: −0.057) and the other ensemble models under the same protocol. Considering both temporal generalizability and ecological interpretability, CatBoost under Scheme II was selected as the final model for subsequent mechanistic analyses. Consistent with the component-based modeling strategy, Ra remained the least predictable component across all models and validation protocols, reflecting the influence of unmeasured plant physiological processes on autotrophic respiration.

Model interpretability and dominant environmental drivers

SHAP analysis of the final CatBoost model revealed clear differences in the dominant predictors among the three soil respiration components and showed that the feature scheme strongly influenced model interpretation (Figure 3). Under Scheme I, which included temporal variables, seasonal predictors dominated the importance ranking for Rs and Rh (temporal importance exceeded 405% of the total SHAP importance). In particular, the cosine-transformed monthly variable (month_cos) had the highest SHAP importance for both components, followed by the deep-layer temperature–moisture interaction (F3_TxWC) for both components. This pattern indicates that, when temporal variables were included, the models for Rs and Rh relied heavily on seasonal structure, despite the additional contribution of subsurface hydrothermal coupling. In contrast, Ra showed only weak dependence on temporal variables under Scheme I: although month_cos ranked first, its importance was much lower (0.038), and moisture ratio (WC_ratio), daytime temperature variables, and hour all contributed to a relatively diffuse importance structure, suggesting that autotrophic respiration was less controlled by simple seasonal periodicity and more closely associated with short-term environmental variation (Figure 3A).

After temporal variables were excluded (Scheme II), the importance structure changed markedly and environmental predictors emerged as the primary drivers (Figure 3B). For both Rs and Rh, the interaction between deepest monitored soil temperature and moisture (F3_TxWC) was the most important predictor, followed by temperature gradient (TEM_gradient), mid-layer temperature–moisture interaction (F2_TxWC), and temperature ratio (TEM_ratio). For Rh, the shallow temperature–moisture interaction (F1_TxWC) also ranked among the top predictors. This pattern indicates that the predictive behavior of the models for Rs and Rh was dominated by subsurface hydrothermal conditions and their vertical organization within the soil profile. By comparison, Ra was primarily associated with deep-layer daytime temperature (F3_TEM_day), temperature gradient (TEM_gradient), moisture ratio (WC_ratio), and temperature ratio (TEM_ratio), while the importance values of individual predictors were more evenly distributed. Relative to Rs and Rh, Ra therefore showed a weaker dependence on any single predictor and a more diffuse importance structure, highlighting the roles of soil thermal stratification and moisture balance in regulating autotrophic respiration.

Overall, these results indicate that the removal of temporal variables substantially improved ecological interpretability. Rs and Rh were mainly associated with hydrothermal coupling in deeper soil layers and with vertical temperature structure, whereas Ra was more strongly related to deep-layer diurnal temperature variation and soil moisture balance.

Nonlinear responses and interaction effects of key environmental drivers

SHAP dependence analysis and two-dimensional partial dependence plots revealed pronounced nonlinear responses and interaction effects among hydrothermal variables, with clear contrasts among respiration components (Figures 4, 5). For total soil respiration (Rs), the deep-layer temperature–moisture interaction term (F3_TxWC) exhibited a distinct threshold-like pattern. At low values of this hydrothermal index, SHAP contributions were generally negative, indicating suppression of predicted respiration. Once the index exceeded approximately 6, SHAP values increased rapidly, indicating a strong positive contribution to model output. The temperature gradient (TEM_gradient) showed a turning point around −2, above which SHAP values declined sharply, suggesting that stronger vertical thermal stratification was associated with reduced predicted Rs. The mid-layer temperature–moisture interaction (F2_TxWC) exhibited a threshold at approximately 4, with SHAP contributions increasing markedly once this value was exceeded (Figure 4A).

Heterotrophic respiration (Rh) showed response patterns broadly similar to those of Rs (Figure 4B). In contrast, autotrophic respiration (Ra) showed weaker dependence on hydrothermal interaction terms and was more strongly associated with deep-layer daytime temperature (F3_TEM_day) (Figure 4C). When F3_TEM_day exceeded approximately 8 °C, SHAP values increased sharply, indicating a persistent positive contribution. The temperature gradient showed a similar response pattern to that observed for Rs and Rh, with negative contributions above the turning point. The moisture ratio (WC_ratio) exhibited a more complex response: SHAP values declined sharply when WC_ratio exceeded approximately 0.6, remained relatively stable over the intermediate range, and then increased sharply as WC_ratio approached 1.1, suggesting nonlinear moisture regulation of autotrophic respiration.

The two-dimensional partial dependence plots further highlighted joint effects among key predictors (Figure 5). For Rs and Rh, predicted respiration increased strongly with increasing deep-layer temperature–moisture coupling (F3_TxWC), but this positive effect was attenuated by high temperature gradients. When TEM_gradient increased, the promoting effect of F3_TxWC on predicted respiration was substantially suppressed, indicating that strong vertical thermal stratification constrained the positive hydrothermal response (Figures 5A,B). A similar joint pattern was observed for Ra, with deep-layer daytime temperature (F3_TEM_day) acting as the primary promoting factor. High F3_TEM_day enhanced predicted Ra, whereas increasing TEM_gradient again exerted a downward constraint. Overall, the response surface for Ra was weaker and flatter than those for Rs and Rh (Figure 5C).

Local SHAP analysis also demonstrated that the direction and magnitude of predictor effects varied among environmental clusters (Supplementary Figure 1). For Rs, representative samples from Cluster 1 and Cluster 2 (low and intermediate respiration) were dominated by strong negative contributions from F3_TxWC. In Cluster 3 (high respiration), TEM_gradient and Tem_ratio exerted prominent positive effects, while F3_TxWC still contributed negatively (Supplementary Figure 1A). Rh exhibited a similar cluster-dependent pattern: Cluster 3 showed strong positive contributions from TEM_gradient and F1_TxWC, alongside negative contributions from F3_TxWC (Supplementary Figure 1B). In contrast, Ra displayed markedly different SHAP patterns across clusters. For Cluster 1, predictors including F2_T2 and F3_TEM_day were predominantly negative. In Cluster 2, F3_TEM_day had a large positive contribution, whereas Tem_ratio showed negative effects. Cluster 3 was characterized by positive contributions from TEM_gradient and WC_ratio. Notably, the interpretation for Ra should be treated cautiously because the model achieved poor prediction performance. Taken together, these results show that different soil respiration components were associated with distinct nonlinear response patterns and interaction structures. Rs and Rh were most strongly related to deepest monitored soil hydrothermal coupling and vertical temperature organization, whereas Ra was more closely associated with lagged thermal conditions and moisture balance, though the weak predictive power limits robust inference for Ra.

Implications of modeling strategy and machine learning framework for soil respiration prediction

This study demonstrates that machine learning performance in predicting soil respiration in subtropical mountainous forests depends strongly not only on the algorithm used, but also on how data heterogeneity is handled and how model generalization is evaluated. The comparison among modeling strategies showed that component-based modeling provided the most effective framework. Grouping data by respiration component substantially improved predictive performance relative to both unified modeling and forest type-based modeling, while avoiding the high structural complexity of the fully grouped approach. Although the fully grouped strategy achieved a slightly higher mean test R 2 than the component-based strategy, it required 12 independent sub-models and exhibited a larger training–test performance gap. The component-based strategy therefore provided a more favorable balance between predictive accuracy, generalization, and model parsimony. These results suggest that heterogeneity in the present dataset was more strongly associated with differences among respiration components than with differences among forest types. Partitioning Rs, Rh, and Ra before model fitting may therefore allow the algorithms to better resolve component-specific response patterns that would otherwise be obscured in pooled datasets.

This advantage is also ecologically meaningful. Rs and Rh showed similar predictive behavior and were strongly linked to hydrothermal conditions, whereas Ra exhibited lower predictability and a more diffuse driver structure. Such differences are consistent with the distinct biophysical controls on respiration components: heterotrophic respiration is more directly constrained by environmental conditions and decomposition processes, whereas autotrophic respiration is additionally regulated by plant carbon allocation, root activity, and physiological status (Kuzyakov and Gavrichkova, 2010; Tang et al., 2019). The superior performance of component-based modeling therefore reflects not only a statistical improvement, but also a better correspondence between model structure and ecological process heterogeneity. Importantly, the particularly low predictive performance for Ra under temporally blocked validation further indicates that its temporal variation cannot be adequately explained by the measured soil environmental variables alone.

The comparison with ecologically interpretable benchmark models further refined our assessment of machine-learning performance. The simple exponential Q 10 model showed limited predictive ability, indicating that temperature alone was insufficient to represent soil respiration dynamics in these subtropical mountainous forests (Davidson et al., 2006). This is consistent with increasing evidence that respiration responses to temperature are strongly conditioned by water availability and nonlinear hydrothermal interactions (Moyano et al., 2013). The GAM, which explicitly accommodated nonlinear temperature and moisture responses and their interaction, substantially improved upon the Q 10 model, but remained less accurate than the ensemble tree-based models, particularly for Rs and Rh. Under repeated random splits, ensemble machine-learning models achieved high predictive accuracy and clearly outperformed both ecological benchmarks and linear models. However, this advantage narrowed considerably under leave-one-month-out blocked validation. The decline in performance under temporal blocking indicates that part of the high accuracy obtained from random validation arose from the ability of machine-learning models to exploit recurrent seasonal and diel structures rather than exclusively from their ability to represent transferable environmental response relationships. This distinction is important for both model evaluation and ecological interpretation. Temporal variables produced little improvement in random-split predictive accuracy, but became highly influential in SHAP rankings when they were included, particularly for Rs and Rh. Their high attribution importance therefore did not necessarily imply that they provided substantial independent predictive information; instead, the models could use them as efficient proxies for recurrent seasonal environmental conditions. Consistent with this interpretation, excluding temporal variables generally improved temporal generalization under blocked validation and shifted model attribution toward measured hydrothermal predictors. Among the tested algorithms, CatBoost trained without temporal variables (Scheme II) showed the most consistent performance under temporally blocked validation, achieving the highest blocked R 2 for Rs, Rh, and Ra. It was therefore selected as the final model for subsequent ecological interpretation. This choice prioritizes transferable environmental relationships over model performance derived from explicit temporal cues and provides a more conservative basis for SHAP and partial-dependence analyses.

Nonlinear regulation mechanisms of hydrothermal conditions and their vertical structure on soil respiration

Beyond model performance, the interpretability analyses revealed that soil hydrothermal conditions and their vertical organization were the dominant controls on respiration dynamics, particularly for Rs and Rh. After temporal variables were excluded, the interaction between deep-layer soil temperature and moisture emerged as the most important predictor for both components, together with variables describing vertical temperature ratio and gradient. This indicates that soil respiration in these subtropical forests was regulated not simply by surface conditions, but by coupled hydrothermal processes across the soil profile. The importance of deepest monitored soil hydrothermal coupling is ecologically plausible because soil CO 2 efflux integrates biological activity from multiple depths, including root respiration, microbial decomposition, and gas transport processes (Tang et al., 2005; Davidson et al., 2006; Subke and Bahn, 2010). In the present study, both SHAP dependence plots and 2D partial dependence analyses showed pronounced nonlinear responses. For Rs and Rh, the deep-layer temperature-moisture interaction exhibited a threshold-like transition, with its contribution to predicted respiration increasing rapidly once the hydrothermal index reached a relatively high range. This suggests that temperature and moisture acted jointly rather than independently, consistent with previous evidence that soil respiration is constrained by hydrothermal coupling rather than by single environmental factors alone (Moyano et al., 2013; Zhao et al., 2025).

Vertical thermal structure also contributed substantially to Rs and Rh predictions. Changes in temperature gradient modified the positive association between deep-layer hydrothermal conditions and predicted respiration, indicating that conditions at a single depth may not fully represent the response of the whole soil profile. Importantly, however, the high predictive importance of variables measured at 20–30 cm does not imply that CO 2 produced at this depth dominates surface efflux. Deep-layer temperature may also provide a smoother signal of whole-profile thermal conditions or seasonal processes. To test whether deep-layer variables contained independent predictive information, we conducted a depth-wise feature-ablation analysis under temporally blocked validation (Supplementary Table S3). Removing predictors related to the 20–30 cm layer reduced pooled R 2 by 0.167 for Rs and 0.143 for Rh and increased RMSE by approximately 11.8% and 10.7%, respectively, whereas removing shallower-layer predictors did not reduce model performance. These results indicate that deep-layer measurements provided additional predictive information beyond the shallower layers, but they should not be interpreted as evidence for the vertical source of respired CO 2.

In contrast, Ra showed a weaker and more diffuse predictor structure. Deep-layer daytime temperature, temperature gradient, moisture ratio, and temperature ratio were among its important predictors, but their contributions were more evenly distributed. Moreover, removing F3-related predictors reduced pooled R 2 by only 0.008, indicating little independent predictive contribution from the deep layer for Ra. This weak environmental predictability is consistent with the more complex regulation of autotrophic respiration by plant physiological processes, including photosynthate transport, root metabolism, phenology, and carbon allocation (Bahn et al., 2008; Kuzyakov and Gavrichkova, 2010; Guan et al., 2022). Overall, Rs and Rh were more strongly associated with whole-profile hydrothermal organization, whereas Ra was less well captured by soil environmental variables alone.

Limitations and future perspectives

Several limitations should be acknowledged. First, the lower predictive accuracy for autotrophic respiration (Ra) indicates that soil environmental variables alone were insufficient to capture the biological controls of this component. Important drivers such as canopy photosynthesis, root phenology, and carbon allocation were not explicitly included (Ryan and Law, 2005; Tang et al., 2019). Second, respiration partitioning based on trenching may introduce uncertainty because root exclusion can alter soil moisture, substrate availability, and microbial processes (Hanson et al., 2000; Subke et al., 2006). Third, the study was restricted to four forest types within a single subtropical mountainous region and covered only a little more than one year, which limits the evaluation of interannual variability, extreme events, and broader model transferability. Moreover, one of the selected types was an intensively managed citrus orchard, and its management practices (e.g., fertilization, irrigation) may influence soil respiration independently of tree species or climate. Consequently, the observed differences among forest types may partly reflect management effects, which could not be fully separated from vegetation or environmental effects in this study. Future studies should directly quantify management intensity to address this issue. In addition, because soil respiration was measured approximately once per month and only under rain-free conditions, the present results primarily characterize hydrothermal regulation during the sampled rain-free periods and may not capture post-rain pulses, dry-to-wet transitions, or high-moisture diffusion limitation. Finally, although interpretable machine learning improved ecological inference, the identified relationships remain associational rather than strictly causal. Future studies should therefore incorporate longer-term observations, broader spatial gradients, and additional plant physiological variables, and should combine machine learning with experimental or process-based approaches to strengthen mechanistic understanding.