Technical note: A validated correction method to quantify organic and inorganic carbon in soils using Rock-Eval® thermal analysis

Soils contain large amounts of carbon stored as organic carbon and carbonates. These carbon pools can contribute to climate regulation and are of primary importance in ensuring proper soil functioning. However, their accurate quantification remains a complex task. Rock-Eval^® thermal analysis has emerged as an alternative to classic dry combustion and wet methods due to its ability to simultaneously provide organic and inorganic carbon measurements on the same subsample. However, it has been observed that Rock-Eval^® systematically underestimates the soil organic carbon (SOC) while overestimating the soil inorganic carbon (SIC). In this technical note, we propose a validated correction of both SOC and SIC based on a machine-learning model and using a diverse dataset of 240 soil samples. We show that the proposed correction significantly increases the accuracy of the Rock-Eval^® method when compared to reference SOC and SIC values and applied to the dataset used for training and testing and that it can be successfully applied to data originating from different Rock-Eval^® machines without changing the routine analytical protocol. The transferability of the model allows for its future implementation in the Geoworks software so that Rock-Eval^® machines can routinely provide accurate SIC and SOC measurements.

1 Introduction

Soils contain large amounts of carbon stored as organic carbon and carbonates (Batjes, 1996; Zamanian et al., 2018). These carbon pools are dynamic and can contribute to climate regulation (Chenu et al., 2019; Zamanian et al., 2021). Moreover, organic carbon is also of primary importance in ensuring proper soil functioning (Hoffland et al., 2020). Global stocks of soil inorganic carbon (SIC) and soil organic carbon (SOC) are of comparable size (∼ 2500 Pg) when the soil is considered down to a depth of 2 m (Zamanian et al., 2021). Not all soils contain SIC, but the presence of carbonates is frequent, particularly in arid or semi-arid regions (Zamanian et al., 2018; Pfeiffer et al., 2023). The presence of SIC influences soil pH and therefore nutrient availability (e.g., Mkhonza et al., 2020). This explains the usual agricultural practice of liming to reduce soil acidity. In addition, recent studies have shown that the amount of SIC can vary significantly over short timescales (∼ 10 years) due to soil acidification resulting from certain farming practices such as nitrogen fertilization or irrigation (Zamanian et al., 2021; de Soto et al., 2017), suggesting that it may be worthwhile to assess the potential role of carbonates as a source or sink of C in the context of climate change. Moreover, land-based enhanced rock weathering is increasingly discussed as a strategy to contribute to atmospheric CO₂ removal. It consists of the spreading of silicate powder, whose weathering in soils leads to atmospheric CO₂ sequestration in carbonates (e.g., Kelland et al., 2020). Measuring the soil carbonates produced can provide information on the effectiveness of adding silicate minerals to sequester atmospheric CO₂ (Haque et al., 2020). As a result, the accurate quantification of soil organic and inorganic carbon (SOC and SIC, respectively) content is recognized as an increasingly important, yet not straightforward, task, especially in carbonated soils. There exist a wide range of methods for quantifying SOC and SIC; the two main approaches are dry and wet combustion methods, and they both come with their advantages and drawbacks.

Dry combustion methods involving CHN (carbon, hydrogen, nitrogen) elemental analyzers are increasingly used to measure total carbon. If the soil sample is not carbonated, the total carbon therefore represents only SOC. In carbonated soils, however, the amount of SIC is quantified using a separate analysis on a second subsample and the quantity of SOC is then determined as the difference between total carbon and SIC. SIC is generally measured using the amount of CO₂ produced by a known quantity of soil after acidification with HCl (Allison and Moodie, 1965). Alternatively, SOC can be directly quantified using CHN if carbonates have been previously removed using acid fumigation (Harris et al., 2001). The standard ISO 10694 method (ISO, 1995, 1999) describes how organic, inorganic, and total carbon can be measured in soil samples using dry methods. The main disadvantages of dry methods are that (1) SIC and SOC quantities are not measured on the same subsamples and (2) both SOC and SIC measurements can be inaccurate if the decarbonation is incomplete.

Wet methods are also largely used, especially when CHN analyzers are not available. The most commonly used wet oxidation method is the Walkley–Black method, which allows for measuring organic carbon directly in carbonated and non-carbonated soils (Walkley and Black, 1934). However, the Walkley–Black method has an average SOC yield of 77 %, and therefore the results have to be corrected to account for non-recovered SOC (Food and Agriculture Organization of the United Nations, 2020). The correction factor can be significantly different from one soil to another (Neal and Younglove, 1993), representing a limitation on the accuracy of wet oxidation methods (ISO 14235:1998 Soil Quality, 1998).

More recently, thermal analysis methods have been proposed as an interesting alternative to quantifying both SOC and SIC on a single subsample (Disnar et al., 2003; Vuong et al., 2016; Apesteguia et al., 2018). Among thermal methods, Rock-Eval^® (RE) thermal analysis has received particular attention in the recent years (Koorneef et al., 2023; Hazera et al., 2023), suggesting that if suitable corrections are applied, the Rock-Eval^® method could provide accurate quantification of both SIC and SOC. Moreover, a recent study has shown that the results provided by the Rock-Eval method are highly repeatable and reproducible, with relative errors on the order of 2 % to 4 % for the measurements of organic and inorganic C contents (Pacini et al., 2023). The aim of this technical note is to provide appropriately validated corrections to enable quantification of both SIC and SOC by Rock-Eval^® and to test the transferability of corrections between Rock-Eval^® machines.

2 Materials and methods

2.1 Soil samples

The study mobilized 240 soil samples taken from the surface layer (0–30 cm) of French agricultural soils. A stratification procedure was used to select the samples and ensure that each combination of texture class, SIC quantity, and SOC quantity was represented with at least 7 samples for low-carbonated soils and 14 samples for carbonated soils. As a result, the selected samples cover a wide range of particle-size distribution and SIC and SOC content, as shown in Fig. 1. SIC and SOC quantification was conducted at the Laboratory of Soil Analyses (Laboratoire d'Analyses des Sols, LAS) of Saint-Laurent-Blangy, accredited by Cofrac (French accreditation committee). LAS is a public laboratory operated by INRAE (French National Research Institute for Agriculture, Food and Environment) that analyzes samples for the French and European research community. At LAS, the total carbon is determined using a CHN analyzer (Thermo Fisher FLASH 2000) and inorganic carbon is quantified using a Bernard calcimeter. In carbonated samples, organic carbon is then determined as the difference between total carbon and inorganic carbon. SIC and SOC measured at LAS can be considered highly accurate and are therefore suitable to constitute reference measurements. Their associated uncertainties are calculated as

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\text{uncertainty SOC}& =\mathrm{0.02}\times \left(\left[\text{SOC}\right]+\mathrm{0.12}\times \left[{\text{CaCO}}_{\mathrm{3}}\right]\right)\\ & +\mathrm{0.49}+\mathrm{0.12}\times (\mathrm{0.016}\times \left[{\text{CaCO}}_{\mathrm{3}}\right]\\ & +\mathrm{0.63}),\mathrm{in}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}},\end{array}\end{array}$$

and

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\text{uncertainty SIC}& =\mathrm{0.12}\times (\mathrm{0.016}\times \left[{\text{CaCO}}_{\mathrm{3}}\right]\\ & +\mathrm{0.63}),\mathrm{in}\mathrm{g}{\mathrm{kg}}^{-\mathrm{1}}.\end{array}\end{array}$$

Figure 1 Distributions of SOC (a) and SIC (b) and particle size as a proportion of sand, silt, and clay (c), in the dataset of 240 samples, as measured by LAS, Saint-Laurent-Blangy. The SIC range is 0.5 to 96.8 g C kg⁻¹, and the SOC range is 5.37 to 44.64 g C kg⁻¹.

2.1.1 Rock-Eval^® measurements

Each of the 240 soil samples was ground using a Fritsch PULVERISETTE 6 tungsten carbide planetary mill and divided into three subsamples. One subsample was sent to LAS, Saint-Laurent-Blangy, for SIC and SOC determination. The two other subsamples were analyzed using Rock-Eval^® at ISTeP (Earth Science Institute, CNRS, Sorbonne University, Paris, France), and at Vinci Technologies (Nanterre, France). At ISTeP, samples were analyzed using a Rock-Eval^® 6 Turbo (RE6 Turbo), whereas at Vinci Technologies samples were analyzed using an RE6 standard instrument and a Rock-Eval^® 7 (RE7) instrument. The RE7 machine is the newest version of Rock-Eval^® machines, and we therefore considered it useful to check whether a correction designed on RE6 can be transferred to the new RE7 version.

The Rock-Eval^® thermal analysis consists of a pyrolysis step followed by an oxidation step. During pyrolysis, the sample is gradually heated in an inert environment (N₂) and evolved hydrocarbons, CO₂, and CO effluents are continuously monitored with time. Then the sample is cooled down before the start of the oxidation stage. During oxidation, the sample is gradually heated in purified laboratory air and the quantity of evolved CO₂ and CO is continuously monitored. The analytical setup used was the “SOIL” routine. The pyrolysis step started with a 3 min isotherm at 200 °C followed by a heating ramp of 30 °C min⁻¹ from 200 to 650 °C (Disnar et al., 2003). The heating routine for the oxidation step started with a 1 min isotherm at 300 °C followed by a heating ramp of 20 °C min⁻¹ until 850 °C and terminated with a 5 min isotherm at 850 °C (Baudin et al., 2015, adapted from Behar et al., 2001). The correct determination of SIC using the original settings of the widely used SOIL routine, even for carbonate-rich samples (see for instance Delahaie et al., 2023), suggests that there is no need to modify the Rock-Eval^® acquisition parameters for carbonated soils as suggested by Hazera et al. (2023). Pacini et al. (2023) detail the RE6 thermal analysis, as well as the usage of the SOIL routine in Geoworks.

All samples (ca. 60 mg) on all machines were analyzed using the SOIL routine. A reference sample (IFPEN_160000) is measured regularly in Rock-Eval® batch analyses for quality control. The target values are 3.28 ± 0.06 wt % and 3.26 ± 0.07 wt % for SOC and SIC, respectively. The machine is recalibrated when SOC and SIC values are outside this narrow range, which means that the acceptable measurement error is ca. 2 %. The thermograms were integrated to provide RE parameters and indices using the Geoworks v1.7 software. Notably, TOCre6 and MinC (corresponding to Rock-Eval^® SOC and SIC measurements) were provided by Geoworks. The thermograms showed that the carbonate soil samples contained neither dolomite nor siderite.

2.2 Design of the correction methods

It has long been known that TOCre6 and MinC tend to slightly (by ca. 10 %) underestimate SOC and overestimate SIC, respectively (Disnar et al., 2003). The data we collected from LAS, Saint-Laurent-Blangy, and the Rock-Eval^® machines (ISTeP and Vinci Technologies) also respect this general tendency, as shown in Fig. 2. The goal of this study was to design and validate a correction method for soils that aligns the TOCre6 and MinC values measured by Rock-Eval^® as closely as possible with reference SOC and SIC measurements provided by LAS. To do so, we considered several correction methods, each one implementing a different machine-learning (ML) model. We began by learning a correction method using the data obtained with the RE6 Turbo machine from ISTeP and pairing it with the soil organic and inorganic carbon data from LAS, Saint-Laurent-Blangy. We then tested the generalizability of the correction by applying it to the data obtained using other machines, i.e., the RE6 and RE7 data from Vinci Technologies. In the interest of brevity and legibility, certain modeling choices and techniques are omitted from the main text and are instead included in the Supplement in the form of a Jupyter Notebook (Kluyver et al., 2016) using open-source Python 3 code (Van Rossum and Drake, 2009).

Figure 2 Comparison of the SIC and SOC measurements before correction from the RE6 analysis by ISTeP and reference measurements by LAS, Saint-Laurent-Blangy. The red line represents the 1:1 diagonal.

After having tested several different approaches, the corrective model we propose here is a support vector machine (SVM) (Cortes and Vapnik, 1995) regression as implemented in the scikit-learn Python library (Pedregosa et al., 2011). SVMs are a collection of supervised ML models that are well suited to our usage as they adapt to both linear and non-linear data. SVMs using the RBF (radial basis function) kernel offer two different regularization hyper-parameters, which allows us to avoid overfitting the training data. Other models we have tested are a regular least-squares linear regression (LR), a constrained ridge regression (Ridge), and a random forest (RF) model. The LR model is an ordinary-least-squares linear regression that minimizes the residual sum of the differences between the observed and the predicted data. Ridge is a special case of LR where regularization is added to the coefficients using the ℓ² norm so as to avoid overfitting the training data. Lastly, RF is an ensemble modeling method that uses averaging of a collection of randomized decision trees, i.e., a random subset of features and of samples, to predict the target variable. Like SVMs, RFs are capable of capturing linear and non-linear relationships in the data. After correcting the SIC and SOC, it can happen that the corrected value is a small negative number close to zero, especially for low-carbonated soil samples. In this case, we simply set these negative SIC and SOC values to zero.

The dataset contains all the RE parameters (such as the hydrogen index (HI), the oxygen index (OI), pyrolyzable carbon (PC), or temperature stability parameters; see Pacini et al., 2023, and Baudin, 2024, for an exhaustive list) for the 240 soil samples including TOCre6 (namely SOC_RE6) and MinC (namely SIC_RE6), as well as SOC and SIC provided by LAS.

We develop two separate models for predicting the required correction of the SOC and the SIC. Our target variable is the difference between the SOC (or SIC) measured by LAS and the RE6 analysis:

$$\begin{array}{}\text{(3)}& {\displaystyle }{y}_{\text{SOC}}={\text{SOC}}_{\text{CHN}}-{\text{SOC}}_{\text{RE6}},\text{(4)}& {\displaystyle }{y}_{\text{SIC}}={\text{SIC}}_{\text{CHN}}-{\text{SIC}}_{\text{RE6}}.\end{array}$$

All the models are cross-validated during the training phase using two-thirds of the available data and then tested on the remaining one-third of samples. After a correlation analysis between the target variables and the RE6 parameters, as well as running several multivariable models, we concluded that the Rock-Eval^® SOC (or SIC) alone is sufficient to accurately estimate the required correction. This procedure began by calculating the Spearman correlation, whereby it was evident that the SOC_RE6 and SIC_RE6 are the most highly correlated with the target variables. Next, we trained the four ML models mentioned above (SVM, LR, Ridge, RF) using all RE parameters or subsets of them that were highly correlated with the target. We used two different thresholds of 0.3 and 0.5 in Spearman correlation for detecting the highly correlated parameters, allowing us to have a smaller or larger set of parameters. We also removed any parameters that were highly correlated among each other, this time using a threshold of 0.9 in terms of Spearman correlation. The results of all these tests were overwhelmingly in favor of using only the SIC_RE6 and SOC_RE6 in the correction models. Therefore, in the remainder of this note, all the proposed corrections are single-feature models. Finally, the corrected SIC and SOC of RE can be calculated as

$$\begin{array}{}\text{(5)}& {\text{SOC}}_{\text{corrected}}={\text{SOC}}_{\text{RE6}}+{\widehat{y}}_{\text{SOC}}\end{array}$$

and

$$\begin{array}{}\text{(6)}& {\text{SIC}}_{\text{corrected}}={\text{SIC}}_{\text{RE6}}+{\widehat{y}}_{\text{SIC}},\end{array}$$

where ${\widehat{y}}_{\text{SOC}}$ and ${\widehat{y}}_{\text{SIC}}$ are the predicted corrections to the ML model.

Figure 3 Correction accuracy comparison of the four tested ML models, the Disnar et al. (2003) corrections Disnar_BRS and Disnar_soil, and the original data before correction.

Figure 4 Precision of the proposed SVM correction applied to the test dataset.

For the sake of comparison, we also implemented the corrective approach described by Hazera et al. (2023) derived from a seminal study by Disnar et al. (2003). In this study, two different corrections are proposed for SOC in calcareous samples, depending on the state of degradation of organic matter. The proposed correction for samples enriched in fresh organic matter is hereafter named Disnar_BRS for biopolymer-rich samples, whereas the correction for samples containing more processed organic matter is hereafter named Disnar_soil. As the degradation state of soil organic matter sensu (Disnar et al., 2003) is not easy to determine, we decided to consider both corrections in our study for carbonated and non-carbonated samples. The correction for carbonated biopolymer-rich samples (Disnar_BRS) is

$$\begin{array}{}\text{(7)}& {\text{SOC}}_{\text{corrected}}=\mathrm{1.17}\times {\text{SOC}}_{\text{RE6}}\end{array}$$

and

$$\begin{array}{}\text{(8)}& {\text{SIC}}_{\text{corrected}}={\text{SIC}}_{\text{RE6}}-\mathrm{0.092}\times {\text{SOC}}_{\text{RE6}}.\end{array}$$

For non-carbonated samples, the same authors propose the following correction (Disnar_soil):

$$\begin{array}{}\text{(9)}& {\text{SOC}}_{\text{corrected}}=({\text{SOC}}_{\text{RE6}}+{\text{SIC}}_{\text{RE6}})\times \mathrm{1.068}\end{array}$$

and

$$\begin{array}{}\text{(10)}& {\text{SIC}}_{\text{corrected}}=\mathrm{0}.\end{array}$$

For non-BRS soils, they propose the following correction, Disnar_soil, for carbonated soils:

$$\begin{array}{}\text{(11)}& {\text{SOC}}_{\text{corrected}}={\text{SOC}}_{\text{RE6}}+\mathrm{1.092}\times {\text{SIC}}_{\text{RE6}}\end{array}$$

and

$$\begin{array}{}\text{(12)}& {\text{SIC}}_{\text{corrected}}={\text{SIC}}_{\text{RE6}}-\mathrm{0.092}\times {\text{SIC}}_{\text{RE6}}.\end{array}$$

They propose the following for non-carbonated soils:

$$\begin{array}{}\text{(13)}& {\text{SOC}}_{\text{corrected}}={\text{SOC}}_{\text{RE6}}+{\text{SIC}}_{\text{RE6}}\end{array}$$

and

$$\begin{array}{}\text{(14)}& {\text{SIC}}_{\text{corrected}}=\mathrm{0}.\end{array}$$

As the Disnar et al. (2003) model divides the soil samples into carbonated and non-carbonated soils, we also divide our data using a similar threshold of 2 g C kg⁻¹ of SIC. This approach results in three different datasets: the non-carbonated (80 samples), carbonated (160 samples), and all data (240 samples). In the accompanying Notebook, we have detailed a second threshold used to distinguish carbonated and non-carbonated samples based on the contribution to the S5 signal in the RE6 analysis. Due to the similarity of the two approaches, we have only presented the former in the following text of this publication.

Figure 5 Results of applying the proposed SVM correction to data obtained using Vinci Technologies' RE6 and RE7 machines.

2.3 Metrics to evaluate the correction methods

We use several different metrics to quantify the quality of the provided correction. The R² coefficient of determination is the proportion of the variable in the dependent variable that can be explained by the independent variable. Its upper bound and best score is 1, while its lower bound is minus infinity. An R² score of 0 means that none of the dependent variable's variability is explained and the model only predicts the mean value. It is calculated as

$$\begin{array}{}\text{(15)}& R^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{\sum {\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{\sum {\left(y_i-\stackrel{\mathrm{‾}}{y}\right)}^{\mathrm{2}}}},\end{array}$$

where ${\widehat{y}}_i$ is the predicted value of the ith sample and $\stackrel{\mathrm{‾}}{y}$ is the average value over all samples. The root mean square error (RMSE) and relative RMSE (RRMSE) are calculated as follows:

$$\begin{array}{}\text{(16)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\sum {\left(y_i-{\widehat{y}}_i\right)}^{\mathrm{2}}}{N}}},\text{(17)}& {\displaystyle }\text{RRMSE}={\displaystyle \frac{\text{RMSE}}{\stackrel{\mathrm{‾}}{y}}},\end{array}$$

where N is the number of samples. We also calculate the bias as the difference between the mean of the measured values and the mean of the predicted values. The bias allows us to know if the model tends to systematically over- or underestimate the target variable.

3 Results and discussion

3.1 Comparison of performances of the corrective methods

Having introduced several models for both the SIC and the SOC, we proceed with assessing their corrective performance using the test dataset. The metric we are mainly interested in is the RMSE, and to a lesser extent we are interested in the bias, the R² coefficient of determination, and the ordinary-least-squares linear equation describing the relationship between the reference CHN values provided by LAS and the (corrected) RE values.

Figure 3 shows the RMSE for the original data and for the corrected data after using one of the four different tested ML models and the correction factors provided by Disnar et al. (2003). Concerning the SIC correction, the overall best performance is achieved by the SVM model, as it is the only one to significantly reduce the error in all three datasets (carbonated, non-carbonated, and all samples). With the exception of the RF model for non-carbonated data, all the models we propose offer a considerable decrease in RMSE. The RF model, due to its intrinsic complexity and bootstrapping, is probably more impacted than other ML models by the smaller dataset size, as only the 80 non-carbonated samples were considered and one-third of those were set aside for testing. Regarding the SOC correction, the models we propose strongly increase the similarity between LAS and Rock-Eval^® ISTeP SOC values. They have virtually equal performance, resulting in a 2- to 3-fold decrease in RMSE. A more detailed comparison of the behavior of the models as a function of the input data is available in the Supplement.

The correction factors proposed by Disnar et al. (2003) significantly decreased the RMSE for SIC and SOC for all datasets. However, Fig. 3 clearly shows that this correction is not optimal as the RMSE is on average about two times higher compared to the ML-based corrections. This is not very surprising as, contrary to the other methods, the correction proposed in Disnar et al. (2003) is not derived from an optimized model trained and tested on separate data. Even though the four ML-based corrections render highly similar performance, the overall best-performing model is SVM. A paired t-test comparison between the Disnar and SVM models is available in the Notebook and shows the two corrections provide significantly different results. The results of the SVM correction are shown in Fig. 4. The RMSE is divided by a factor of more than 3 for the SIC and more than 2 for the SOC with respect to the raw data before correction. For both the SIC and the SOC, the equation describing the relationship between corrected RE6 and LAS data is almost identical to the 1:1 ratio, showing that our proposed correction closes the gap between the two methods.

3.2 Transferability of the corrective method to other Rock-Eval^® machines

In order to test the transferability of our results to other Rock-Eval^® machines, we apply our SVM model learned on ISTeP's RE6 Turbo machine to the data obtained using Vinci Technologies' RE6 and RE7 machines. The reference SIC and SOC data stay the same, i.e., provided by LAS, Saint-Laurent-Blangy. Figure 5 shows a comparison for the RE6 (a and b) and RE7 data (c and d), respectively. The correction quality is comparable to that of the ISTeP data, with, not surprisingly, only slightly higher RMSEs. The overall bias is not only of the same order but also sometimes lower, in part due to the larger sample size as the original correction metrics are calculated on only one-third of the dataset used as testing data. The possibility of transferring a correction model for SIC and SOC learned on one Rock-Eval^® machine to another was somehow expected, as Pacini et al. (2023) observed that TOCre6 and MinC are very reproducible across different RE6 instruments and for RE6 and RE7 instruments.

4 Conclusions

Accurately estimating the quantity of soil organic carbon has important implications for the monitoring of soil health as well as climate regulation techniques and policies. The classic dry combustion and wet methods, though standardized decades ago, come with a series of drawbacks when considering carbonated soils. Previous studies have suggested that, on top of providing information on soil organic matter biogeochemical stability (Barré et al., 2016), thermal analyses could be an accurate means of determining SOC and SIC quantities in carbonated soil with lower experimental uncertainty. Our work confirms this hypothesis and proposes the first validated correction method to accurately determine the quantities of SOC and SIC from a Rock-Eval^® thermal analysis. This correction method, based on SVM machine learning, can be transferred to different RE6 and RE7 instruments and will be implemented in the Geoworks software so that Rock-Eval^® machines can routinely provide accurate SIC and SOC measurements, at least in the SOC and SIC value ranges investigated in this study (up to ca. 5 % and ca. 8 % for SOC and SIC, respectively).