1 Introduction

Forests are critical components of the global carbon cycle, acting as both sinks and sources of atmospheric CO\(_2\) through their carbon pools (e.g., above-ground biomass, soil organic carbon) and fluxes (e.g., respiration, net primary productivity). Disturbances–whether natural or anthropogenic, including wildfire, windthrow, pest outbreaks, or harvesting–can substantially alter ecosystem carbon dynamics, with effects that depend on both disturbance intensity and environmental context. With increasing frequency and severity of disturbances under climate change, understanding potential non-linear responses of disturbance severity and environmental moderators of recovery is essential for informing forest management and climate mitigation strategies [1,2,3].

Large-scale empirical experiments provide crucial insight into these dynamics. The Forest Resilience Threshold Experiment (FoRTE), a manipulative forest study in the Upper Great Lakes region, offers open-access measurements of carbon pools and fluxes across replicated plots subjected to varying disturbance severities [4]. FoRTE measurements include soil CO\(_2\) flux, leaf-level gas exchange, above-ground biomass, coarse woody debris, and environmental covariates such as soil temperature, volumetric water content, and light availability. Complementary modeling studies using the Ecosystem Demography model version 2.2 (ED−2.2), a cohort-based dynamic vegetation simulator, suggest that high disturbance severity can reduce immediate flux resistance while enhancing recovery, highlighting trade-offs between short-term vulnerability and long-term ecosystem resilience [5]. Similarly, moderate logging disturbances may alter species composition and turnover rates without necessarily reducing diversity, with implications for productivity and nutrient cycling [6]. Broader analyses across China, Africa, and Mediterranean regions underscore that land-use change, disturbance regime, and intercropping systems influence carbon budgets in context-specific ways [7,8,9].

Disturbance and recovery dynamics extend beyond biophysical processes into socio-ecological contexts. Post-disaster studies in East Asia and Australia illustrate that government intervention, institutional capacity, and socio-economic factors can shape recovery trajectories, emphasizing the importance of contextual moderators in resilience studies [10, 11]. These insights highlight the need for models that integrate ecological heterogeneity, disturbance gradients, and environmental controls.

Threshold detection frameworks offer powerful tools for analyzing non-linear responses in noisy and heterogeneous ecological data. Hierarchical Bayesian change-point regression has been successfully applied to identify ecological thresholds in temperate forests [12], aquatic ecosystems [13], and other longitudinal or spatially structured studies [14,15,16,17,18,19,20,21]. Applied ecological studies using thresholds include modeling grouse nesting along disturbance gradients [22], hydrological disturbance thresholds in British Columbia watersheds [23], and satellite-based vegetation recovery in US drylands [24]. Advances in covariate selection, predictive modeling, and uncertainty quantification provide additional methodological guidance [25,26,27], while simulation frameworks in carbon capture and biomass supply chains demonstrate the value of integrating statistical and computational approaches [28, 29].

Long-term forest studies indicate that disturbance legacy effects can persist for decades, with full recovery to net carbon sinks taking one to two decades in North American temperate forests [30], and demographic structure strongly influencing ecosystem productivity [31]. Intensifying disturbance regimes may buffer biodiversity losses in some contexts but complicate static conservation planning [32, 33], creating the so-called “disturbance paradox,” whereby disturbances simultaneously reduce ecosystem services and enhance species diversity [34].

In summary, (i) disturbances fundamentally alter forest carbon dynamics at multiple scales, (ii) threshold effects are prevalent but challenging to quantify without robust statistical frameworks, and (iii) hierarchical or frequentist change-point models provide a principled approach to detect non-linear responses under environmental heterogeneity. Building on these insights, this study applies a frequentist threshold mixed-effects model to FoRTE data, integrating environmental covariates to link disturbance severity with recovery dynamics.

1.1 Contributions

This paper makes the following contributions:

  • We analyze the FoRTE dataset [4], applying threshold mixed-effects models to explore potential disturbance severity thresholds affecting soil CO\(_2\) fluxes and carbon pools.

  • Environmental covariates (soil moisture, temperature, light availability) are incorporated to examine differential recovery rates under varying conditions.

  • We provide interpretable, conditional insights for hypothesis generation, identifying disturbance severities associated with slower or faster recovery trajectories and highlighting environmental contexts that modulate these responses.

The remainder of the paper is structured as follows. Section 2 describes the dataset, Sect. 3 outlines the statistical methodology, Sect. 4 presents the results, Sect. 5 interprets the findings, and Sect. 6 concludes.

2 About the dataset

We analyze ecological observations from the Forest Resilience Threshold Experiment (FoRTE), a large-scale manipulative forest experiment established in northern temperate forests of the Upper Great Lakes region of the USA. FoRTE was designed to quantify how forest ecosystems respond to varying disturbance severities, representative of windthrow, insect outbreaks, and harvesting, under changing climatic conditions.

The experiment consists of replicated forest plots subjected to controlled disturbance treatments spanning a gradient from undisturbed reference conditions to high-severity canopy removal. Individual experimental plots are approximately 0.1 ha in size and are distributed across mixed temperate forest stands dominated by species such as Acer saccharum (sugar maple), Populus tremuloides (trembling aspen), and Betula papyrifera (paper birch), together with associated hardwood and conifer species characteristic of the Upper Great Lakes region.

Plot layout and tree sampling: Experimental plots are spatially separated within forest stands, typically by distances ranging from several tens to hundreds of meters depending on site configuration and treatment allocation. This spacing minimizes direct treatment interference while preserving comparable climatic and edaphic conditions among plots. Within each plot, tree inventories follow standard forest sampling protocols in which all trees exceeding a minimum diameter-at-breast-height threshold are measured. Recorded attributes include tree diameter, mortality status, and regeneration dynamics, allowing characterization of vegetation structure and disturbance impacts through time.

Temporal monitoring and measurements: Vegetation structure and ecosystem functioning are monitored through repeated post-disturbance measurement campaigns conducted at regular revisit intervals following treatment implementation. Observations include ecosystem carbon pools (e.g., above-ground biomass and coarse woody debris) and carbon flux measurements such as soil respiration and leaf-level gas exchange. Environmental drivers including soil temperature, volumetric water content, and light availability are simultaneously recorded using standardized field instrumentation across all experimental plots.

The spatial replication and repeated temporal measurements enable assessment of ecosystem resistance and recovery dynamics across disturbance gradients while accounting for natural environmental variability.

All analyses presented in this study use publicly available observations distributed through the fortedata R package, which provides harmonized plot-level measurements, disturbance treatment metadata, species composition summaries, and spatial identifiers derived from the FoRTE experimental network. In particular, the model fitting and validation analyses reported here are based on the FoRTE soil respiration dataset, together with the corresponding treatment and environmental variables used after preprocessing.

Details regarding data indexing, preprocessing procedures, validation design, and statistical implementation are provided in Sect. 3, where the hierarchical modeling framework is introduced.

3 Methodology

This study develops a hierarchical statistical framework to quantify how disturbance severity and environmental conditions regulate soil CO\(_2\) flux dynamics following forest disturbance. The modeling strategy integrates (i) dominant environmental controls, (ii) nonlinear ecosystem responses to disturbance intensity, and (iii) temporal recovery processes after disturbance.

The methodological formulation explicitly reflects the replicated spatial plot design and repeated temporal measurements of the FoRTE experiment described in Sect. 2. In particular, the hierarchical structure and validation strategy are designed to account for potential spatial dependence among plots while preserving statistically independent model evaluation.

3.1 Data structure and preprocessing

Let

$$ \mathcal {I}=\{1,\dots ,N\} $$

denote experimental plots spatially distributed across study sites (where \(N = 4\) in this study), and

$$ \mathcal {T}_i=\{1,\dots ,T_i\} $$

measurement occasions for plot i.

For plot i at time t:

  • \(Y_{i,t}>0\) represents observed soil CO\(_2\) flux,

  • \(d_i\in [0,1]\) denotes assigned disturbance severity,

  • \(\textbf{x}_{i,t}\) contains standardized environmental covariates,

  • \(t_{0,i}\) is disturbance timing,

  • \(s_{i,t}=\max (0,t-t_{0,i})\) denotes time since disturbance.

Observations affected by instrument malfunction were removed prior to analysis. Continuous predictors were centered and scaled to improve numerical stability and parameter identifiability. Flux observations containing near-zero values were stabilized using a small additive constant prior to logarithmic transformation.

Because measurements within the same plot are temporally repeated and nearby plots may share environmental similarities, model validation is conducted at the plot level. Entire plots are held out during evaluation using a leave-one-plot-out cross-validation scheme. Given that the experiment incorporates \(N = 4\) distinct plots, this validation process involved exactly 4 separate iterations. At each iteration, a single plot constitutes the independent hold-out validation dataset, while the remaining \(N-1\) (i.e., three) plots are used exclusively for model fitting.

Although the validation is performed on the same processed analytical dataset derived from the FoRTE soil respiration resource, the held-out plot is completely excluded from model fitting in each fold; thus, validation is strictly conducted on unseen observations at the plot level rather than on the training data itself. This design is appropriate here because the study is based on a single coherent FoRTE soil respiration dataset with repeated measurements across plots, and the principal goal of validation is to assess out-of-sample predictive stability across experimental units while avoiding within-plot information leakage.

Furthermore, validation using an external dataset was not performed because our modeling framework is tightly coupled with the specific, controlled manipulative disturbance gradient (0% to 85% severity) and the paired high-frequency environmental tracking unique to the FoRTE experimental layout. An independent external dataset mirroring this precise experimental manipulation, temporal frequency, and covariate structure is currently unavailable. Consequently, this reliance on internal cross-validation is acknowledged as a structural constraint regarding broader geographical generalizability, shifting the focus of our validation toward ensuring the internal stability and robustness of the inferred ecological relationships.

This strategy prevents spatial and temporal information leakage between training and testing data and provides an assessment of predictive performance for previously unseen experimental units. Additional leave-one-severity-band-out validation evaluates extrapolation across disturbance gradients.

3.2 Outcome distribution

Soil CO\(_2\) fluxes are strictly positive and typically right-skewed. We therefore assume a log-normal observation model:

$$\begin{aligned} Y_{i,t}\mid \eta _{i,t},\sigma ^2 \sim \textrm{LogNormal}(\eta _{i,t},\sigma ^2), \end{aligned}$$
(1)

where

$$ \eta _{i,t}=\log \mu _{i,t}. $$

This formulation represents multiplicative ecological variability and allows regression coefficients to be interpreted as proportional changes in expected flux.

3.3 Environmental drivers

Soil respiration responds to interacting ecological mechanisms including temperature regulation of microbial metabolism and moisture control of substrate diffusion. Accordingly, soil temperature and volumetric water content are modeled as primary mechanistic predictors:

$$ \textbf{x}_{i,t} = (\texttt {SoilTemp}_{i,t}, \texttt {VWC}_{i,t}, \textbf{z}_{i,t}), $$

where \(\textbf{z}_{i,t}\) represents additional ecological influences such as vegetation condition, substrate availability, or microsite heterogeneity that are not directly observed.

Rather than introducing numerous correlated predictors that may induce instability, these secondary influences are incorporated through hierarchical random effects and recovery dynamics, allowing dominant environmental mechanisms to remain interpretable while maintaining model parsimony.

3.4 Hierarchical mixed-effects formulation

The baseline linear predictor for log flux is

$$\begin{aligned} \eta _{i,t} = \beta _0 +\beta _d d_i +\varvec{\beta }_x^\top \textbf{x}_{i,t} +b_{0i}+b_{1i}s_{i,t}, \end{aligned}$$
(2)

where fixed effects quantify average environmental and disturbance responses and random effects capture plot-level ecological heterogeneity.

Random effects follow

$$\begin{aligned} \begin{bmatrix} b_{0i}\\ b_{1i} \end{bmatrix} \sim \mathcal {N} \left( \textbf{0}, \varvec{\Sigma }_b \right) . \end{aligned}$$
(3)

The inclusion of plot-level random effects absorbs persistent spatial variability arising from differences in soil properties, vegetation composition, and local microclimate, thereby reducing bias associated with spatial autocorrelation among nearby plots.

3.5 Disturbance threshold modeling

Forest ecosystems may exhibit nonlinear responses once disturbance exceeds a critical intensity. To detect such structural transitions, we introduce a disturbance threshold \(\delta \in (0,1)\) using a hinge function:

$$ H_\delta (d_i)=(d_i-\delta )_+. $$

The threshold model becomes

$$\begin{aligned} \eta _{i,t} = \beta _0+\beta _d d_i +\gamma H_\delta (d_i) +\varvec{\beta }_x^\top \textbf{x}_{i,t} +b_{0i}+b_{1i}s_{i,t}. \end{aligned}$$
(4)

Below the threshold, disturbance effects follow slope \(\beta _d\), whereas beyond \(\delta \) the response shifts to \(\beta _d+\gamma \), representing nonlinear ecosystem degradation under severe disturbance.

3.6 Post-disturbance recovery dynamics

Temporal ecosystem adjustment following disturbance is modeled using an exponential recovery function:

$$\begin{aligned} R(s_{i,t};\phi ,\rho _i) = \phi \left( 1-e^{-\rho _i s_{i,t}}\right) , \end{aligned}$$
(5)

yielding the full predictor

$$\begin{aligned} \eta _{i,t} = \beta _0+\beta _d d_i+\gamma H_\delta (d_i) +\varvec{\beta }_x^\top \textbf{x}_{i,t} +R(s_{i,t};\phi ,\rho _i) +b_{0i}+b_{1i}s_{i,t}. \end{aligned}$$
(6)

Recovery rates vary with long-term environmental conditions:

$$\begin{aligned} \rho _i= \textrm{softplus} \left( \rho _0+\varvec{\rho }_x^\top \bar{\textbf{x}}_i \right) , \end{aligned}$$
(7)

allowing ecosystem resilience to depend on site characteristics.

3.7 Parameter estimation

The disturbance threshold \(\delta \) is estimated via profile likelihood:

  1. 1.

    evaluate candidate \(\delta \) values over a dense grid,

  2. 2.

    fit mixed-effects models using REML,

  3. 3.

    select \(\widehat{\delta }\) minimizing AIC/BIC,

  4. 4.

    refit the final model at the optimal threshold.

3.8 Model evaluation

Model adequacy is assessed using residual diagnostics, cross-validation performance, calibration analysis, and parametric bootstrap checks. Plot-level validation together with hierarchical random effects ensures that predictive assessment remains robust to potential spatial autocorrelation within the experimental network.

3.9 Ecological interpretation

Model parameters provide mechanistic ecological interpretation. Fixed environmental effects quantify dominant short-term controls on soil respiration, threshold parameters \((\delta ,\gamma )\) characterize ecosystem vulnerability to disturbance intensity, and recovery parameters \((\phi ,\rho _i)\) describe resilience and adjustment rates. Random effects represent aggregated influences of unmeasured ecological processes operating at the plot scale.

4 Results

To improve interpretability, the results are organized around the principal empirical findings rather than sequential descriptions of figures and tables. Each subsection highlights a major research conclusion supported by statistical and graphical evidence.

4.1 Dataset characteristics and experimental structure

We begin by summarizing the FoRTE soil CO\(_2\) flux dataset to establish the empirical context for subsequent modeling. Tables 1–3 describe overall variability, spatial heterogeneity, and treatment balance across observations.

Table 1 summarizes the principal characteristics of the FoRTE soil CO\(_2\) flux observations. Soil respiration exhibits moderate mean flux values accompanied by substantial variability, reflecting heterogeneous environmental conditions and disturbance treatments across plots. Logarithmic transformation reduces dispersion and stabilizes variance, supporting its suitability as the response variable in subsequent hierarchical modeling.

The disturbance severity index spans a wide range of treatment intensities, indicating representation of both low- and high-impact disturbance scenarios. Similarly, time since disturbance varies considerably among observations, enabling analysis of post-disturbance recovery dynamics. Standardized soil temperature and volumetric water content describe relative environmental gradients influencing microbial respiration processes across experimental conditions.

Spatial variability across plots is summarized in Table 2. Mean flux values differ modestly among plots, although variability is notably higher in plot 4. The relatively stable distribution of log-transformed flux across plots suggests comparable underlying ecological processes, whereas differences in disturbance intensity motivate inclusion of plot-level random effects.

Treatment-level summaries presented in Table 3 indicate broadly similar flux magnitudes between baseline (B) and treatment (T) conditions. Comparable disturbance exposure and sampling duration confirm balanced experimental design, implying that observed response differences primarily arise from environmental and disturbance dynamics rather than sampling bias.

Table 1 Summary statistics of the FoRTE soil CO\(_2\) flux dataset
Full size table

Disturbance severity (d) represents the experimentally assigned treatment intensity within the FoRTE design and is expressed as a normalized index ranging from 0 (undisturbed reference condition) to 1 (maximum disturbance severity). The severity metric reflects proportional canopy and structural disturbance imposed during treatment implementation, corresponding to increasing levels of tree removal and mortality simulated to represent natural disturbance processes such as windthrow or harvesting.

Environmental variables (soil temperature and volumetric water content) were standardized to zero mean and unit variance prior to analysis; therefore, reported values are unitless and represent relative deviations from average site conditions.

Table 2 Summary statistics of soil CO\(_2\) flux by plot
Full size table
Table 3 Summary statistics of soil CO\(_2\) flux by treatment
Full size table

4.2 Exploratory patterns of disturbance and environmental controls

Exploratory visualizations provide insight into dominant ecological gradients motivating the threshold modeling framework.

Figure 1 demonstrates a systematic decline in soil CO2 flux with increasing disturbance severity. Undisturbed conditions display both higher median flux and greater variability, suggesting reduced biological activity under intensified disturbance regimes.

Temporal dynamics following disturbance are illustrated in Fig. 2. Soil CO\(_2\) flux exhibits greater variability during early post-disturbance periods, followed by a gradual decline and subsequent stabilization as days since disturbance increase. This pattern suggests progressive adjustment of ecosystem respiration toward more stable post-disturbance conditions rather than an initially elevated flux response.

Environmental controls exhibit contrasting bivariate relationships. Figure 4 shows a strong positive association between soil temperature and flux, indicating enhanced microbial and root respiration under warmer conditions. In contrast, Fig. 3 suggests a negative relationship between soil moisture and flux. However, this apparent effect reflects confounding with soil temperature, as demonstrated by the negative correlation between moisture and temperature in Fig. 7. After adjustment in multivariate models, moisture exhibits a weak positive contribution.

Distributional comparisons across treatments and plots (Figs. 5 and 6) reveal approximately unimodal and symmetric log-flux distributions, supporting modeling assumptions of approximate normality. The correlation structure displayed in Fig. 7 further confirms soil temperature as the dominant environmental driver while highlighting interactions among covariates.

Fig. 1
Fig. 1
Full size image

Soil CO2 flux across disturbance severities (0%, 45%, 65%, and 85%)

Fig. 2
Fig. 2
Full size image

Soil CO2 flux as a function of days since disturbance, colored by disturbance severity

Fig. 3
Fig. 3
Full size image

Soil CO2 flux as a function of scaled soil moisture (VWC), colored by disturbance severity

Fig. 4
Fig. 4
Full size image

Soil CO2 flux as a function of scaled soil temperature, colored by disturbance severity

Fig. 5
Fig. 5
Full size image

Density of log Flux by Treatment, colored by treatment type (B and T)

Fig. 6
Fig. 6
Full size image

Density of log Flux per Plot, colored by factor(plot)

Fig. 7
Fig. 7
Full size image

Pairs plot illustrating correlations and distributions among variables Y, logY, SoilTempScaled, vwcScaled, d, s, p, a, and w

4.3 Threshold model estimation and interpretation

We next quantify disturbance effects using threshold regression models.

Frequentist estimation results reported in Table 4 identify soil temperature as the strongest positive predictor of soil CO\(_2\) flux, followed by smaller positive contributions from soil moisture and time since disturbance. Disturbance intensity exhibits nonlinear behavior, with the threshold component (\(H_{\delta }\)) indicating a structural change in ecosystem response beyond a critical disturbance level.

Although the model explains a moderate proportion of variability (\(R^{2}=0.463\)), such explanatory power is typical for ecological flux observations characterized by substantial spatial heterogeneity and biological variability. Importantly, the statistically significant and directionally consistent effects of disturbance severity, environmental drivers, and recovery time indicate that the model captures the dominant mechanisms governing soil CO\(_2\) dynamics. Thus, the interpretability of estimated parameters supports mechanistic inference regarding disturbance thresholds despite remaining unexplained variability.

The optimal disturbance threshold is determined using the AIC profile shown in Fig. 8, which identifies \(\delta = 0.9\) as providing the best model fit.

Results from the final mixed-effects threshold model are presented in Table 5. Disturbance displays a negative linear effect below the threshold, while the positive threshold adjustment indicates altered system dynamics at higher disturbance severity. Soil temperature remains highly significant, confirming its dominant ecological influence across plots.

Model fit statistics and random effects summarized in Table 6 reveal moderate between-plot variability in baseline flux and recovery trajectories, supporting the hierarchical modeling framework.

Parameter interpretations provided in Table 7 clarify ecological implications. The negative \(\texttt {beta\_d}\) reflects declining flux responses under increasing disturbance below the threshold, whereas the positive \(\gamma \) captures an additional response mechanism activated beyond the critical disturbance level. Environmental covariates contribute positively but remain secondary relative to temperature effects.

Exploratory threshold behavior visualized in Fig. 9 aligns with model-based evidence showing declining flux followed by stabilization at higher disturbance levels. While threshold location is robust, slope magnitude remains sensitive to model specification.

Table 4 Frequentist model with threshold for soil CO\(_2\) flux (dependent variable: logY)
Full size table
Table 5 Summary of Final Mixed-Effects Threshold Model
Full size table
Table 6 Model Fit and Random Effects Summary
Full size table
Table 7 Interpretation of Final Model Parameters
Full size table
Fig. 8
Fig. 8
Full size image

AIC Profile for Threshold (\(\delta \)), with minimum AIC at \(\delta = 0.9\)

Fig. 9
Fig. 9
Full size image

Threshold effect of disturbance on flux, colored by different factors

4.4 Model validation and predictive performance

Model generalization is evaluated using leave-one-plot-out cross-validation (Table 8). Given that the experimental design incorporates exactly four distinct plots, this cross-validation routine involved exactly four independent iterations. Prediction errors remain broadly consistent across all four validation folds, indicating stable predictive performance when individual plots are excluded from model fitting.

For validation, we used the soil respiration dataset from the Forest Resilience Threshold Experiment (FoRTE) data repository, accessed through the fortedata R package (https://fortexperiment.github.io/fortedata/). This dataset contains soil CO\(_2\) efflux observations together with the associated soil temperature and volumetric water content measurements used in the present analysis. After preprocessing, the resulting analytical dataset was used both for model fitting and for validation through held-out plot-based assessment.

For clarity, benchmarking in the present study refers to comparative evaluation of candidate model specifications on the same processed analytical dataset, rather than comparison against a separate external validation dataset. Specifically, candidate threshold formulations were compared using information-criterion-based model comparison during threshold selection, and the predictive performance of the selected model was evaluated using leave-one-plot-out cross-validation. For each of the four folds, one full plot was held out; the model was fitted to the remaining three plots, and predictive error on the held-out plot was summarized using RMSE and MAE.

Because the FoRTE soil respiration dataset provides a unified observational basis for both model development and evaluation, benchmarking is conducted via internal held-out validation rather than a separate external dataset; accordingly, each validation fold assesses performance on plots completely withheld from estimation for that specific iteration.

Validation using an entirely independent, external dataset was not performed because our threshold-based mixed-effects framework is intrinsically coupled with the unique manipulative experimental design of the FoRTE project, which features an engineered, controlled disturbance gradient spanning 0% to 85% severity alongside high-frequency covariate tracking. Finding an external, independent ecosystem study that matches this exact multi-treatment infrastructure, temporal revisit scale, and covariate structure is unfeasible.

Consequently, the absence of an independent external validation dataset represents a structural limitation of this study, restricting our capacity to assert population-level generalizability across different forest systems. Because disturbance severity is assigned at the plot level and only four experimental plots are available, the effective replication for severity and threshold inference is necessarily limited. Model estimates must therefore be interpreted as conditional on the observed FoRTE experimental units.

Under this constraint, the four-fold cross-validation primarily evaluates the internal robustness and stability of inferred relationships rather than the independent replication of disturbance effects. The relatively similar RMSE and MAE values across the held-out plots suggest that estimated environmental responses and recovery dynamics are not dominated by any single experimental unit. Accordingly, the identified disturbance threshold should be interpreted as an exploratory indication of nonlinear ecosystem response within the studied experimental framework rather than definitive evidence of a universal critical threshold.

Residual diagnostics shown in Fig. 10 exhibit an approximately symmetric distribution centered near zero, suggesting reasonable adherence to model assumptions and acceptable calibration performance.

Finally, Fig. 11 illustrates modeled recovery trajectories following disturbance under contrasting moisture regimes. Soil CO\(_2\) flux increases progressively over time since disturbance under both conditions; however, recovery occurs more rapidly under high soil moisture (red line) than under low moisture (blue line), indicating environmental modulation of ecosystem recovery rates.

Table 8 Leave-One-Plot-Out Cross-Validation Metrics
Full size table
Fig. 10
Fig. 10
Full size image

Histogram of in-sample residuals (proxy for PIT), where a well-calibrated model should produce uniform residuals

Fig. 11
Fig. 11
Full size image

Post-disturbance recovery dynamics of soil CO\(_2\) flux under contrasting moisture conditions. The blue line represents low soil moisture conditions, while the red line represents high soil moisture conditions. Flux increases over time since disturbance, with higher moisture associated with faster recovery trajectories

5 Discussion

This study examined how disturbance severity and environmental conditions regulate soil CO\(_2\) flux dynamics using hierarchical threshold modeling applied to FoRTE experimental data. Because disturbance treatments are assigned at the plot level and only four experimental plots are available, all interpretations presented below should be viewed as conditional on these experimental units rather than as universal forest responses. Accordingly, results emphasize mechanistic interpretation and internally consistent patterns rather than confirmatory inference.

Before interpreting ecological implications, it is important to distinguish between exploratory visual patterns and adjusted model estimates. Several exploratory figures (e.g., Figs. 3 and 9) display marginal relationships that differ from regression coefficients obtained under multivariate adjustment. These differences arise from covariate interactions and highlight the importance of conditional inference when multiple environmental drivers operate simultaneously.

5.1 Key findings

  1. 1.

    Soil temperature represents the dominant concurrent control on soil CO \(_2\) flux. Across all model specifications, soil temperature exhibits the largest and most precisely estimated effect (Table 5). A one standard deviation increase in temperature corresponds to an approximate multiplicative increase of \(\exp (0.44)\approx 1.55\) in expected flux. This result is consistent with established ecological understanding that microbial metabolism and root respiration respond strongly to thermal conditions. Although the overall model explains a moderate proportion of variability (\(R^2\approx 0.46\); Table 4), the stability and interpretability of temperature effects support the study’s primary mechanistic conclusions.

  2. 2.

    Soil moisture exerts a secondary but ecologically meaningful influence. Bivariate visualization (Fig. 3) suggests decreasing flux with increasing moisture; however, after accounting for temperature, the adjusted regression coefficient becomes small and positive. This apparent sign reversal reflects negative correlation between temperature and moisture (0.. 7). Once thermal effects are controlled, additional moisture slightly enhances respiration, likely through improved substrate diffusion and microbial activity. Thus, moisture primarily modulates rather than governs short-term flux variability.

  3. 3.

    Disturbance responses indicate a possible high-severity threshold, interpreted cautiously. Model comparison using the AIC profile (Fig. 8) identifies a candidate disturbance threshold near \(\hat{\delta }=0.9\). Estimated slopes suggest declining flux with increasing disturbance up to high severity, followed by a structural change in response. However, differences between frequentist and mixed-effects specifications (Tables 4 and 5) indicate sensitivity of threshold direction to model structure. Given that severity contrasts occur among only four plots, threshold detection should therefore be interpreted as exploratory evidence of nonlinear ecosystem behavior within the FoRTE experiment rather than confirmation of a universal ecological breakpoint.

  4. 4.

    Post-disturbance recovery occurs gradually and is environmentally mediated. Time since disturbance shows a small but consistent positive association with flux, indicating gradual ecosystem adjustment. Recovery trajectories illustrated in Fig. 11 demonstrate faster increases in flux under higher soil moisture conditions, suggesting that favorable environmental conditions accelerate recovery processes even when mean temporal trends remain modest.

  5. 5.

    Predictive performance is moderate but robust across plots. Leave-one-plot-out cross-validation (Table 8) yields similar RMSE and MAE values across held-out plots, indicating that inferred relationships are not driven by any single experimental unit. Because validation excludes entire plots, these results primarily demonstrate robustness within the experimental system rather than independent replication across landscapes. Residual diagnostics (Fig. 10) suggest approximate normality with minor deviations, indicating remaining unexplained ecological variability.

5.2 Ecological and management implications

The results highlight several practical insights. First, strong temperature sensitivity implies that warming periods may substantially increase short-term soil carbon emissions following disturbance. Second, although interpreted cautiously, the emergence of nonlinear responses near high disturbance severity suggests that avoiding extreme canopy removal may reduce ecosystem vulnerability during early recovery stages. Third, maintenance of soil moisture through canopy retention or debris management may promote faster post-disturbance recovery even when instantaneous moisture effects appear modest after statistical adjustment.

5.3 Limitations and scope of inference

Four primary limitations warrant emphasis to properly contextualize the scope and applicability of our findings. First, the study is subject to limited experimental replication, as the disturbance severity varies among only four plots within the FoRTE experimental design. This restriction inherently limits the effective sample size for independent macro-level estimation of severity effects and necessitates a cautious, exploratory interpretation of the threshold parameters. Second, the absence of external validation datasets stands as a clear structural limitation regarding the geographic scalability and generalizability of our findings. Because our framework is tightly coupled with the specialized manipulative treatment configuration and high-frequency tracking unique to the FoRTE infrastructure, testing the model against fully independent external data was unfeasible, leaving the absolute threshold bounds unverified outside this specific regional baseline.

Third, our ecological conclusions exhibit explicit model dependence, meaning that alternative hierarchical random-effects or fixed-effects specifications can influence estimated disturbance responses and parameter boundaries, indicating that findings remain conditional upon specific statistical formulations. Fourth, the framework is constrained by unobserved ecological drivers, where key underlying biological factors–such as fine-root substrate availability, canopy-derived carbon dynamics, and localized microsite heterogeneity–are represented indirectly through spatial-temporal random effects rather than explicit, mechanistic predictors. Consequently, the results presented herein should be interpreted as localized mechanistic insights and hypotheses derived from a highly controlled, specific experimental system rather than as generalized regional population estimates.

5.4 Future methodological directions

Future work could strengthen inference through: (a) smooth nonlinear severity functions replacing a single hinge threshold; (b) state-space or seasonal process models separating observation and ecological dynamics; and (c) joint modeling of carbon pools and fluxes to better connect disturbance impacts with ecosystem recovery mechanisms across spatial scales.

6 Conclusion

Using FoRTE experimental measurements, we developed a threshold mixed-effects framework to examine how disturbance severity and environmental conditions relate to soil CO\(_2\) flux variability and post-disturbance recovery. Because disturbance treatments are assigned at the plot level and inference relies on four experimental plots, conclusions should be interpreted as conditional on this experimental system rather than as broadly generalizable ecosystem responses.

Within this context, three main findings emerge. First, concurrent soil temperature consistently represents the dominant control on short-term soil CO\(_2\) flux variability, exceeding the influence of other measured covariates. Second, model comparison identifies evidence for a potential high-severity threshold (\(\hat{\delta }\approx 0.9\)) at which disturbance–flux relationships change; however, the estimated direction and magnitude of this effect vary across model specifications, and threshold identification should therefore be regarded as exploratory. Third, time since disturbance exhibits gradual recovery trends, with scenario analyses suggesting that higher soil moisture conditions may facilitate faster recovery dynamics when temperature is favorable.

From a management perspective, these results suggest that limiting extreme disturbance intensity and maintaining post-disturbance microclimatic buffering–particularly soil moisture retention–may support recovery processes within systems similar to the FoRTE experiment. However, given the limited spatial replication, these implications should be viewed as hypothesis-generating rather than prescriptive recommendations.

Methodologically, the analysis demonstrates the value of hierarchical threshold models for integrating disturbance gradients with environmental variability while accounting for plot-level dependence. Diagnostic results further indicate opportunities for improvement through smoother nonlinear severity responses, dynamic process-based recovery models, and expanded spatial replication. Future work will focus on incorporating additional FoRTE sites and measurement campaigns, jointly modeling carbon pools and fluxes, and conducting validation using fully independent temporal or spatial holdouts to strengthen robustness and causal interpretation.