Explainable Machine Learning Reveals Persistent Carbon Sink in Xishuangbanna Tropical Forests Under Future Climate Scenarios

1. Introduction

Tropical rainforest ecosystems serve as natural buffers against global climate change, playing a dual core role in the global carbon and water cycles Almeida. On one hand, tropical rainforests sequester approximately 1.6±0.5 Pg of carbon annually, constituting a critical carbon sink for the planet Xu. On the other hand, the intense evapotranspiration from tropical trees acts as a biological “pump”, transporting deep soil water to the atmosphere, which is irreplaceable for sustaining regional precipitation patterns and the global hydrological cycle Luiz. However, global warming and deforestation are diminishing the carbon sink strength and evapotranspiration of major rainforests, such as the Amazon Doughty. Rising temperatures not only increase plant respiratory carbon loss Du, and alter species composition Zhu, but also exacerbate atmospheric aridity and trigger extreme climatic events like heatwaves and droughts, thereby impacting rainforest ecological functions. For instance, soil moisture depletion due to excessive evaporation can induce water stress, potentially reducing photosynthetic rates Guan, and altering precipitation regimes and increased frequency of extreme droughts can lead to tree mortality Hubau. Extensive observational and modeling studies warn that if current trends persist, the functional degradation of rainforests driven by warming could transform them from carbon sinks into carbon sources by the mid-21st century Katharyn.

The Xishuangbanna tropical rainforest shares core characteristics with typical tropical rainforests, such as community structure and species composition Zhu. However, located at the northernmost edge of the Indo-Malayan rainforest region Zhou, it exhibits significant marginal and transitional properties. Consequently, it is more sensitive to climate change and is often regarded as a critical indicator for global rainforest ecosystems Alistair. Notably, against the backdrop of widespread functional decline in global tropical rainforests, the carbon sequestration capacity of the pristine Xishuangbanna rainforest ecosystem has maintained a steady annual increase of approximately 3.4% over the past two decades Liu. This observed phenomenon presented by Xishuangbanna offers a unique case for re-examining the responses of global tropical rainforests to future climate change.

The response of tropical forest carbon and water fluxes to environmental variables is highly complex. Previous studies have attributed the enhanced carbon sink in the Xishuangbanna rainforest to improved carbon use efficiency, highlighting increased dry-season solar radiation as a primary driver Liu. However, the enhancement of carbon sequestration capacity likely stems from more intricate physical processes driven by the synergistic effects of multiple environmental factors. Beyond the direct regulation of forest carbon uptake via photosynthesis by solar radiation Aguilos, meteorological factors such as air temperature, humidity, and wind speed also significantly influence ecosystem CO₂ exchange Wang. For instance, prolonged dry seasons Esquivel may induce drought stress, yet simultaneously lengthen the carbon absorption window, leading to divergent potential impacts on carbon sink capacity Dai. Warmer soils may also enhance soil carbon emissions Zhang. Furthermore, under global climate change, increased evaporative demand from a warmer atmosphere may accelerate soil moisture depletion, leading to surface drying Leander. Deforestation drastically reduces land surface evapotranspiration (ET), decreasing regional atmospheric moisture content and subsequently diminishing rainfall Qin. These changes not only affect the water balance of rainforests but may also fundamentally compromise rainforest functionality by altering vegetation structure Xu.

Although long-term in situ monitoring data from eddy covariance (EC) systems provide a foundation for analyzing the environmental responses of carbon and water fluxes Chang, disentangling the complex, non-linear influences of environmental variables on these fluxes remains challenging, specifically net ecosystem exchange (NEE) and evapotranspiration (ET). This is due to the high-dimensional interactions among variables, which obscure the identification of key drivers. Moreover, while physics-based models like those from CMIP6 offer crucial frameworks for projecting future ecosystem trajectories, their simulations of key ecological parameters (e.g., NEE, ET) at regional scales carry substantial uncertainties Shen, For example, CMIP6 models generally overestimate global ET, showing only a 55% consensus rate with the ERA5 dataset Sun, Direct reliance on uncorrected CMIP6 outputs can thus lead to biased conclusions.

Therefore, this study proposes an integrated analytical framework that combines explainable machine learning (ML), structural equation modeling (SEM), and CMIP6 data. We apply this framework to a two-decade-long dataset of carbon-water fluxes and multi-source environmental data from the Xishuangbanna tropical rainforest. The framework proceeds in three steps: First, we employ Pearson correlation coefficients to filter parameters within similar groups, then use explainable ML to quantitatively parse the independent contributions and interactive effects of multiple environmental drivers on NEE and ET. Second, based on the identified key variables, we construct SEMs to solve for the path coefficients representing their directed influences on NEE and ET. Finally, we predict the future trajectories of NEE and ET under different climate scenarios by applying the validated ML models to bias-corrected CMIP6 environmental variable data.

2. Materials and Methods

2.1. Study Site and Data Sources

The study was conducted in the Xishuangbanna Tropical Rainforest Nature Reserve, Yunnan Province, China (Figure 1), The reserve harbors four primary forest types: tropical rainforest, riparian tropical monsoon forest, tropical low-mountain evergreen broad-leaved forest on slopes, and tropical seasonal moist forest Zhu. A comprehensive soil-rainforest-atmosphere flux observation station was established within the reserve. This station integrates an eddy covariance (EC) flux system, a conventional meteorological measurement system, and a soil temperature and moisture monitoring network Liu, The station is situated in a relatively flat area within an approximately 2-km-long valley Zhang.

Carbon-water flux data and associated meteorological and soil data for Xishuangbanna from 2003 to 2015 were obtained from a publicly available dataset hosted by the National Ecosystem Science Data Center, China (http://www.nesdc.org.cn). This dataset comprises 28 variables, including evapotranspiration (ET), net ecosystem exchange (NEE), solar radiation, and air temperature above the canopy, with a temporal resolution of 30 minutes, amounting to over 200,000 records.

Meteorological data for the period 2015 to 2022 were acquired from the same center under a data use agreement . This dataset includes six categories of variables: radiation, air temperature, and relative humidity, among others, with a daily temporal resolution.

Climate projection data were sourced from the publicly accessible Coupled Model Intercomparison Project Phase 6 (CMIP6) archive (https://aims2.llnl.gov/search). We selected two global climate models (GCMs): MPI-ESM1-2-HR (at~1.875°×1.865°resolution) and BCC-CSM2-MR (at~1.125°×1.12°resolution). For each model, we utilized historical simulation data (covering 2003–2014) and future projections under four Shared Socioeconomic Pathways (SSPs): SSP1-2.6, SSP2-4.5, SSP3-7.0, and SSP5-8.5.

Appendix Table1 presents a cross-reference list of the abbreviations and full meanings for the variables used from the flux observation station data.

2.2. Study Framework

As illustrated in Figure 2, the analytical framework was divided into four parts: data collection and preprocessing, comparison of machine learning models, SHAP analysis, and prediction of future carbon sink and evapotranspiration capacity. Firstly, in situ carbon and water flux data, soil environmental data, and atmospheric environmental data were integrated. Data cleaning was performed to remove null values and anomalous data (e.g., -9999). After standardization, Pearson correlation coefficients (AppendixA(1)), were calculated among all variables. From variables of the same type, the one with the highest correlation with the target variable (ET or NEE) was selected as the representative feature for subsequent analysis; Subsequently, six machine learning models—RF, ELM, SVR, MLP, XGBoost, and KNN—were used to fit the dataset. The representative features served as input, with ET and NEE as separate outputs. All models employed Bayesian hyperparameter optimization to ensure a fair comparison. Model performance was evaluated using the root mean square error (RMSE), mean absolute error (MAE), and coefficient of determination (R²) (AppendixA formula $\left(3\right)$, (2), (4)); The optimal model was selected to calculate sample SHAP values. Stratified sampling was used to accelerate the computation, and Jensen-Shannon (JS) divergence was applied to measure the consistency between the distribution of the sampled subset and the full dataset (JS < 0.05 indicated that the sample could represent the overall population). SHAP analysis and Partial Dependence Plots (PDP) were used to analyze the local contribution of individual variables to the prediction results and the interactions between variables. Feature importance ranking plots were used to screen key variables and analyze their monthly importance changes. Structural Equation Modeling (SEM) was employed to solve for the path coefficients of the impact of key variables on ET and NEE, clarifying the mechanism of environmental influences on the carbon sink and water cycle; Finally, based on environmental data from 2003 to 2022, probability distribution correction was applied to the environmental variables in CMIP6. Data after 2030 were input into the machine learning model to predict future ET and NEE, and the results were compared and evaluated with the original CMIP6 outputs.

2.3. SHAP Analysis and SEM-Based Causal Analysis

SHapley Additive exPlanations (SHAP) is a widely used analytical tool for interpreting the prediction mechanism of machine learning models. It calculates the marginal contribution (the SHAP value) of each variable to the prediction outcome for every individual sample Scott. The formula for calculating the SHAP value is as follows:

$$\mathrm{SHAP}\left({\mathrm{x}}_{\mathrm{i}}\right)=\sum _{\mathrm{S}\subseteq \left\{{\mathrm{x}}_1,...,{\mathrm{x}}_{\mathrm{n}}\right\}\backslash {\mathrm{x}}_{\mathrm{i}}}{\displaystyle \frac{\left|\mathrm{S}\right|!(\mathrm{n}-|\mathrm{S}|-1)!}{\mathrm{n}!}}\left[\mathrm{f}\right(\mathrm{S}\cup \{{\mathrm{x}}_{\mathrm{i}}\left\}\right)-\mathrm{f}\left(\mathrm{S}\right)]$$

(1)

Where S is a subset of all features excluding xi, n is the total number of features, and f(x) is the model output. ${\displaystyle \frac{\mid S\mid !(n-\mid S\mid -1)!}{n!}}$ used to adjust the weights of different features to ensure fairness in the model.

The statistical characteristics of SHAP values were analyzed through visualization to interpret interactions among variables. Common graphical tools include the SHAP Summary Plot (SP), which ranks features by their global importance, and Partial Dependence Plots (PDPs), which illustrate the marginal effect of one or two variables with different values on the model predictions.

Structural Equation Modeling (SEM) is a causal analysis method commonly employed in ecology Grace, It constructs candidate models based on theoretical understanding, which are then iteratively validated and optimized to derive a path diagram that reflects causal relationships. This technique integrates factor analysis and path analysis to examine and estimate complex causal relationships among multiple variables Kenneth. SEM evaluates whether the proposed network of influence paths fits the observed data through fit indices, moving beyond merely testing individual coefficients. Explainable machine learning provides powerful, data-driven prior insights for establishing causal testing models by pre-screening and prioritizing the most critical driving factors Sun. Therefore, we combined SHAP analysis with SEM to uncover the key response mechanisms of environmental variables on ET and NEE.

2.4. Stratified Sampling and Jensen-Shannon Divergence

The computation of SHAP values for the entire large data size imposes a substantial computational burden. Therefore, a stratified sampling approach was adopted to alleviate this load. Following the rule of randomly sampling 100 data points per month, a subset of 1,200 samples was selected from the massive dataset to form a representative sub-dataset. A schematic of the sampling procedure is provided in Appendix Figure 1a.

The Jensen-Shannon (JS) divergence is a method commonly used in machine learning to measure the similarity between two probability distributions Pedro, It is derived from the Kullback-Leibler (KL) divergence and addresses the asymmetry issue inherent in the KL divergence (KL(P||Q)≠KL(Q||P), providing a more accurate reflection of distributional similarity. The formulas are defined as follows:

$$\mathrm{JS}\left(\mathrm{P}\right|\left|\mathrm{Q}\right)={\displaystyle \frac{1}{2}}\mathrm{KL}\left(\mathrm{P}\right(\mathrm{x}\left)\right||{\displaystyle \frac{\mathrm{P}\left(\mathrm{x}\right)+\mathrm{Q}\left(\mathrm{x}\right)}{2}}+{\displaystyle \frac{1}{2}}\mathrm{KL}\left(\mathrm{Q}\right(\mathrm{x}\left)\right||{\displaystyle \frac{\mathrm{P}\left(\mathrm{x}\right)+\mathrm{Q}\left(\mathrm{x}\right)}{2}})$$

(2)

$$\mathrm{KL}\left(\mathrm{P}\right|\left|\mathrm{Q}\right)=\sum \mathrm{P}\left(\mathrm{x}\right)\log {\displaystyle \frac{\mathrm{P}\left(\mathrm{x}\right)}{\mathrm{Q}\left(\mathrm{x}\right)}}$$

(3)

Here, Q represents the original sample population, with Q(x) denoting its probability distribution. P represents the stratified sampling subset, with P(x) denoting its probability distribution. The value range of JS divergence is [0, 1]. A value closer to 0 indicates greater similarity between the two distributions, while a value of 1 indicates complete dissimilarity.

2.5. Probability Distribution Correction of CMIP6 Data

In climate change impact studies, direct use of raw outputs from CMIP6 global climate models is subject to systematic biases. Therefore, applying bias correction to the data is a necessary step for generating reliable regional climate information Syukri. The core objective is to adjust the probability distribution of the CMIP6-simulated variables to align with the distribution of observed data.

In this study, quantile delta mapping was employed to achieve this goal, with the procedure illustrated in Appendix Figure 1b. A transfer function was constructed by matching quantiles between historical observed data and CMIP6-simulated historical data. For any future CMIP6 data point, its quantile position τ within the future simulated data series is first determined. Based on τ, the corresponding value from the historical CMIP6 simulations is identified. This value is then mapped to its observational counterpart. Finally, the difference between the future simulated value and its corresponding historical simulated value at quantile τ is added Qin. The formula is expressed as:

$$\begin{array}{rr}& {\mathrm{x}}_{\mathrm{corr}}={\mathrm{F}}_{\mathrm{obs}}^{ -1}\left({\mathrm{F}}_{\mathrm{fut}}\left(\mathrm{x}\right))+[\mathrm{x}-{\mathrm{F}}_{\mathrm{hist}}^{ -1}\right({\mathrm{F}}_{\mathrm{fut}}\left(\mathrm{x}\right)\left)\right]\\ & \end{array}$$

(4)

This method corrects the probability distribution of the simulated values to match the observations while preserving the relative climate change signal predicted by the model. Here, xcorris the bias-corrected future value, x is the raw future CMIP6 simulated value, F denotes the cumulative distribution function, and F−1 denotes the quantile function. The subscripts obs, hist, and fut denote historical observations, historical CMIP6 simulations, and future CMIP6 simulations, respectively.

3. Results and Analysis

3.1. Representative Feature Selection

Pearson correlation coefficients were calculated among the 23 observed variables. Figure 3 presents the correlation matrix heatmap. The 21 variables (excluding ET and NEE) were categorized into nine types: air temperature, air humidity, vapor pressure, wind speed, atmospheric pressure, radiation, soil temperature, soil moisture, and precipitation. The heatmap reveals high correlation coefficients among variables within the same group, indicating strong linear dependencies.

To reduce multicollinearity and model complexity, a representative feature was selected from each variable type for subsequent modeling. The selection criterion was based on the feature's absolute correlation coefficient with the target variable (ET or NEE). The variable exhibiting the strongest correlation within each group was chosen. A correlation threshold of |ρ| < 0.1 was applied to exclude features considered uncorrelated with the target. Based on this process, the following representative features were selected:

For ET: ATATC (ρ = 0.35), AHATC (ρ = -0.36), NGVP (ρ = 0.20), WSATC (ρ = 0.35), SR (ρ = 0.50), and PR (ρ = 0.08).

For NEE: ATATC (ρ = -0.31), AHATC (ρ = 0.29), NGVP (ρ = -0.20), WSATC (ρ = -0.40), AP (ρ = -0.15), SR (ρ = -0.54), STL2 (ρ = -0.17), and PR (ρ = -0.02).

The ET dataset and NEE dataset are constructed based on the representative features and target variables.

3.2. Optimal Model Selection for ML

The ET and NEE datasets were split into training and test sets with an 8:2 ratio. The training set was used for model fitting and optimization, while the test set was reserved for evaluating model performance. As shown in Figure 4, XGBoost outperformed the other models for predicting both ET and NEE. For ET prediction, XGBoost achieved a mean absolute error (MAE) of 0.025, a root mean square error (RMSE) of 0.037, and a coefficient of determination (R²) of 0.403. For NEE prediction, it achieved an MAE of 0.171, an RMSE of 0.240, and an R² of 0.464, indicating a reasonably good fit. To ensure a fair comparison and avoid bias in model selection due to parameter settings, Bayesian hyperparameter optimization was employed during training to automatically search for the optimal performance of each model. The search spaces for the hyperparameters are detailed in Appendix A Table 2. Based on its superior performance, XGBoost was selected as the optimal model for all subsequent explainability analysis. The optimized hyperparameter set identified for XGBoost was used for its final training.

3.3. Explainable ML Analysis

Based on the stratified sampling method, a subset of 1,200 samples was selected, and their SHAP values were computed to analyze the marginal contributions of environmental variables to the predictions of ET and NEE. Figure 5 illustrates the distribution of the sampled data compared to the full dataset. The distribution of ET is markedly right-skewed with a long tail, exhibiting a skewness of 2.17. The Jensen-Shannon (JS) divergence between the sampled subset and the full ET dataset is 0.05, indicating a high degree of distributional similarity. In contrast, the distribution of NEE approximates a normal distribution, with a skewness of only -0.07. The JS divergence for the NEE distributions is 0.0426, also confirming a high degree of similarity. Therefore, the stratified sample of 1,200 instances is representative of the original dataset.

Based on the stratified sample, SHAP values were calculated to rank the marginal contributions of each variable to the model predictions, Figure 6 presents this ranking. For ET prediction, the four most influential variables are Solar Radiation (SR), Air Temperature Above the Canopy (ATATC), Near-Ground Vapor Pressure (NGVP), and Wind Speed Above the Canopy (WSATC). For NEE prediction, the top four contributors are SR, NGVP, ATATC, and Soil Temperature at 10 cm depth (STL2). The key variables for both fluxes show significant overlap. Notably, the marginal contribution of SR to both ET and NEE predictions far exceeds that of all other variables, while Precipitation (PR) shows the smallest marginal contribution (Figure6b, d). Furthermore, Figure 6a, bindicate that, aside from SR and PR which exhibit a generally positive contribution to ET, other environmental variables do not show a consistent directional bias (positive or negative) in their marginal contributions to the target predictions. To elucidate the ecological significance of this ranking, we interpret the results in the context of established ecophysiological processes. In rainforest ecosystems, SR serves as the primary energy source, driving both the water and carbon cycles Yu. This driving force operates through two principal pathways: i) the direct effect of SR on vegetation physiology, and ii) indirect effects mediated through SR-induced changes in atmospheric and soil thermal conditions.

Firstly, SR directly governs evaporation and photosynthesis, thereby influencing ET and NEE Marcos. The atmosphere acts as a medium for radiation transfer; its temperature and humidity are directly affected by radiation Wazwaz, Vapor pressure (VP), a latent variable derived from temperature and relative humidity, influences stomatal aperture in rainforest plants. This, in turn, affects transpiration and photosynthesis, thereby impacting ET and NEE Qie, forming two indirect influence pathways. While SR influences wind speed at larger spatial scales Lin, its effect on local wind speed variations is minimal. Nevertheless, wind facilitates the diffusion of water vapor and CO₂ in the atmosphere, affecting both ET and NEE Li. Soil temperature (ST) increases both by directly absorbing radiation and through heat transfer from the warmed atmosphere Wu, ST subsequently influences plant growth and soil microbial respiration, thereby affecting NEE Zhu. Furthermore, ET itself influences NEE. Transpiration cools leaf surfaces and regulates plant water status, which can promote carbon assimilation Pauline. ET also provides feedback to atmospheric humidity and temperature, and ST can be lowered through evaporative cooling effects, representing another indirect pathway affecting NEE Wang. These key variables thus form a complex network of direct and indirect interactions.

Figure 7 further illustrates the marginal contributions of the four most critical variables to the predictions of ET and NEE under specific conditions. The 8 subplots reveal that, with the exception of STL2 for NEE which shows no distinct threshold, all other variables exhibit contribution thresholds in their predictions of the target variables. The threshold for SR's contribution to rainforest ET prediction is approximately 400 W m⁻². Below this value, its contribution is negative, turning positive above it. In contrast, the threshold for SR's contribution to NEE prediction is around 200 W m⁻², which aligns with the light saturation point commonly observed in plant photosynthesis Zhang. Beyond this point, SR's marginal contribution to NEE becomes negative, highlighting a unique aspect of the light-response mechanism in this tropical rainforest ecosystem. The patterns of marginal contribution for ATATC to both ET and NEE are similar, with a threshold near 25°C. Contributions are positive below and negative above this temperature. The marginal contributions of NGVP to ET and NEE predictions show opposite trends, with a threshold around 2.5 kPa. When NGVP > 2.5 kPa, its contribution to ET prediction is positive but negative to NEE; the opposite is true when NGVP < 2.5 kPa. For WSATC's contribution to NEE prediction, the threshold lies within an interval of 0.25 to 0.75 m s⁻¹. Contributions are generally positive when WSATC < 0.25 m s⁻¹ and negative when WSATC > 0.75 m s⁻¹. While a threshold for STL2's contribution to NEE is indicated at 25°C, the threshold effect is extremely subtle. It is noteworthy that the contribution thresholds of these variables are largely consistent for both ET and NEE predictions. This finding serves as key evidence for the coupled nature of carbon and water cycles in rainforests Naoya. The changes in ET and NEE are likely driven by a shared set of fundamental biophysical processes that operate within common optimal environmental ranges and share critical stress points. This observation reinforces the tight coupling of photosynthesis and transpiration mediated by stomatal behavior [Wang0,51]; Furthermore, these shared thresholds may represent critical tipping points where the functional state of the rainforest ecosystem undergoes transition. If future climate conditions persistently exceed these thresholds, it could simultaneously trigger a decline in carbon sink strength and an exacerbation of water loss through evapotranspiration, underscoring the importance of studying rainforest carbon-water coupling mechanisms in a changing environment.

Analysis of feature interactions indicates that the interaction between SR and ATATC is the strongest. High SR attenuates the positive contribution of high ATATC, while low ATATC mitigates the negative contribution of low SR (Figure 7a, e), interactions between other variable pairs are less pronounced. Nevertheless, the partial dependence plots align with established physical relationships, showing that high SR is associated with high ATATC (Figure 7a, e), and high ATATC coincides with high VP (Figure 7b, d, f).

Climatic conditions vary across months of the year. To investigate the seasonal dynamics in the environmental responses of ET and NEE, Figure 8 presents the monthly marginal contributions of each variable to their predictions. The marginal contributions exhibit pronounced seasonal variations. The marginal contributions of ATATC, AHATC, and NGVP to both ET and NEE show broadly similar seasonal patterns. The contribution of ATATC is lower in winter than in summer, peaking for ET in July and for NEE in September. The contribution of AHATC is greater during February, March, and April. NGVP shows higher contributions in winter and summer compared to spring and autumn. Notably, the marginal contribution of SR to ET prediction peaks in March, September, and October. These months coincide with the critical dry-to-wet and wet-to-dry transition periods in the Xishuangbanna rainforest climate Zafar, characterized by rapid rebalancing of energy, water, and plant physiological activity. In contrast, the influence of SR on NEE is greater in January, February, July, September, and November. This indicates that, despite the overall strong coupling, seasonal carbon-water decoupling occurs in the Xishuangbanna rainforest. The marginal contribution of STL2 to NEE prediction is greater from June to October, which aligns closely with the local rainy season (May to October). This provides important evidence that soil microbial respiration becomes a dominant factor in the carbon cycle under the warm and humid conditions prevalent during the rainy season. The marginal contributions of WS and PR to both ET and NEE predictions show relatively minor seasonal variation.

4. Discussion

4.1. Pathway Analysis of Environmental Impac(SEM)

While machine learning is effective in identifying key environmental variables influencing ET and NEE and capturing complex nonlinear relationships among them, it has limitations in establishing explicit causal links between external energy sources, atmospheric properties, soil characteristics, and the fluxes Sun. SEM addresses this limitation by quantifying the direct and indirect effects of environmental variables on ET and NEE. Figure 9 illustrates the pathways and standardized strengths of the effects of SR, AHATC, ATATC, ST, NGVP, and WS on ET and NEE. The analysis reveals that the energy provided by SR exerts a predominantly direct effect on NEE, with a path coefficient of -0.5034***, Other influential pathways on NEE are primarily indirect, mediated either through temperature or via the sequence from ATATC and AHATC to NGVP. The indirect pathway mediated through soil is relatively weaker, with a coefficient of only -0.4094***. For ET within the water cycle, the dominant influence occurs through the indirect pathway from ATATC and AHATC to NGVP. The direct effect of SR on ET is comparatively smaller, with a coefficient of 0.2043***. The feedback from evapotranspiration to environmental variables is most substantial in regulating ATATC, showing a strong negative coefficient of -1.4774***, The feedback to AHATC is weaker (-0.3156***), and the feedback to ST is minimal (0.0240***). The effects of ET and WSATC on NEE are relatively minor compared to other pathways, with coefficients of -0.1016*** and -0.1094*** respectively, The influence of WSATC on ET is of a similar magnitude (0.1069***). Overall, ATATC emerges as the most critical environmental variable through which SR indirectly regulates both ET and NEE.

4.2. Reliability of the Explainable Machine Learning Framework: Bridging Data Patterns and Biophysical Mechanisms

A key insight from this study is that the explainable machine learning (ML) framework successfully captured the complex response patterns of tropical rainforest carbon-water fluxes to environmental drivers—patterns that are intricately coupled with seasonal phenology. For instance, the model identified that the period of peak solar radiation influence on evapotranspiration coincided precisely with the transitional phases between dry and wet seasons. Similarly, the enhanced influence of soil temperature during the rainy season aligned with the period when soil microbial activity typically dominates carbon cycling dynamics. The high degree of consistency between the model-derived patterns and established, independent ecological knowledge demonstrates that the ML approach is not an impenetrable “black box.” Rather, it serves as a powerful tool for effectively revealing and quantifying known ecophysiological processes. More importantly, by using ML to identify key variables and their effects, and subsequently integrating these insights with confirmatory methods like Structural Equation Modeling (SEM), we were able to further disentangle the intrinsic quantitative relationships and causal pathways underlying these processes. Therefore, this work not only deepens the understanding of ecosystem functional dynamics in the Xishuangbanna tropical rainforest but also provides a robust paradigm for applying and validating data-driven methodologies in the study of complex ecosystems. The integrated framework presented here offers a replicable approach to generate mechanistically interpretable insights from high-dimensional observational data, bridging the gap between pattern detection and causal understanding.

4.3. Forecasting of ET and NEE Under Different Carbon Emission Scenarios in the Future

Utilizing the selected key variables, we retrained the optimal machine learning model and combined it with the monthly CMIP6 future key variable data corrected via quantile mapping to predict future ET and NEE. Figure 10(a) shows the historical actual values and CMIP6 simulated values for each variable. Using these values, we constructed mapping functions to correct the CMIP6 simulations. Compared to the uncorrected CMIP6 simulations, the corrected CMIP6 values show a higher degree of overlap with the actual observed curves, with Pearson correlation coefficients reaching: ET-0.5037, NEE-0.2374, AHATC-0.7349, SR-0.6689, WSATC-0.4100, ATATC-0.9527, STL2-0.2609, and NGVP-0.9597. Key input variables AHATC, SR, ATATC, and NGVP all achieved good correction results (with NGVP calculated based on the Magnus formula). We used the mapping functions to correct all variable simulations from CMIP6 under four future climate scenarios. Based on the corrected input variable data and the trained XGBoost model, we output the machine learning predictions for future ET and NEE. Figure 10(b) shows the changes in machine learning-predicted and CMIP6-corrected ET and NEE from 2026 to 2100.

It can be observed that, based on CMIP6 data, the Xishuangbanna rainforest will complete the transition from a carbon sink to a carbon source (NEE turning from negative to positive) between 2046 and 2066, depending on the future climate scenario; the higher the carbon emissions, the earlier this transition occurs. This result is consistent with predictions from previous studies Katharyn. However, the XGBoost predictions indicate that under all four scenarios, the carbon sink capacity of the Xishuangbanna rainforest will gradually decline (negative NEE increases), but it will not transition to a carbon source. This outcome aligns more closely with the flux observation trends over the past 20 years (increasing carbon sink) Liu, likely because the strict conservation policies implemented over the past few decades have introduced positive effects into the data Wang. It also suggests that, under the current level of protection, the Xishuangbanna rainforest ecosystem will continue to stably maintain its carbon sink function until the end of this century. For ET, the XGBoost predictions show a similar trend to the CMIP6 results, though CMIP6 values are slightly higher than those from XGBoost.

Overall, for human-intervened ecosystems like the Xishuangbanna rainforest, which have been under long-term strict protection, their future trajectory may not simply follow the path predicted by purely climate-driven models. Proactive and sustainable local management actions should be considered a critical, non-negligible variable in determining the fate of global ecosystems. Future research needs to further quantify the boundaries and mechanisms of this protective effect and more effectively integrate it into Earth System Models to enhance the completeness of global change predictions and the precision of policy responses.

5. Conclusions

This study elucidates the complex environmental regulation of carbon and water fluxes in the Xishuangbanna tropical rainforest, a unique ecosystem under long-term strict protection, and projects their future trajectories of ET and NEE under various climate scenarios. Our integrative analysis reveals three layers of understanding:

Firstly, at the mechanistic level, we have decoded the core biophysical drivers. Solar radiation (SR) serves as the paramount energy source but regulates carbon uptake (NEE) and water loss (ET) through distinctly different pathways: it directly controls NEE, while indirectly mediating ET via its influence on atmospheric variables. This dual role explains the observed tight coupling yet seasonal decoupling between the forest's carbon and water cycles. Building upon this mechanistic understanding, we then identified the critical environmental thresholds that govern overall ecosystem function, such as for air temperature (~25°C) and vapor pressure (~2.5 kPa). The quantification of these thresholds is essential for diagnosing vulnerability and predicting systemic stability under future climate regimes.

Most significantly, integrated machine learning and bias-corrected CMIP6 framework provides a novel resilience assessment. It reveals a fundamental divergence: traditional climate models project a mid-century sink-to-source transition, while our ML model, trained on conservation-era data, projects a persistent carbon sink. This indicates that sustained protection may buffer ecosystems against a purely climate-driven decline.