the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 12 Jun 2026
Modelling and Interpreting Thermal Stability Indices to Understand Soil Carbon Stabilization Using Soil Properties Data
Abstract. Soil organic carbon (SOC) sequestration and nutrient cycling are related to the susceptibility of soil organic matter to biological decomposition. Several studies have demonstrated associations between biological stability and thermal stability, as assessed using programmed pyrolysis. We sought to develop parsimonious machine learning (ML) models to predict SOC stability indices from measured soil properties. The study analyzed indices such as S1, S2, and S3; the oxygen index (OI); the hydrogen index (HI); and T50, which reflect SOC composition, thermal behaviour, and stability. A total of 203 soil samples collected at 0–15cm depth increments from agricultural, forest, and wetland landscapes in New Brunswick, Canada, were analyzed. Feature selection techniques optimized predictive models, and a random forest (RF) was used to develop one. Correlation results revealed that HI was negatively associated with pH (r = –0.35) and bulk density (r = –0.33), whereas OI showed a positive correlation with pH (r = 0.34). Thermal indices were more strongly related to soil chemistry and texture, with S1 closely correlated with total carbon (r = 0.88) and nitrogen, and S2 negatively associated with sand (r = –0.66) but positively associated with clay (r = 0.24) and POXC (r = 0.84). T50 showed positive correlations with both pH (r = 0.48) and bulk density (r = 0.36), indicating greater thermal stability in higher pH, compacted soils, though these patterns varied by land use. Random Forest (RF) model predicted S1 S3 indices with high accuracy (CCC = 0.83–0.86), while HI and OI were more difficult to model (CCC = 0.44–0.48), suggesting missing biological or environmental predictors; NH₄⁺ and POXC emerged as key predictors. Structural equation modeling (SEM), after addressing multicollinearity, supported a hypothesis driven model that explained ~54% of T50 variation. Clay dissolved organic carbon, pH, and aluminum showed significant direct associations with T50 (β for pH = 0.44), whereas bulk density showed no meaningful relationship. Our study demonstrates that ML and SEM can reveal patterns and associations between soil properties and thermal stability indices, offering insight into understanding the SOC stability under a changing climate as well as presenting a framework for rapid estimation of SOC stability proxies.
- Preprint
(4813 KB) - Metadata XML
-
Supplement
(1372 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-2849', Anonymous Referee #1, 19 Jun 2026
-
AC1: 'Reply on RC1', Kingsley John, 29 Oct 2026
Thank you for the opportunity to revise our manuscript. We are grateful to Reviewer #1 for the thorough and constructive assessment. Below, we provide a point-by-point response to each comment.
General Assessment
"The manuscript represents a scientific work of both good quality and significance. The presentation of this work, however, needs to be revised and, at times, restructured."
Response: We thank the reviewer for this positive assessment. We will undertake a comprehensive revision to improve clarity, reduce redundancy, and strengthen methodological reporting.
Comment 1: RF Modeling Section — Detailing Parameters and Libraries
"In the modeling section, more effort should be put in detailing the exact parameters/libraries used, and explanations of common methods (like PDPs and ICEs) should be minimal and concise. A significant effort has been put into making the paper accessible to a wider audience, however the presentation of many modeling choices was omitted in doing so."
Response: We agree and will substantially revise Sections 2.3.2 and 2.4. We will report: the caret package (Kuhn, 2007) with the underlying randomForest implementation; repeated 10-fold cross-validation with 20 repetitions; grid search for mtry from 1 to the total number of predictors at intervals of two variables; default ntree and tree-growing settings; and the specific R functions used for variable importance (caret::varImp), PDPs/ICEs (iml::FeatureEffects, method = "pdp+ice"), and SEM (lavaan::sem). We will condense the general explanations of PDPs and ICEs and remove the duplicated description of how PDPs are computed.
Comment 2: Code and Data Availability
"Regarding code and data availability, authors have declared that the data will be made available upon reasonable request (due to privacy concerns). Will the code be accessible online? If not, even more precisions need to be provided."
Response: We will revise the Code and Data Availability section to clarify:
- Data: Not publicly available due to privacy concerns of participants (landowner agreements and Canadian Soil Biodiversity Observatory program restrictions). Available upon reasonable request from the corresponding author.
- Code: The R code used for all analyses (VIF analysis, RF modeling, RFE, PDP/ICE generation, SEM, and cross-validation) will be available upon reasonable request from the corresponding author.
Revised text (lines ~835–838):
"Data will be shared upon reasonable request, though no online repository is used to safeguard participant privacy. The R code for all analyses may also be obtained from the corresponding author upon request."
Comment 3: Metrics — Systematic Introduction, References, and Equations
"Metrics have to be systematically introduced and, when appropriate, references and/or equations should be provided. Certain metrics are never introduced, others change their nomenclature throughout the text, and some will depend on the code implementation used (thus, rendering code availability even more important)."
Response: We will systematically revise all metric reporting:
Metric
Action to Be Taken
R²
Will be defined in Section 2.3.2 as the squared Pearson correlation between observed and predicted values, equivalent to R² = 1 − RSS/TSS. Will be capitalized consistently throughout (R²).
Adjusted R²
Equation will be added: R²adj = 1 − [(1 − R²)(n − 1)]/(n − k − 1). A note will explain why values are identical to R² (large n = 203, small k = 4–8; difference < 0.005 when rounded to two decimals).
CCC
Will be defined in Section 2.3.2 as Lin's concordance correlation coefficient.
RMSE
Will be defined in Section 2.3.2 as root mean square error, used as the primary metric for predictive accuracy in cross-validation.
RMSEA
Will be defined in Section 2.4.3 as root mean square error of approximation, used exclusively for SEM model fit. An explicit statement will distinguish RMSE (prediction) from RMSEA (fit).
Bias
Will be defined as mean of (predicted − observed).
MAE
Will be defined as mean absolute error.
AIC/BIC/CFI/SRMR
Will be defined in Section 2.4.3 with threshold criteria and references (Hu & Bentler, 1999; Shi et al., 2019).
Pearson r
Will be explicitly named in Section 2.2.3 (correlation analysis) and Section 3.5 (residual correlations).
Comment 4: Results Section — Repetition and Redundancy
"In general, the results section (3.1, 3.2 …) are very repetitive and contain lots of redundant information. The authors should consider rephrasing and shortening their result explanations so as to not repeat certain conclusions many times."
Response: We will substantially condense Sections 3.1–3.3:
- Section 3.1: Will be reduced from ~3 paragraphs to 2. Repeated descriptions of land-use patterns already covered in the Supplementary Material will be removed.
- Section 3.2 (S₁–T₅₀): Duplicated correlation values and repeated interpretive statements will be removed. Each subsection will report the key correlations once and interpret them concisely.
- Section 3.3 (RF models): The six model subsections (3.3.1–3.3.6) will be shortened by removing repeated descriptions of the same PDP/ICE patterns. Where patterns are consistent across indices (e.g., POXC positive, BD negative), we will refer back to earlier subsections rather than repeating the full interpretation.
Example of redundancy to be removed:
In Section 3.3.5 (HI Model), the sentence "In Fig. 7 for the PDP plots HI model, C:N exhibited a relatively stable trend…" will be replaced with the concise: "The PDP for the HI model (Fig. 7) shows that C:N exhibited a relatively stable trend with minor fluctuations, indicating a minimal effect on HI." Duplicated sentences in Section 3.3.2 (S₂ Model) describing variable importance will be consolidated into a single statement.
Comment 5: Distinguishing Own Results from Consensus
"In the results section, and particularly in the summary results of the thermal indicators of stability, authors need to make clearer distinctions between the conclusions from their own results and the overall consensus regarding the relationships between thermal indicators and SOM stability."
Response: We will revise the Results and Discussion to explicitly separate our findings from established literature. Specifically:
- In Section 3.1, statements about land-use patterns will begin with "In our dataset…" or "Our results show…", followed by "This is consistent with…" when linking to literature.
- In Section 3.2, each correlation will be presented as "We observed…" before any literature-based interpretation.
- In Section 3.5 (SEM), we will strengthen the existing distinction by explicitly labelling hypothesis tests as "our a priori predictions" and comparing them to "consistent with established theory (e.g., Six et al., 2004)."
Comment 6: Interpretation of PDP and ICE Results
"Authors need to interpret PDP and ICE results with a grain of salt, as these can't always be taken at face value due to the marginal effects they represent (as the authors have well explained in the modeling section)."
Response: We agree and will add a caveat in Section 3.3 and in the Conclusions:
"We caution that PDP and ICE curves represent marginal effects — that is, the average effect of a predictor while holding all other predictors constant. These plots therefore illustrate model behavior rather than direct causal relationships, and interactions among predictors may not be fully captured. We interpret these patterns as indicative of associations consistent with known biogeochemical mechanisms, not as definitive evidence of causal effects."
This caveat will be added at the start of Section 3.3 and echoed in the Conclusions (lines ~815–820).
Comment 7: Consistency in Introduction — Repetition, Capitalization, Acronyms
"I urge the authors to pay attention to having consistency in the introduction, repetition, and capitalization of acronyms and initialisms. Same goes for the introduction of new terms, chemical solutions, and references."
Response: We will conduct a full consistency check:
- Acronyms: All acronyms will be defined at first use and used consistently thereafter (e.g., SOC, SOM, OM, RF, SEM, PDP, ICE, VIF, RFE, CEC, POXC, DOC, BD, HI, OI, T₅₀).
- Capitalization: "Random Forest" → "random forest" (except when referring to the proper noun algorithm name), "Structural Equation Modeling" → "structural equation modeling" (except in headings).
- Chemical notation: Standardized subscripts and superscripts (e.g., S₁, S₂, S₃, T₅₀, NH₄⁺, NO₃⁻, CO₂, P₂O₅, K₂O).
- Repetition: Duplicated sentences in the Introduction will be removed (e.g., the repeated description of Fang et al., 2021 and An et al., 2023 in lines 78–84 of the tracked manuscript).
- References: Formatting inconsistencies will be corrected (e.g., An et al., 2023 missing period; Kögel-Knabner spelling; Kasraei vs. Kasrae).
Comment 8: Specific Technical Corrections from the Annotated PDF
"Detailed corrections are provided in a commented pdf."
Response: We will address all technical corrections marked in the annotated PDF. Below is a summary of the key corrections to be made (with line references to the tracked-changes manuscript):
PDF Comment
Location
Correction to Be Made
"Reads like an oral presentation" — PDP plots
Section 3.3.5
Will be rephrased to "The PDP for the HI model (Fig. 7) shows…" (removing "PDP plots" redundancy)
"How does it indicate this?" — S₂ variable importance
Section 3.3.2
Mechanistic explanation will be added linking POXC/DOC to labile C pools and BD/CEC to physical/mineral properties
"Is this a LOWESS fit?" — Figure 9
Figure 9 caption
Will be clarified that the left panel shows OLS linear regression (not LOWESS); the right panel uses LOWESS for residuals
"If you use both RMSE and RMSEA it needs to be justified"
Sections 2.3.2 & 2.4.3
Explicit distinction will be added: RMSE = predictive accuracy; RMSEA = SEM model fit
"The hypotheses were not introduced as H1–H5"
Section 2.4.3
H1–H5 will be confirmed as introduced in Methods (lines 326–341); a clarifying list format will be added
"Which correlations? Spearman, Pearson…?"
Section 3.5
"Pearson residual correlations" will be specified
"Saturated nature of our model (df = 0)?"
Section 3.5
The "saturated" reference will be removed; replaced with an explanation that excellent fit reflects theoretical grounding, not causation
"Consider replacing with relative RMSE"
Table 2
RRMSE column and note will be added
"R² and Adjusted R² values are the same"
Table 2
A note will be added explaining rounding due to large n and small k
"R² often capitalized"
Throughout
Will be standardized to R²
"Provide equation for R²"
Section 2.3.2
R² equation and explanation of caret implementation will be added
"Further emphasized that SOC stability…"
Section 3.3.6
Will be rephrased as a complete sentence with correct citation format
Closing Statement
We thank Reviewer #1 for their thorough and constructive feedback, which will substantially improve the clarity, rigor, and readability of our manuscript. All comments will be addressed, and we believe the revised manuscript will be suitable for publication in SOIL.
-
AC1: 'Reply on RC1', Kingsley John, 29 Oct 2026
-
CC1: 'Comment on egusphere-2026-2849', Pierre Barré, 28 Jul 2026
This study presents the results of thermal analyses of more than 200 topsoil samples collected in the New Brunswick region of Canada. The authors attempted to predict some selected thermal variables from soil characteristics using ML approaches and to identify the most influential predictors of each thermal variable. Finally, they used SEM to explore causal pathways linking soil properties to thermal stability. This study is original, particularly with regard to the causal models. I recommend its publication following a few minor comments that can nonetheless take some time to address.
Previous studies have linked biogeochemical stability observed in the field to thermal indicators (Barré et al., 2016; Cécillon et al., 2021). These studies could have served as a basis for selecting the indicators to be considered in the present draft. Furthermore, the T50 indicator is not described in sufficient detail. Is this the T50 calculated from S2 or from the entire pyrolysis signal? If it is the T50 measured from the S2 signal, that is unfortunate because it is really not the best indicator of biogeochemical stability (Cécillon et al., 2021). If it is the T50 calculated from the entire pyrolysis signal, it is necessary to better explain how this signal is constructed so that it incorporates both S2 and S3 (in particular, how were the CO2 and HC signals converted into carbon in order to determine the temperature at which half of the pyrolyzable C was emitted?)
The authors considered numerous soil characteristics to predict thermal indicators. Why didn’t they consider total organic carbon? Total organic carbon is a good indicator of biogeochemical stability (Zhao et al., 2026). I believe that adding TOC should improve the models’ metrics.
Delahaie et al. (2023) presented a study reporting the results of more than 1,700 thermal analyses and linking them to pedoclimatic variables. Delahaie et al. (2023) is cited in the introduction but is not referenced again in the discussion of the results. This is regrettable given the close relationship between the two studies.
Given the model metrics, at this stage, I would continue to recommend actually conducting the thermal measurement rather than predicting the result based on soil characteristics. It would be interesting for the authors to discuss this point.
Cited references
Barré et al. (2016) The energetic and chemical signatures of persistent soil organic matter. Biogeochemistry, 130 :1-12.
Cécillon et al. (2021) Partitioning soil organic carbon into its centennially stable and active fractions with machine-learning models based on Rock-Eval® thermal analysis. Geoscientific Model Development, 14 :3879-3898.
Zhao et al. (2026) Declines in organic carbon persistence with increased soil carbon. Nature communications. https://doi.org/10.1038/s41467-026-75185-4
Citation: https://doi.org/10.5194/egusphere-2026-2849-CC1 -
AC3: 'Reply on CC1', Kingsley John, 29 Oct 2026
Response to Reviewer CC1 (Pierre Barré)
We thank Dr. Barré for the positive assessment and for recognizing the originality of our causal modeling approach. Below we address each comment.
Comment 1: Indicator Selection — Link to Biogeochemical Stability Literature
"Previous studies have linked biogeochemical stability to thermal indicators (Barré et al., 2016; Cécillon et al., 2021). These studies could have served as a basis for selecting indicators."
Response: We agree. We will add citations to Barré et al. (2016) and Cécillon et al. (2021) in the Introduction and will explicitly justify our indicator selection (S₁, S₂, S₃, HI, OI, T₅₀) based on their established relationships with SOM persistence.
Revised text (Introduction; lines 68–72): We will add:
"Thermal analysis has emerged as a robust proxy for biogeochemical stability, with Rock-Eval indices shown to correlate with long-term carbon persistence (Barré et al., 2016; Cécillon et al., 2021; Gillespie et al., 2014; Mesgar et al., 2024b). Building on this foundation, we selected indicators capturing complementary dimensions of stability: S₁, S₂, S₃ (thermal pools of varying lability), HI and OI (OM quality and diagenetic state), and T₅₀ (continuous thermal stability)."
Comment 2: T₅₀ — Insufficient Detail
"Is this the T₅₀ calculated from S₂ or the entire pyrolysis signal? If S₂-derived, it is not the best indicator (Cécillon et al., 2021)."
Response: We confirm that our T₅₀ is calculated from the HC part of the pyrolysis signal (S1 and S2), not from S₂ alone. We will add explicit methodology.
Revised text (Section 2.2.3):
"Using the methodology of Gregorich et al. (2015), T₅₀ was calculated by integrating the area under the S1 and S2 curves (i.e., total HC). This approach follows the methodology of Gregorich et al. (2015). The combined cumulative carbon release was then normalized to total pyrolyzable carbon, and T₅₀ was interpreted as the temperature at which 50% of this cumulative carbon signal had been emitted (Barros et al., 2011; Gregorich et al., 2015). This approach follows the methodology of Gregorich et al. (2015) and Cécillon et al. (2018) and captures the thermal stability of the full pyrolyzable organic matter pool, not just the HC fraction, thereby providing a more robust proxy for biogeochemical stability (Cécillon et al., 2021)."
Comment 3: Why Not Total Organic Carbon (TOC)?
"Why didn't they consider TOC? TOC is a good indicator of biogeochemical stability (Zhao et al., 2026). Adding TOC should improve model metrics."
Response: TOC (measured as TC) was removed due to collinearity with POXC, DOC, TN, and C:N. We deliberately exclude it for two reasons:
- Data leakage: HI = S₂/TOC and OI = S₃/TOC are mathematically derived from TOC. Including TOC would inflate accuracy metrics unrealistically.
- Skewed interpretation: TOC would dominate variable importance, obscuring mechanistic predictors (POXC, DOC, CEC, Al, BD) central to our SEM.
- The selection process was computational, using variance inflation factor (VIF) analysis and recursive feature elimination (RFE).
Revised text (Section 2.3.3):
"TOC is a well-established stability indicator (Barré et al., 2016; Zhao et al., 2026) but was excluded to avoid data leakage and preserve interpretability."
Comment 4: Delahaie et al. (2023) Not Discussed
"Delahaie et al. (2023) is cited in the introduction but not referenced again in the discussion. This is regrettable given the close relationship between the two studies."
Response: We agree and will add a dedicated passage in Section 3.4 comparing our findings with Delahaie et al. (2023).
Revised text (Section 3.4):
"In addition, Delahaie et al. (2023) observed a positive relationship between T₅₀ and pH under acidic conditions, attributed to reduced microbial activity, as it slows organic matter mineralization and favors the accumulation of labile, thermally unstable SOC — a mechanism that aligns with our interpretation of reduced decomposition rates in lower-pH soils. However, we caution that our single time-point sampling cannot establish temporal precedence, a key criterion for causal inference (Hill, 1965; Wade et al., 2026). Thus, while our results are consistent with the hypothesized directional relationships, they represent interaction evidence that contributes to, but does not definitively establish, causal understanding. For example, Delahaie et al. (2023) observed weak but significant positive correlations between thermal stability and mean annual temperature (MAT), whereas our study did not include climatic covariates, which we identify as a priority for future work."
Comment 5: Recommend Direct Measurement Over Prediction
"At this stage, I would continue to recommend actually conducting the thermal measurement rather than predicting from soil characteristics. It would be interesting to discuss this point."
Response: We agree that direct thermal measurement remains the gold standard. We will add a paragraph clarifying that our models are a screening tool, not a replacement.
Revised text (Recommendation section):
"We note that our predictive models are not intended to replace direct thermal measurement. With the current accuracy levels, for example, with HI and OI (CCC = 0.44–0.48), modelled values cannot substitute for Rock-Eval or HAWK comprehensive analysis when precise stability estimates are required. Rather, our framework is presented to show the relationship between thermal indicators and soil characteristics and serve as a screening and standard tool: it can identify soils or regions where thermal analysis is most informative, and it could provide approximate stability proxies in large-scale monitoring where direct measurement is infeasible. We recommend a hybrid approach: direct thermal measurement on a representative subset of samples, followed by predictive extrapolation to larger sample sizes using the validated models developed here, consistent with established soil transfer models (Bouma, 1989; Van Looy et al., 2017)."
We thank Dr. Barré for these constructive comments, which will improve the manuscript's clarity and rigor.
Regards,
Kingsley John (on behalf of all authors)Citation: https://doi.org/10.5194/egusphere-2026-2849-AC3
-
AC3: 'Reply on CC1', Kingsley John, 29 Oct 2026
-
RC2: 'Comment on egusphere-2026-2849', Teneille Nel, 08 Oct 2026
General Comments
Thank you for the opportunity to review this manuscript. This is a comprehensive and well-designed study investigating the potential drivers of soil organic matter stability using thermal techniques across soils and land-use types in Eastern Canada. The manuscript addresses an important topic, and the combination of thermal indicators, soil physicochemical properties, correlation analyses, and machine-learning approaches provides a potentially valuable framework for identifying the factors associated with SOM persistence.
I particularly appreciated the thorough statistical analyses and the attention given to model diagnostics. The manuscript also makes a useful contribution by considering multiple indicators of SOM stability rather than relying on a single metric. However, I think the manuscript would benefit from a stronger mechanistic interpretation of the observed relationships. Soil stabilization mechanisms are frequently mentioned throughout the manuscript, but they are not always explained in sufficient detail or explicitly connected to the soil properties and soil types represented in this dataset.
In particular, I encourage the authors to consider whether the relative importance of the identified drivers differs among soil types and/or land-use types. Relationships identified across the entire dataset may reflect differences among soil groups rather than consistent relationships operating within each soil type. Exploring this possibility, where statistically feasible, would substantially strengthen the interpretation and help distinguish general controls on SOM stability from relationships that are specific to particular soil environments.
Specific Comments
Introduction
L22: Please define the S1–S3 parameters when they are first introduced.
L28: Please clarify: more strongly than what?
L73: Regarding the statement about “integrating their outputs,” I suggest mentioning Bayesian approaches as a potential avenue for future research, which could may provide a useful framework for simultaneously evaluating several complementary indicators of SOM stability.
Materials and Methods
Figure 1: Please ensure that there is at least one city or other recognizable location that is named legibly on the map.
L161–175: Is a summary of the physical and chemical properties of the soils provided elsewhere in the manuscript or supplementary material? If not, I strongly recommend including one. A table reporting the mean, range, and standard deviation of the major soil physical and chemical properties for the overall dataset, and preferably by soil type and/or land use, would substantially improve the manuscript.
Information about how the absolute range over which those properties vary, is important for interpreting the subsequent statistical relationships between soil properties and SOM stability indices. These summary statistics would also help establish whether the proposed physical and chemical processes are realistically plausible within the soils investigated.
L191–193: Please clarify whether the use of HI and OI as predictors in the machine-learning models introduces circularity in the machine-learning models, since HI and OI are derived from S1–S3 measurements and TOC content?
Results and discussion
Figure 3: Please indicate in the figure caption that these plots represent the full dataset, including all land-use types. This is important for interpreting the relationships shown. Please also add the appropriate units to the x-axis labels of the CDF plots.
Section 3.2: This section is currently quite lengthy and contains many individual correlation results. I suggest focusing the Results section on the clearest and most biologically or mechanistically meaningful relationships rather than reporting all correlations in detail.
An alternative approach would be to move some of the interpretation of individual correlations into the Discussion and integrate the correlation results with the machine-learning results. This would reduce repetition and allow the manuscript to develop a clearer mechanistic narrative. For example, correlations could be used primarily to establish broad relationships, while the discussion could then consider whether these relationships are supported by the predictive models and what mechanisms might explain them.
L382: The phrase “protecting the release of CO₂” is unclear in terms of functional interpretation. I suggest rephrasing this in terms of the established concept of clay-associated protection of organic matter from mineralization or decomposition.
L391: The statements regarding pH appear potentially contradictory. The manuscript states that “alkaline conditions favor oxidative transformations,” whereas later (L430) it states that higher pH promotes SOC preservation. Please reconcile these statements and clarify whether the effects of pH differ depending on the SOM fraction, the specific thermal stability indicator, or the mechanism being considered. It may be that pH has opposing effects on different processes. If so, this distinction should be made explicit rather than presenting the two statements as generally applicable.
L404–405: Please clarify the mechanism by which organo-mineral interactions are proposed to enhance oxidative transformation of SOM. If the intended mechanism involves preferential stabilization or transformation of particular organic compounds, altered microbial accessibility, catalytic effects of minerals, or selective preservation of oxidation-resistant compounds, please explain this more explicitly.
L413: Please specify how bulk density (BD) is expected to alter SOM protection mechanisms. For example, is the proposed relationship related to changes in pore structure, aggregate formation, oxygen diffusion, physical accessibility of SOM, or simply the association of BD with soil texture/mineralogy? It would be helpful to distinguish direct effects of BD from its potential role as a proxy for other soil properties.
L427: The relationship between CEC and SOM stability could potentially have an alternative explanation. Higher CEC may be associated with greater soil fertility and therefore greater plant productivity or fresh organic matter inputs. Thus, the observed relationship may not necessarily indicate a direct protective effect of CEC on SOM. If possible, discuss whether the available data allow these mechanisms to be distinguished.
L494: Please provide more detail regarding the statement that “interactions from N content” influence SOM stability. What specific interactions are being referred to? A more explicit mechanistic explanation would strengthen this section.
L510: Please clarify the proposed mechanisms through which trace metals influence SOM stability. Are the authors referring to direct complexation with organic compounds, association with mineral surfaces, catalytic effects on oxidative reactions, effects on microbial activity, or another mechanism? Please provide a more specific mechanistic explanation and appropriate references.
Organization of the Discussion
Sections 3.3.1, 3.3.2, and 3.3.3: These sections contain considerable repetition and discuss several similar mechanisms across the different models and thermal stability indices. I suggest restructuring these sections to provide a more integrated discussion. Rather than discussing each model or SOM fraction independently, the authors could first identify the major mechanisms consistently supported across multiple models/fractions. The authors could then discuss the exceptions and differences among individual indices. This type of organization would make the discussion more concise while also making the overarching conclusions easier to identify. It would also help avoid repeatedly explaining the same mechanism when it appears as an important predictor of several stability indicators.
Section 3.3.6: Please reduce repetition in this section. Consider focusing this section on the unique findings or implications that have not already been addressed.
Minor Editorial Comment
L593–594: This is an incomplete sentence. Please revise for grammatical completeness.
Overall Recommendation
Overall, I found this to be a strong and potentially valuable study. The experimental/data framework and statistical analyses are substantial, and the manuscript has the potential to make a useful contribution to our understanding of the controls on SOM stability as assessed through thermal techniques.
The main opportunity for improvement is to move beyond identifying statistical associations and develop a more integrated mechanistic interpretation of the results. Greater attention to the ranges of soil properties, differences among soil types and land uses, and the specific mechanisms underlying the observed relationships would considerably strengthen the manuscript.
I also recommend reducing repetition in the Results and Discussion and synthesizing the common mechanisms across SOM fractions and predictive models. Doing so would make the manuscript more concise while highlighting the most important findings and improving the overall narrative.
Citation: https://doi.org/10.5194/egusphere-2026-2849-RC2 -
AC2: 'Reply on RC2', Kingsley John, 29 Oct 2026
Response to RC2
We thank the reviewer for their thorough and constructive assessment. We will carefully address each comment, revise the manuscript accordingly, and add supporting material where needed. Below, we provide point-by-point responses. Line numbers refer to the revised manuscript unless otherwise stated.
General Comments
"The manuscript would benefit from a stronger mechanistic interpretation, and the authors should consider whether the relative importance of drivers differs among soil types and/or land-use types. Relationships identified across the entire dataset may reflect differences among soil groups rather than consistent relationships operating within each soil type."
Response: We agree that distinguishing general controls from land-use-specific relationships is essential. We will take the following steps to address this, including new statistical tests to formalize the comparison:
- We will expand the mechanistic explanation of mineral protection in Section 3.2: "The mineral protection works through adsorption, compartmentalization, aggregate-scale isolation, and reduced accessibility to decomposers rather than by literally protecting C itself (Kleber et al., 2021; Islam et al., 2022; Xu et al., 2022)." We will also add: "Furthermore, we must note that clay content is a bulk metric. SOC stabilization depends on the specific mineralogy, surface reactivity, and degree of microbial processing (Islam et al., 2022; Xu et al., 2022), highlighting that mineral protection is highly context-dependent."
- We agree that pooled correlations could, in theory, mask within-group differences. While our study is defined within the New Brunswick landscape, where land-use classes are point attributes rather than independent strata, and the pooled matrix remains our primary landscape-level result, we will carefully inspect the within-land-use matrices (see Figures S10–S13 in the Supplementary Materials) to address this concern. The indicators' relationships discussed in the manuscript (e.g., S₁/S₂ with TC, DOC, and POXC; S₃ with TN and TC; and textural relationships with clay and sand) are strongly consistent across agriculture, forest, pasture, and wetland. However, we will identify relationships that are visibly different, primarily involving dissolved nitrogen species and metals; for example, the negative correlations of T₅₀ with NH₄⁺ and TIN, the internal nitrogen pool correlations (e.g., TDN–NH₄⁺), and the Zn relationships (e.g., S₁–Zn, DOC–Zn) appear noticeably weaker or inconsistent in the pasture and wetland subsets compared to agriculture and forest. We will interpret these divergences cautiously, as the dataset is imbalanced across land-use types (agriculture and forest dominate; wetland and pasture are much smaller), meaning apparent differences in correlation strength are highly likely to reflect unequal sample sizes and within-group variance rather than consistent biological mechanisms operating differently across land covers.
Specific Comments — Introduction
L22: "Please define the S₁–S₃ parameters when they are first introduced."
Response: We will add the following definitions to the Abstract: "The study analyzed indices such as S₁ (volatile hydrocarbons), S₂ [hydrocarbons from thermal cracking of non-volatile organic matter (OM)], S₃ (CO₂ from oxygen-containing functional groups); the oxygen index (OI); the hydrogen index (HI); and T₅₀ (temperature at which 50% of the organic matter has been emitted)."
L28: "Please clarify: more strongly than what?"
Response: We will revise the sentence to read: "Thermal indices were more strongly related to soil chemistry and texture than to other measured soil properties, with S₁ closely correlated with total carbon (r = 0.88) and nitrogen, and S₂ negatively associated with sand (r = −0.66) but positively associated with clay (r = 0.24) and POXC (r = 0.84)."
L73: "Regarding 'integrating their outputs,' I suggest mentioning Bayesian approaches as a potential avenue for future research."
Response: We will add: "Bayesian hierarchical models may offer a promising approach for incorporating several complementary SOM stability indices while explicitly propagating measurement uncertainty and accommodating each technique's different calibrations and error structures (Jia et al., 2025)."
Materials and Methods
Figure 1: "Please ensure that there is at least one city or other recognizable location named legibly on the map."
Response: The map will be updated to include legible location labels.
L161–175: "Is a summary of the physical and chemical properties of the soils provided? A table reporting mean, range, and standard deviation, preferably by soil type and/or land use, would substantially improve the manuscript."
Response: We will provide a comprehensive set of boxplots in the Supplementary Material (Figures S6–S9) showing the distributions and ranges of the predictors (soil physical and chemical properties). We will add an explicit reference to these figures in Section 3.1 to make this information easier to locate.
In addition, we will provide full within-land-use correlation matrices; see Supplementary Materials (Figures S10–S13).
L191–193: "Please clarify whether the use of HI and OI as predictors in the machine-learning models introduces circularity, since HI and OI are derived from S₁–S₃ measurements and TOC content."
Response: We will clarify that HI and OI were used only as target variables in their respective models and never as predictors for S₁, S₂, S₃, or T₅₀. All predictors were independent soil physicochemical properties measured separately from the thermal analysis, precluding circularity. TOC will be excluded from the predictor set (Section 2.3.1) via multicollinearity analysis to avoid the mathematical dependence of HI = S₂/TOC and OI = S₃/TOC on TOC and data-leakage issues in modelling.
Results and Discussion
Figure 3: "Please indicate in the figure caption that these plots represent the full dataset, including all land-use types. Please also add appropriate units to the x-axis labels of the CDF plots."
Response: The figure and caption will be updated accordingly.
Section 3.2: "This section is currently quite lengthy and contains many individual correlation results. I suggest focusing on the clearest and most biologically or mechanistically meaningful relationships, or moving some interpretation into the Discussion and integrating correlation results with the machine-learning results."
Response: We agree and will substantially revise Section 3.2. We will condense the reporting of individual correlations, retain the most mechanistically meaningful relationships, and move broader interpretation into the Discussion, where correlation results will be integrated with the machine-learning findings. This will reduce repetition and strengthen the mechanistic narrative. All cited works will be retained.
L382: "The phrase 'protecting the release of CO₂' is unclear. I suggest rephrasing in terms of clay-associated protection of organic matter from mineralization or decomposition."
Response: We will replace "protecting the release of CO₂" with a description of clay-associated protection of organic matter from microbial decomposition and mineralization. The revised text will distinguish the retention and reduced accessibility of organic matter from the gases released during programmed pyrolysis. We will also clarify that the positive correlation between clay and S₃′ is consistent with differences in organic matter retention but does not independently demonstrate a protective mechanism. The supporting references will be retained.
L391 and L430: "The statements regarding pH appear potentially contradictory ('alkaline conditions favor oxidative transformations' vs. 'higher pH promotes SOC preservation'). Please reconcile."
Response: We agree that the previous wording was overly general. We will remove the statements that alkaline conditions necessarily favor oxidative transformation and that higher pH promotes SOC preservation by reducing microbial decomposition. The revised text will distinguish the thermal indicators: OI describes an oxygen-related signal normalized to TOC, HI describes hydrogen-related composition, and T₅₀ describes thermal resistance under the analytical conditions. In our results, pH was positively correlated with OI and T₅₀ but negatively correlated with HI. These relationships can coexist because the indicators describe different characteristics of SOM; they do not establish that pH has a uniform effect on biological decomposition or preservation. We will now restrict interpretation to the observed pH gradient and acknowledge possible covariation with organic matter composition, mineral properties, and land use. The available correlations cannot determine whether pH has opposing effects on specific processes.
L404–405: "Please clarify the mechanism by which organo-mineral interactions are proposed to enhance oxidative transformation of SOM."
Response: We will remove the implication that the positive OI–Al correlation demonstrates enhanced oxidative transformation through organo-mineral interactions. The revised interpretation will identify preferential association or retention of oxygen-bearing organic compounds as one possible explanation for the observed relationship; such compositional selection would not necessarily require increased oxidation. We will also distinguish mineral-associated retention from catalytic transformation. The correlation data do not identify metal speciation, mineral-surface reactions, or changes in individual organic compounds and therefore cannot discriminate among selective retention, transformation, and indirect environmental associations. The relevant cited works will be retained, while the mechanistic wording will be qualified accordingly.
L413: "Please specify how bulk density (BD) is expected to alter SOM protection mechanisms, and distinguish direct effects from its role as a proxy."
Response: We agree that BD should not be interpreted as direct evidence of a specific protection mechanism. We will revise the text to distinguish its observed correlations from possible explanations involving pore structure, oxygen diffusion, aggregation, and microbial accessibility. BD may also act as a proxy for differences in organic matter concentration, texture, or land use. Because the correlation analysis does not resolve these alternatives, we will no longer describe the BD relationships as demonstrating that compaction enhances SOM protection or promotes mineralization. Detailed consideration of these alternatives will be integrated with the model results in Section 3.3.
L427: "The relationship between CEC and SOM stability could have an alternative explanation (fertility/productivity rather than direct protection)."
Response: We will revise the interpretation to acknowledge that CEC may reflect nutrient retention and conditions associated with plant productivity and fresh organic matter inputs, rather than a direct protective effect on SOM. Organic matter itself can also contribute to CEC, making the direction of the relationship uncertain. The contrasting associations of CEC with pyrolysis yields and T₅₀ further caution against treating CEC as a uniformly protective property. The reported correlations cannot distinguish mineral-associated protection from productivity-related inputs or covariation with organic matter abundance. We will therefore retain the observed relationships while removing the inference that they establish a direct stabilization mechanism.
L494: "Please provide more detail regarding 'interactions from N content' and the specific interactions referred to."
Response: We will expand this passage to specify the mechanisms considered — namely, the role of nitrogen availability in constraining microbial decomposition, the formation of N-bearing organic compounds with differing thermal lability, and the covariation of N with labile carbon pools (POXC, DOC) and C:N stoichiometry. We will state explicitly which interactions are inferred from the correlation structure and which remain speculative.
L510: "Please clarify the proposed mechanisms through which trace metals influence SOM stability."
Response: We will clarify that we refer to (i) direct complexation of metals with organic functional groups, (ii) association with mineral surfaces and co-precipitation, (iii) catalytic or redox effects on oxidative reactions, and (iv) indirect effects on microbial activity. We will specify that the correlation data cannot distinguish among these mechanisms and will qualify the interpretation accordingly, with appropriate references retained.
Organization of the Discussion
Sections 3.3.1–3.3.3: "These sections contain considerable repetition. I suggest restructuring to provide a more integrated discussion — first identifying mechanisms consistently supported across models/fractions, then discussing exceptions."
Response: We will restructure these sections as recommended. The revised Discussion will first synthesize the major mechanisms supported consistently across multiple models and thermal indices, then discuss the exceptions and index-specific differences. This will reduce repetition, make the overarching conclusions easier to identify, and avoid re-explaining the same mechanism for each indicator.
Section 3.3.6: "Please reduce repetition and focus on unique findings or implications."
Response: This section will be condensed to focus on findings and implications not already addressed elsewhere, with repeated material removed.
Minor Editorial Comment
L593–594: "This is an incomplete sentence. Please revise for grammatical completeness."
Response: Thank you. This will be corrected.
Overall Recommendation
"Overall, I found this to be a strong and potentially valuable study… The main opportunity for improvement is to move beyond identifying statistical associations and develop a more integrated mechanistic interpretation of the results."
Response: We thank the reviewer for recognizing the strength of the study and for the constructive suggestions. The revisions will strengthen the approach, clarity, and interpretation of the manuscript.
Please see the attached for the proposed figure updates.
Regards,
Kingsley John (on behalf of all authors)
-
AC2: 'Reply on RC2', Kingsley John, 29 Oct 2026
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 280 | 122 | 43 | 445 | 66 | 32 | 30 |
- HTML: 280
- PDF: 122
- XML: 43
- Total: 445
- Supplement: 66
- BibTeX: 32
- EndNote: 30
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Using a dataset of n=203 samples from the New Brunswick province in Eastern Canada, the authors proposed a detailed modeling and analysis framework linking thermal stability indicators with soil organic matter (SOM) stability. The framework includes both a predictive model (using a Random Forest (RF) implementation in R), and a more analytical model that aims at removing the black-box quality so often associated with machine learning predictions (using Structural Equation Modeling (SEM)). The manuscript represents a scientific work of both good quality and significance. The presentation of this work, however, needs to be revised and, at times, restructured.
The RF modeling section presents a serious effort of finding the most parsimonious set of predictors, by both eliminating (multi)collinearity in the feature set before model training, and by applying recursive feature elimination to further eliminate less informative features. During the RF model training and selection, a grid search was used to find the best-fitting hyperparameters. The authors also included a posteriori analysis using partial dependence plots and individual conditional expectation plots. A detailed analysis of the relationships between the main thermal indicators (S1, S2, S3, OI, HI, T50) is introduced before presenting the modeling results, confirming pre-existing hypotheses about the links between SOM stability and thermal stability indicators.
In the modeling section, more effort should be put in detailing the exact parameters/libraries used, and explanations of common methods (like PDPs and ICEs) should be minimal and concise. A significant effort has been put into making the paper accessible to a wider audience, however the presentation of many modeling choices was omitted in doing so.
Regarding code and data availability, authors have declared that the data will be made available upon reasonable request (due to privacy concerns). Will the code be accessible online? If not, even more precisions need to be provided.
Metrics have to be systematically introduced and, when appropriate, references and/or equations should be provided. Certain metrics are never introduced, others change their nomenclature throughout the text, and some will depend on the code implementation used (thus, rendering code availability even more important).
In the results section, and particularly in the summary results of the thermal indicators of stability, authors need to make clearer distinctions between the conclusions from their own results and the overall consensus regarding the relationships between thermal indicators and SOM stability. In general, the results section (3.1, 3.2…) are very repetitive and contain lots of redundant information. The authors should consider rephrasing and shortening their result explanations so as to not repeat certain conclusions many times. Similar to the modeling section, the results are at times very elaborate while significant statistical/modeling choices are omitted.
Authors need to interpret PDP and ICE results with a grain of salt, as these can’t always be taken at face value due to the marginal effects they represent (as the authors have well explained in the modeling section).
Detailed corrections are provided in a commented pdf. I urge the authors to pay attention to having consistency in the introduction, repetition, and capitalization of acronyms and initialisms. Same goes for the introduction of new terms, chemical solutions, and references.