Exploring the Temperature Response of Nighttime Ecosystem Respiration in a Deciduous Broadleaf Forest Using Long-Term Observations and Causal Shapley Value

Abstract

Robust estimation of the temperature sensitivity of ecosystem respiration (ER) is critical for projecting terrestrial carbon dynamics under climate change. While ER is traditionally modeled as an exponential function of temperature, recent evidence suggests a unimodal response, with respiration peaking at an optimum temperature (T_opt) and declining thereafter. However, confounding variables like soil moisture, vegetation status, and vapor pressure deficit complicate isolating the effect of temperature. In this study, we combined long-term eddy covariance data from the flux tower in Sorø, Denmark (DKSor) with an explainable machine learning framework to disentangle the direct effects of air temperature (T_air) on nighttime ER. Using an optimized machine learning model and Causal Shapley value, we quantified the marginal contributions of environmental predictors while accounting for known causal dependencies. Then we employed a new concept of T_opt, which is a temperature that maximizes the Causal Shapley value of T_air (T_opt^SHAP). Our results revealed a clear unimodal temperature response (T_opt^SHAP = 15.99℃), suggesting that ER may begin to decline at a cooler temperature than previously assumed (T_opt = 18.25℃). A binning analysis furtherrevealed that unimodalresponse could be observed in normalized difference vegetation index ranges where high temperature (≥15℃) occurred, suggesting that unimodal temperature responsewas not biased by vegetation phenology. These findings underscore the importance of incorporating causal inference and interpretable machine learning in capturing the temperature response of ER.

1. Introduction

Understanding how ecosystem respiration (ER) responds to rising temperature is essential for accurately assessing global terrestrial carbon sequestration. Terrestrial ecosystems play a critical role in regulating global carbon fluxes, absorbing approximately 30% of annual anthropogenic CO₂ emissions (Intergovernmental Panel on Climate Change, 2023). Several studies project that rising temperature could significantly enhance ER, potentially causing terrestrial ecosystems to transition from carbon sinks into carbon sources by the mid-21st century (Cox et al., 2000; Duffy et al., 2021; Green et al., 2019). Since projections of future carbon dynamics largely depend on the balance between ER and gross primary production (GPP), understanding the temperature responses of these two fluxes is essential (Yuan et al., 2025). While GPP’s temperature dependence is well described by biochemical models such as the Farquhar–von Caemmerer–Berry (FvCB) model (Bernacchi et al., 2001, 2013; Medlyn et al., 2002), no comparable mechanistic framework exists for ER (Atkin et al., 2017; Niu et al., 2024). This gap in understanding two major carbon fluxes introduces uncertainties in carbon flux projections, which can misinform global climate policies and ecosystem management strategies. Therefore, robust estimation of the temperature sensitivity of ER is critical for improving carbon cycle assessments and ensuring reliable predictions of terrestrial carbon uptake.

Respiration is traditionally modeled as an exponential function of temperature (Carey et al., 2016; Enquist et al., 2003; Gillooly et al., 2001; Lloyd and Taylor, 1994; Tjoelker et al., 2001), but recent studies suggest that respiration, ranging from soil microbes and plants to entire ecosystems, often exhibits a unimodal response to temperature (Liu et al., 2018; Michaletz and Garen, 2024; Scafaro et al., 2021). In this framework, respiration rate increases with temperature up to an optimum temperature (T_opt) where ER is maximized and subsequently declines. This unimodal response results from multiple interacting mechanisms, including thermal adaptation of organisms (Sun et al., 2023), enzyme inactivation at higher temperature (Fanin et al., 2022), and substrate limitation (Scafaro et al., 2021). Emerging evidence indicates that even ER, comprising both autotrophic and heterotrophic components, follows this unimodal temperature pattern. Chen et al. (2023) analyze ER response to air temperature (T_air) using the FLUXNET 2015 dataset across 212 global sites (Baldocchi, 2020; Pastorello et al., 2020).

They construct ER-T_air response curves by binning growing-season T_air and daily ER for each site-year and then averaging the values within each temperature bin, following the methodology established by Niu et al. (2012). The result indicates that 183 out of 212 sites exhibit a unimodal ER-T_air relationship, with site-specific T_opt values mostly influenced by maximum temperature and the T_opt of GPP. Likewise, Meng et al. (2024) analyze the half-hourly nighttime net ecosystem exchange (NEE_night) response to air nighttime T_air and soil temperature (T_soil) across 196 FLUXNET2015 sites (Baldocchi, 2020; Pastorello et al., 2020), using comparable methods. The study finds that 66 sites displayed unimodal temperature response for both T_air and T_soil. Furthermore, analysis demonstrates that relying on exponential rather than unimodal equations leads to significantly higher estimates of NEE_night compared to observed values. These findings suggest that ER may not increase as rapidly with warming as previous projections, indicating that current Earth system models (ESMs) that incorporate exponential temperature-ER equations might be overestimating terrestrial carbon emissions.

However, many studies struggle to disentangle the mixed effects of multiple factors influencing the temperature-ER relationship, thereby introducing biases (Chen et al., 2023; Meng et al., 2024). Specifically, Chen et al. (2023) average daily T_air and corresponding carbon flux values within temperature bins to construct ER–T_air response curves. To quantitatively evaluate the shape of these curves for each site-year, they fit a quadratic function, where the parameters T_opt and maximum of ER represent the vertex of the parabola. Subsequently, to assess potential confounding effects, the study employs a linear mixed-effects model with random effects at the site level, testing for significant linear relationships between each potential confounding variable and the residuals of the T_air–ER regression. The absence of significant relationships leads the researchers to conclude that the temperature response of ER is not influenced by these confounding factors. However, the method used to rule out confounding factors presents notable limitations. Residual analysis might fail to capture non-linear or complex interactions among environmental variables such as soil moisture (SM), vapor pressure deficit (VPD), solar radiation, and leaf area index (LAI). Such complex interactions, including threshold effects, might remain undetected through a simple linear regression of residuals (Becher, 1992; Pearl, 1998).

Additionally, multicollinearity among confounders, which frequently covary with each other and temperature, complicates the isolation of individual factors’ effect on ER. Specifically, when analyzing the residuals of the temperature-ER relationship, the influence of temperature on ER is removed. If another variable like LAI is highly correlated with temperature, its impact may also be masked in the residuals. This can lead to a situation where, despite LAI affecting ER, no significant correlation appears between LAI and the residuals of the temperature-ER relationship. Finally, the quadratic model representing the T_air–ER relationship may oversimplify ecological dynamics, suggesting the need for more sophisticated modeling techniques to effectively capture the actual interactions between temperature and other environmental factors.

Advanced machine learning models, such as XGBoost, combined with explainable artificial intelligence (XAI) techniques like Causal Shapley value, effectively overcome these limitations. XGBoost enables representation of complex, non-linear interactions between temperature and other confounding environmental variables (Grinsztajn et al., 2022). XGBoost handles multicollinearity more effectively than traditional linear models, as it does not rely on estimating linear coefficients. However, correlated predictors may still share or dilute feature importance, which can affect interpretability. Using XAI methods like Causal Shapley value helps explicitly address these interpretability issues by incorporating known causal structures. Shapley Additive Explanations (SHAP), a popular XAI method, quantify each predictor’s contribution to individual model predictions using Shapley value, derived from cooperative game theory (Shapley, 1953). SHAP excels in providing interpretable insights into high-dimensional and nonlinear model predictions (Lundberg and Lee, 2017), clarifying how each predictor (e.g., temperature) specifically influences outcomes (e.g., respiration rate) within a specific context, not just reporting overall feature importance.

While SHAP has been successfully applied to assess the impact of climate variables like temperature and precipitation on ecosystem productivity and carbon dynamics (Xie et al., 2024), SHAP assumes independence among predictors (Lundberg and Lee, 2017). Thus, it can not fully capture causal relationships and sometimes results in counterintuitive explanations (Janzing et al., 2020). To overcome this limitation, Causal Shapley value enhances this approach by incorporating causal reasoning frameworks, notably Pearl’s do-calculus (Pearl, 2012), explicitly accounting for causal relationships within analyses (Pearl, 2022; Runge et al., 2023). This integration allows researchers to isolate and quantify each variable’s direct causal effect transparently, addressing complexities arising from non-linear interactions and multicollinearity. Consequently, combining XGBoost and Causal Shapley value presents a robust alternative to simpler quadratic models, effectively capturing the complexity of ecological dynamics and yielding more accurate and interpretable insights into the temperature response of ER.

In addition to the Causal Shapley value, the binning method provides another effective approach for analyzing complex relationships among variables while controlling for confounding factors. For example, when investigating acclimation of photosynthesis to light intensity, Jiang et al. (2020) construct T_air–fraction of Absorbed Photosynthetically Active Radiation (fAPAR) bins to account for the confounding effects of temperature and foliage density, which covary with light intensity on intermediate timescales. Similarly, Liu et al. (2020) use the binning method to disentangle the high correlation between SM and VPD. More recently, Liu et al. (2024) observed thermal acclimation of photosynthetic rate along gradients of T_air–fAPAR bin pairs, effectively controlling for the simultaneous changes in temperature and foliage quantity. Ultimately, when combined with XAI techniques, such a refined analytical approach provides a robust framework for assessing the unimodal temperature response of ER by rigorously controlling for confounding variables.

This integrated approach, combining XAI techniques and strategic binning, enables us to clearly distinguish the direct causal influence of temperature on nighttime ER (ER_night) from indirect effects and spurious correlations. The primary goals of this research are to (1) rigorously determine whether ER_night exhibits a unimodal temperature response after adequately controlling for confounding variables such as phenological changes and water availability; (2) enhance the understanding of temperature–ER relationships to improve predictions of terrestrial carbon flux dynamics under future climate scenarios, thereby supporting more reliable climate change projections. Unlike previous studies that rely on empirical binning or statistical methods, we introduce a robust causal inference framework using Causal Shapley value analysis. This approach enables the quantification of temperature’s direct effect on respiration, separating it from indirect or spurious influences of correlated variables. By integrating causal reasoning with explainable machine learning, this study provides a new pathway for deriving interpretable and robust estimates of T_opt, offering improved insights into the nonlinear behavior of ER under climate change.

2. Materials and Methods

2.1. Data Collection and Preprocessing

The overall workflow of this analysis is presented in Fig. 1. This study utilized eddy covariance measurements from the FLUXNET 2015 dataset (1996–2014; Baldocchi, 2020; Pastorello et al., 2020), with a primary focus on the DK-Sor site (Pilegaard et al., 2011), a temperate deciduous broadleaf forest (DBF). Among all available sites, DK-Sor provided the largest number of valid half-hourly observations spanning 19 years from 1996 to 2014. The site is located in Sorø, Denmark (55.4859°N, 11.6446°E) at 40 m above sea level with a mean annual temperature of 8.2°C and mean annual precipitation of 660 mm.

Fig. 1. An overview of the analysis workflow. The process begins with the collection and preprocessing of FLUXNET2015 nighttime NEE data and MODIS-derived NDVI. Filtering steps retain nighttime, measured, and growing-season data. Following correlation analysis, an XGBoost model is trained using randomized hyperparameter optimization. Model interpretation is conducted via Causal Shapley value analysis, including PCA-based representative sampling, response curve generation, and binning to derive temperature response patterns.

To estimate ER_night, we employed directly measured growing season NEE_night data (NEE_VUT_REF_QC = 0). Additionally, only positive NEE values were retained under the assumption that nighttime CO₂ fluxes represent respiratory flux from the ecosystem. The growing season was defined as the set of months with a mean monthly T_air exceeding 0°C (Prentice et al., 2011).

To represent vegetation dynamics, we calculated the normalized difference vegetation index (NDVI) using the Moderate Resolution Imaging Spectroradiometer (MODIS) MCD43A4.061 Nadir bidirectional reflectance distribution function (BRDF)-adjusted reflectance 500 m daily product. NDVI was computed as Eq. (1):

\(\begin{align}N D V I=\frac{N I R-R e d}{N I R+R e d}\end{align}\)        (1)

where NIR and Red represent the near-infrared and red reflectance bands, respectively. To ensure consistency in satellite-derived vegetation indices, we excluded all samples collected prior to the operational period of MODIS (i.e., before February 24, 2000).

2.2. Model Development and Hyperparameter Optimization

To identify environmental variables significantly associated with ER_night, we conducted correlation analyses with environmental and meteorological variables in the FLUXNET2015 dataset (Fig. 2). Following these analyses, T_air, VPD, T_soil, soil water content (SWC), and NDVI were selected as predictor variables for the modeling phase based on their significant correlations (absolute correlation coefficient (|r|) > 0.2) with ER_night. As shown in Fig. 2, ER_night (represented by NEE_night) was found to be positively correlated with T_air (r = 0.68), T_soil (r = 0.69), and NDVI (r = 0.58), while moderate correlations were also observed with VPD (r = 0.35) and SWC (r = –0.23). Incoming longwave radiation (LW_in) and outgoing longwave radiation (LW_out) were not used because T_soil can represent both variables. Based on these results, five variables—T_air, VPD, T_soil, SWC, and NDVI—were selected, using a cutoff of |r| > 0.2. The final variable information used in the modelling and subsequent analysis is presented in Table 1.

Table 1. Summary of variables used in the Causal Shapley value analysis

Descriptions, units, temporal coverage, and data sources for the response variable (NEEnight) and five predictors were obtained from the FLUXNET2015 dataset, except NDVI from the MODIS MCD43A4 Collection 6.1 product.

Fig. 2. Lower-triangle Pearson correlation matrix among environmental and meteorological variables and nighttime ecosystem respiration (ER_night) at the DK-Sor site. Warmer colors indicate stronger positive correlations, cooler colors negative ones. T_air, T_soil (1 cm), and NDVI showed strong positive correlations with ER_night, while VPD and SWC were moderately correlated. Radiation variables (SW_in, SW_out, LW_in, LW_out, R_net) and other meteorological factors (RH, U*, WS, P_atm, P, CO₂, PPFD_in, PPFD_out) were generally weakly associated. Variables with | r | > 0.2 were retained for modeling. Shortwave radiation (SW_in), outgoing shortwave radiation (SW_out), incoming and outgoing longwave radiation (LW_in and LW_out), net radiation (R_net), relative humidity (RH), friction velocity (U*), wind speed (WS), atmospheric pressure (P_atm), precipitation (P), atmospheric CO₂ concentration (CO₂), incoming and outgoing photosynthetic photon flux density (PPFD_in and PPFD_out).

To train the XGBoost model, we randomly partitioned the dataset into a training set (80%) and a test set (20%). To optimize model performance, a randomized hyperparameter search was implemented using the caret package in R. A total of 100 candidate hyperparameter combinations were randomly sampled from predefined search ranges: the number of boosting rounds (nrounds: 100–500), maximum tree depth (max_depth: 1–10), learning rate (eta: 0.01–0.30, log-uniformly sampled), row subsampling ratio (subsample: 0.60–1.00) and column subsampling ratio (colsample_bytree: 0.60–1.00). Each configuration was evaluated using five-fold cross-validation on the training set, with model selection based on the highest average coefficient of determination (R²) across folds. The best-performing model configuration was retained and subsequently used for interpretation via Causal Shapley value analysis.

To assess the model’s predictive consistency across environmental conditions, we evaluated relative prediction bias (%) within five equally sized bins of T_air and NDVI. Specifically, the full range of each variable was divided into five bins of equal width. For each bin, we calculated the mean predicted ER_night and the mean observed ER_night, and then computed the relative prediction bias as Eq. (2):

\(\begin{align}\text {Relative Bias}(\%)=\frac{\text { Mean Predicted }- \text { Mean Observed }}{\text { Mean Observed }} \times 100\\\end{align}\)       (2)

This calculation was performed independently for the training and test sets. The resulting bias values allowed us to evaluate the stability and generalization of model predictions across gradients of temperature and vegetation conditions.

2.3. Causal Shapley Value Analysis

To reduce computational cost in the Causal Shapley value analysis while preserving the diversity of predictor combinations, a representative subset of the whole dataset was selected. After centering and scaling, we applied a dimensionality reduction step to the full predictor matrix, comprising T_air, VPD, T_soil, SWC, and NDVI. Principal component analysis (PCA) was performed, and the first three principal components were retained to represent the major variation in the feature space. An initial set of 300 random observations was used to define the convex hull of the full predictor distribution in PC space. Then, a two-tiered resampling strategy was employed over 20 iterations. In each iteration, 2,000 observations were randomly drawn and used to construct a convex hull. The “coverage score” was computed as the fraction of the whole dataset lying within this hull.

Because the first three principal components captured the majority of variance in the original feature space (see Results), the convex hull in PC space served as an effective low-dimensional proxy for the multivariate distribution of the data. Although PCA introduces a linear transformation, the relative positions and combinations of original predictors are preserved along the axes of maximum variation. Thus, maximizing coverage in the PC space ensured that the selected subset retained broad representativeness across the joint distribution of the original features, while also enabling tractable computation for Causal Shapley value estimation. The sample that achieved the highest coverage was selected, and its corresponding original (non-PCA-transformed) predictor values were retained. This subset of 2,000 representative observations was used for all subsequent Causal Shapley value estimations.

Causal Shapley values were then computed for 2,000 representative samples to quantify the marginal contributions of each predictor to the model’s output while accounting for known causal relationships among variables. A domain-informed causal hierarchy was employed based on prior ecological and physiological knowledge (Fig. 3). T_air was treated as the primary driver, with VPD, T_soil, and SWC as intermediate variables influenced by temperature. NDVI was considered a tertiary variable, reflecting vegetation status influenced by environmental conditions.

Fig. 3. Directed Acyclic Graph (DAG) representing hypothesized causal relationships among the five environmental predictors used in the XGBoost model. Air temperature (T_air) is treated as the primary driver, influencing vapor pressure deficit (VPD), soil temperature (T_soil), soil water content (SWC), and Normalized Difference Vegetation Index (NDVI). NDVI is assumed to respond to integrated environmental conditions. Arrows denote assumed direct causal effects informed by ecological and physiological reasoning.

To begin Causal Shapley value analysis, a background dataset of 2,000 observations was randomly sampled from the filtered predictor matrix. This background served dual purposes: it provided a reference distribution for estimating expected model outputs and was used to compute the model’s baseline prediction (φ₀), defined as the mean prediction over the background set. All subsequent Causal Shapley values were interpreted as deviations from this baseline. Using the optimized XGBoost model and the selected representative sample, Causal Shapley values were estimated using the shapr package. To further investigate the response of ER to temperature, we fitted the T_air-ER_night curve and the T_air-Causal Shapley value curve for T_air using a generalized additive model (GAM), visualized with a density scatter plot. We also computed 95% pointwise confidence intervals, and identified the local maximum by detecting a sign change in the first derivative of the predicted curve. We interpreted this local maximum as the optimum temperature that maximizes the Causal Shapley value (T^SHAP_opt), similar to T_opt. It represents a potential ecological threshold at which temperature has the strongest influence on ER, beyond which its effect begins to decline.

To visualize interactions between T_air and its marginal contributions to predictions, we constructed a two-dimensional grid with NDVI and T_air as the two axes, and the mean Causal Shapley values of T_air were computed within each NDVI–T_air bin. This approach enabled the examination of how the contribution of T_air to ERnight varied across different levels of vegetation phenology and temperature. Additionally, NDVI was split into two representative vegetation states (0.7–0.8 and ≥0.8), and GAMs were fitted separately to the T_air-Causal Shapley value response curve within each NDVI range. Local maxima of the fitted curves were identified to verify whether they exhibit a unimodal response under different phenological conditions.

3. Results

3.1. Model Performance

The optimized XGBoost model exhibited strong predictive capability for ER_night, achieving a train R² of 0.94 and test R² of 0.73 (Fig. 4). This indicates that the model effectively captured complex, nonlinear relationships between environmental drivers and ER_night. Across both T_air and NDVI gradients, relative prediction bias remained within ±1% in most bins (Fig. 5). Bias was close to zero throughout the T_air range in the training set (Fig. 5a), while the test set showed slightly larger underestimation at the lowest temperature bin (Fig. 5b). NDVI-based patterns exhibited minimal variation, with bias values consistently low and balanced across all seasonal conditions (Figs. 5c–d). These results indicate stable and unbiased thermal and seasonal patterns in model performance.

Fig. 4. Observed versus predicted nighttime ecosystem respiration (ER_night) for the DK-Sor site using the optimized XGBoost model. (a) Training set (N = 37,514) and (b) test set (N = 9,377) predictions. Each dot represents a half-hourly observation, with color indicating the magnitude of observed ER_night. The red dashed line denotes the 1:1 line. Model performance was high on the training set (R² = 0.94) and remained robust on the independent test set (R² = 0.73), indicating that the model generalizes well to unseen data while capturing the nonlinear relationships between ER_night and environmental predictors.

Fig. 5. Relative prediction bias (%) of the XGBoost model across binned environmental conditions. (a–b) Bias by air temperature (T_air) bins for the (a) training and (b) test datasets. (c–d) Bias by NDVI bins for the (c) training and (d) test datasets.

3.2. Representative Selection and Mean Causal Shapley Values

Across 20 iterations of random sampling, the mean coverage score reached 97.5%, with a minimum and maximum of 97.0% and 97.9%, respectively. These values confirm that the selected subset of 2,000 samples was highly representative of the full data distribution in multivariate space. The average absolute Causal Shapley values from the selected subset revealed that T_air had the greatest impact on model predictions. NDVI ranked second, followed by T_soil and SWC. In contrast, VPD contributed minimally to ER_night (Fig. 6).

Fig. 6. Mean absolute Causal Shapley values for the five environmental predictors of nighttime ecosystem respiration (ER_night). Air temperature (T_air) is the most important variable, followed by normalized difference vegetation index (NDVI), soil temperature (T_soil), soil water content (SWC), and vapor pressure deficit (VPD).

3.3. Temperature Response Patterns and T_opt Estimation

When using the full dataset (n = 46,891), the T_air-ER_night curve fitted by GAM exhibited clear unimodal relationship with a T_opt observed at 18.25°C (Fig. 7a). However, when we examined the same relationship using the 2,000-point representative subset, the response curve lost its unimodal shape and instead showed a monotonically increasing trend (Fig. 7b). To investigate whether feature attribution could recover the original structure, we applied GAM to the normal Shapley values of T_air derived from the same subset. The resulting curve remained monotonic and did not exhibit a distinct peak, suggesting that standard Shapley values may not fully isolate the intrinsic effect of air temperature (Fig. 7c). When GAM was fitted to the Causal Shapley values of T_air, the curve recovered a clear unimodal pattern, with T^SHAP_opt at 15.99°C (Fig. 7d).

Fig. 7. Generalized additive model (GAM) fitted response curve between air temperature (T_air) and nighttime ecosystem respiration (ER_night) with density scatter plots. (a) A GAM fitted to the full dataset (n = 46,891), revealing a clear unimodal response between T_air and ER_night with an optimum temperature (T_opt) of 18.25°C. (b) A GAM fitted to a representative subsample of 2,000 points loses this unimodality, showing a monotonically increasing trend. (c) A GAM fitted to standard Shapley values derived from the same 2,000-point subset does not recover the unimodal shape. (d) A GAM fitted to the Causal Shapley values of T_air, recovering a unimodal response with a lower optimum temperature (T^SHAP_opt) of 15.99°C. All GAM fits include 95% confidence intervals as shaded bands.

Across the NDVI–T_air bins, the causal influence of T_air on ER_night exhibited a consistent unimodal pattern. In all NDVI bins where T_air extended into a higher range (≥15°C), the mean Causal Shapley value of T_air increased with temperature up to a clear local maximum (black dot in Fig. 8a) and then declined thereafter (Fig. 8a). To further investigate the robustness of the unimodal temperature response across vegetation states, we examined the T_air –Causal Shapley value response curve within two NDVI ranges: 0.7–0.8 and ≥0.8 (Figs. 8b–c). In both NDVI divisions, the fitted GAM curves displayed a clear unimodal pattern with distinct T^SHAP_opt, supporting the conclusion that the estimated T^SHAP_opt is not an artifact of vegetation variation but rather reflects a native response of ER_night to T_air.

Fig. 8. Binned heatmap illustrating the variation in the Causal Shapley value of air temperature (T_air) on nighttime ecosystem respiration (ER_night) across normalized difference vegetation index (NDVI) and T_air gradients. (a) Binned heatmap illustrating the variation in the Causal Shapley value of T_air on ER_night across NDVI and T_air gradients. Each cell represents the average Causal Shapley value of T_air within an NDVI–T_air bin, and only bins with at least five observations are outlined in black. Grids with a local maximum are marked with black dots for NDVI values ≥ 0.65, highlighting the emergence of a temperature optimum (T^SHAP_opt). (b–c) Generalized Additive Model fits of the T_air-Causal Shapley value response curve, stratified by NDVI bin. Panel (b) shows the response in NDVI 0.7–0.8, and panel (c) shows the response under NDVI ≥ 0.8. In both cases, a unimodal pattern is evident, and the estimated T^SHAP_opt is indicated by a red dot.

4. Discussion

4.1. Effectiveness of Causal Shapley Values in Examining Unimodal Temperature Response

In this study, we introduced a framework for estimating the T^SHAP_opt of ER_night, defined as the temperature that maximizes Causal Shapley value of T_air in estimating ER_night. This contrasts with previous approaches that estimate temperature sensitivity based on the conventional statistical T_air-ER_night relationship, which is susceptible to confounding effects (Chen et al., 2023; Meng et al., 2024). Causal Shapley values have shown their effectiveness in isolating the direct effect of the independent variable on the dependent variable, even in the presence of correlated or confounding variables. For instance, Liu et al. (2025) successfully applied this approach to disentangle the separate effects of VPD and SM on photosynthesis, despite their strong interdependence. Similarly, while the T_air-ER_night relationship exhibited a unimodal pattern when using the full dataset, this pattern disappeared in a representative subsample of 2,000 observations, where the response appeared monotonically increasing. Interestingly, this monotonic increase persisted even when using normal Shapley values, which reflect feature importance under conditional expectations.

In contrast, the T_air-Causal Shapley value curve clearly showed a unimodal response despite being applied to the same subset. This suggests that the loss of unimodality cannot be fully attributed to the sampling noise, and causality contributes meaningfully to the temperature response relationship. This highlights the robustness of T^SHAP_opt in examining the intrinsic temperature response, not only by addressing confounding factors through causal reasoning but also by maintaining stability under data-limited conditions. While T^SHAP_opt is not intended to represent a physiological optimum directly, it serves as a data-driven indicator of the temperature at which the marginal causal effect on respiration is highest, offering a complementary perspective to traditional empirical methods.

4.2. Implications of a Lower T^SHAP_opt Estimate

The estimated T^SHAP_opt was approximately 15.99°C, notably lower than the 18.25°C derived from the traditional regression method (Meng et al., 2024). This discrepancy may result from the influence of uneven data distribution or confounding variables such as foliage quantity that elevate ER_night and covary with temperature, leading to an overestimation of T_opt. If such early saturation of temperature response is consistent across sites, this suggests that the temperature-driven increase in respiration may level off or even reverse sooner than previously assumed (Chen et al., 2023; Cox et al., 2000; Duffy et al., 2021; Green et al., 2019). This finding has important implications for projections of carbon–climate feedbacks. If terrestrial respiration is less sensitive to warming than current models assume, it may result in lower carbon emissions under future warming scenarios, partially mitigating expected climate impacts.

4.3. Limitations and Future Directions

While our study incorporated a causal structure and XAI on the temperature sensitivity of ER, several limitations should be acknowledged. First, we focused exclusively on ER_night to avoid uncertainties from daytime flux partitioning. Although this improves data reliability, it leaves out potential daytime-specific temperature responses. Future work should integrate both daytime and nighttime ER components to develop a more comprehensive understanding of total ER dynamics.

Second, our analysis was confined to a single temperate deciduous forest site (DK-Sor). Although this site offers a long-term, high-quality dataset, extending this causal framework to diverse biomes and climate zones will be essential to assess the generalizability of our findings. In addition, future work should perform inter-site comparisons and sensitivity analyses to evaluate the robustness of T^SHAP_opt across ecosystem types and modeling choices.

Third, while the GAM response curves include confidence intervals, the derived T^SHAP_opt does not have an associated uncertainty range due to the lack of a standardized approach for quantifying uncertainty in Causal Shapley values. Developing such methods is an important next step for enabling robust model parameterization. And to enhance the credibility of causal inference, future research should consider applying do-intervention simulations or sensitivity analyses on the assumed directed acyclic graph structure. Such approaches would enable formal evaluation of whether the estimated Causal Shapley values remain robust under varying structural assumptions or potential violations of the causal graph.

Finally, although the Causal Shapley value provides a robust framework to measure the direct contributions of individual variables, it does not account for the physiological processes that may underlie these responses, such as microbial activity, substrate availability, or enzyme regulation. As a result, the method identifies correlations between variables, but does not explain the underlying processes. Integrating this approach with process-based models or physiological function could help connect observed patterns to underlying mechanisms and strengthen both the interpretability and predictive capacity of ecosystem models. Also, while the causal T^SHAP_opt identified in this study is derived from in-situ eddy covariance observations, its correspondence with leaf or tissue level physiological thresholds remains to be validated. Future studies that directly compare T^SHAP_opt with experimental saturation points would help bridge this gap.

5. Conclusions

This study suggests a robust framework to examine the unimodal temperature response of ER_night in a temperate forest using long-term flux data and a causal machine learning approach. By isolating the direct effect of T_air, we found that respiration peaked at around 15.99°C, suggesting that respiration may begin to decline at lower temperatures than previously assumed. This pattern persisted across varying levels of vegetation activity, supporting the consistency of thermal sensitivity. Our findings highlight the value of causal explainable models in improving predictions of ecosystem carbon dynamics under climate change.

Author Contributions

Conceptualization: Son MH, Jeong SC; Data curation: Son MH; Methodology, Formal analysis, Validation: All authors; Project administration: Jeong SC; Supervision: Jeong SC, Zhang H; Writing–original draft: Son MH; Writing–review & editing: All authors.

Conflicts of Interest

No potential conflict of interest relevant to this article was reported.

Funding

This research was supported by the Technology Development Project for Creation and Management of Ecosystem-based Carbon Sinks (202300218237) through the Korea Environmental Industry & Technology Institute (KEITI), funded by the Ministry of Environment (MOE).

Data Availability Statement

The dataset used to calculate NDVI is available on the Google Earth Engine platform: https://developers.google.com/earth-engine/datasets/catalog/MODIS_061_MCD43A4?hl=ko#bands. Flux and environment data used in this study can be found on the FLUXNET2015 website: https://fluxnet.org/data/fluxnet2015-dataset/.