Estimating soil carbon sequestration potential with mid-IR spectroscopy and explainable machine learning

Soil carbon sequestration refers to the process of capturing atmospheric carbon through plant photosynthesis and storing it in soil as organic carbon. The primary mechanism for carbon sequestration is the adsorption of organic carbon molecules onto the mineral surfaces of the soil's fine fraction (clay + silt ≤ 20 µm), forming mineral-associated organic carbon (MAOC). Soil has a finite capacity to stabilise and sequester organic carbon, known as carbon saturation capacity, which depends on the proportion of reactive minerals in the soil. The difference between the current MAOC content and the carbon saturation capacity is referred to as the organic carbon saturation deficit (C_def) or sequestration potential. Fourier-transformed (FTIR) mid-infrared (mid-IR) spectroscopy can simultaneously measure soil properties relevant to carbon stabilisation: organic carbon functional groups, clay and iron-oxide mineralogy and particle size. Therefore, we hypothesise that mid-IR spectroscopy can effectively and accurately estimate C_def. Here, we aim to (i) develop spectroscopic models to estimate the MAOC and C_def of 482 Australian topsoil samples, (ii) model MAOC and C_def using mid-IR spectra and an interpretable machine learning algorithm, and (iii) further interpret the MAOC and C_def models using SHapley Additive exPlanations (SHAP). Using frontier line analysis, we fitted a function to the upper envelope of the MAOC vs. clay + silt relationship to derive C_def. We recorded mid-IR spectra of the samples and used the regression trees method CUBIST to model MAOC content and C_def. We interpreted these models by examining the regression trees and using SHAP. The models were unbiased and estimated MAOC content with R² of 0.86 and RMSE of 2.77 (g kg soil⁻¹), and C_def with R² of 0.89 and RMSE of 3.72 (g kg soil⁻¹). Model interpretation showed that C_def estimates relied on negative interactions with absorptions from organic matter functional groups and positive interactions with absorptions from clay minerals. Our results demonstrate that mid-IR spectra can effectively estimate MAOC and soil C_def, providing a rapid, cost-effective method for assessing and monitoring this critical soil function.

1 Introduction

Soil organic carbon (C) sequestration refers to the process by which plants capture atmospheric C through photosynthesis and store it in the soil. The United Nations Framework Convention on Climate Change (UNFCCC) has identified soil C sequestration as a critical, nature-based process for withdrawing atmospheric carbon dioxide (CO₂) (UNFCCC, 2019). Soil organic C sequestration also improves soil health, food and nutritional security, water quality, biodiversity, and elemental recycling (Lal et al., 2015). Thus, it is crucial to estimate the amount of C stored in soil and how much it could store in the future to advance our scientific understanding of C cycling. This understanding will provide the foundation for land managers to develop practices that enhance C sequestration and for policymakers to formulate climate change adaptation strategies. However, rapidly, cost-effectively, and scientifically estimating the soil C saturation deficit remains challenging.

Soil C from plants begins as particulate organic C (POC). Over time, soil microorganisms consume this POC, and some of it is broken down into smaller molecules. Some of these molecules are protected from further decomposition through adsorption onto mineral particles, forming mineral-associated organic carbon (MAOC) and providing protection within soil microaggregates (Hassink and Whitmore, 1997; Six et al., 2002; Beare et al., 2014). Soils with higher silt and clay content have a larger mineral surface area and a greater capacity to adsorb and stabilise C. Hassink (1997) found a positive linear relationship between the proportion of clay and silt (particles <20 µm) and the amount of C in this fraction. This relationship has been used to estimate the soil's maximum capacity to stabilise C (Hassink, 1997), referred to as the C saturation capacity (C_sat). The difference between actual MAOC content and C_sat is known as the C saturation deficit (C_def) or C sequestration potential.

Subsequent studies, such as Six et al. (2002), also found a direct relationship between MAOC and the amount of clay and silt in soil, further recognising that this relationship depends on the reactivity of the soil's clay minerals. Many researchers have since used such linear relationships to estimate C_sat. Feng et al. (2013) found this approach underestimated C_sat and proposed a boundary line method as an alternative, fitting a line to the upper tenth percentile of the data in the MAOC vs. clay and silt relationship. More recently, the relationship has been fitted using quantile regression at the 95th percentile of the data (Georgiou et al., 2022). However, these methods underestimate C_sat by fitting through the data rather than capturing the maximum values. To address this shortcoming, Viscarra Rossel et al. (2024) proposed using a bootstrapped frontier lines analysis that fits an envelope to the maximum values of the relationship between MAOC and the soil's fine fraction, thereby preventing underestimation of the soil's C storage capacity and providing uncertainty estimates. Additionally, considering the maximum attainable C storage under a given environment (C_Amax) approaches maximum asymptotically (Ingram and Fernandes, 2001; Viscarra Rossel et al., 2024), the frontier line approach better reflects C_Amax by showing the asymptotic increase in soil C storage capacity with increasing soil clay and silt content (Viscarra Rossel et al., 2024).

Establishing reliable estimates of C_Amax and C_def using the methods above requires many soil samples with measured MAOC and clay-plus-silt content. Measuring MAOC involves fractionating soil to isolate the C in the ≤50 µm soil fraction and then measuring the organic C content (Poeplau et al., 2018). Fourier-transformed (FTIR) mid-infrared (mid-IR) spectroscopy offers a faster, more cost-effective, and repeatable alternative. It measures soil composition by capturing interactions between mid-IR wavelengths and the vibrations of bonds in soil molecules, providing data on a soil's organic and mineral composition (Viscarra Rossel et al., 2022). These spectra have been used to estimate organic and inorganic C, clay, sand and silt contents, cation exchange capacity and other chemical, physical and biological properties through calibration that relates the measured soil properties to their spectra (Soriano-Disla et al., 2014).

Mid-IR spectra serve as an integrative “molecular fingerprint” of the soil, reflecting its mineralogy, organic matter, and physical properties (Viscarra Rossel et al., 2006), which directly determine a soil's biological activity, soil structure and ultimately the ability to sequester C (Wiesmeier et al., 2019). Baldock et al. (2019) estimated the C saturation deficit (C_def) of New Zealand soils using pedotransfer functions derived from the quantile regression approach, modelling C_def with mid-IR spectra through partial least squares regression (PLSR), showing good predictability. Similarly, Karunaratne et al. (2024) estimated the C_def of Australian soils using a quantile regression approach and modelled it with mid-IR spectra coupled with PLSR, also achieving good predictability. We did not find other research that estimates C_def using soil spectra. We hypothesise that mid-IR spectra, combined with explainable machine learning, can be used to estimate soil MAOC content and C_def while also providing insights into how the model uses spectral absorption features to identify the soil constituents important for prediction. Thus, we aimed to:

1. Develop spectroscopic models to estimate the MAOC content and the C_def of Australian soils using mid-IR spectra with the regression trees algorithm CUBIST;

2. Interpret these models by analysing the CUBIST rulesets and SHapley Additive exPlanations (SHAP) values to understand how the absorptions of soil organic and inorganic constituents affected model prediction.

2 Methods

2.1 Soil samples

We used 488 topsoil samples from 275 sites across Australia (Fig. 1). The soils were sampled from three depth layers (0–10, 10–20 and 20–30 cm). All soil orders in the Australian soil classification were present, except Anthroposol and Organosol (Teng et al., 2018). Kandosols were the most abundant soil type, followed by Tenosols and Calcarosols, Chromosols and Vertosols, while Rudosols, Dermodols, Kurosols, Ferrosols, and Podosols were present in smaller numbers. Three Hydrosols were excluded from further analysis due to the distinct C storage mechanisms in anoxic soils (Six et al., 2024).

Figure 1 Location of sampling points.

The sampling area spans the main Köppen-Geiger climate zones (Beck et al., 2018), with most of the samples collected from arid hot deserts, with smaller proportions from arid hot steppes and tropical savannahs. Samples were primarily collected from areas with minimal human impact, particularly nature conservation sites, native vegetation grazing lands, and other minimally used areas. Only a small proportion of samples came from production or intensive land use. The vegetation at the sampling sites was diverse, comprising 24 major vegetation groups, with eucalyptus woodlands being the most common (Commonwealth of Australia, 2020). Most samples were taken from native vegetation or natural bare land, with the rest from non-native vegetation or cleared land (ABARES, 2022).

2.2 Soil fractionation

Soil samples were fractionated through physical granulometric separation. The samples were dispersed in deionised water using an ultrasonic probe (Sonics VCX 500 Sonicator, Newtown, Connecticut) with an energy output of 500 J mL⁻¹ for 200 s (Walden et al., 2025). After dispersion, the samples were fractionated using an automated wet sieving apparatus (Analysette 3 Pro, Fritsch GmbH, IdarOberstein, Germany) with 250 and 50 µm sieves. The resulting soils were in three size fractions: macroaggregates (2000–250 µm), microaggregates (250–50 µm), and the fine fraction (≤ 50 µm). The fractionated samples were then oven-dried at 60 °C overnight and ground to approximately ≤ 80 µm, before the organic C content of each size fraction was measured using an elemental analyser (SoliTOC Cube, Elementar Analysensysteme, Hanau, Germany). The organic C content of the fine fraction, representing MAOC, was recorded in grams per kilogram of whole soil.

2.3 mid-IR spectroscopy

The whole soils (sieved to ≤ 2 mm) were air-dried before fine grinding to $\approx <\mathrm{80}$ µm. The mid-IR spectra of the finely ground samples were measured with a diffuse reflectance infrared Fourier transform (DRIFT) spectrometer (Bruker Invenio HTS-XT, Massachusetts, United States). Spectra were recorded from 4000–450 cm⁻¹ with a spectral resolution of 4 cm⁻¹ and measuring 64 scans per sample. The spectrometer was calibrated with a gold standard before measuring each sample plate with 23 samples (Bruker, Massachusetts, United States). Reflectance spectra were recorded in $\log \frac{\mathrm{1}}{R}$ (apparent absorbance).

Silt + clay content

The silt and clay content of the whole soil was determined using mid-IR spectroscopic modelling with CUBIST (Hicks et al., 2015). The silt % model has an R² value of 0.84 with a concordance of 0.92, and the clay % model has an R² value of 0.90 with a concordance of 0.95. The estimated silt and clay content in % was combined for further analysis.

2.4 Frontier lines and calculation of C_def

The MAOC content of samples displayed a log-normal distribution. We performed a log_e transformation on the MAOC content and removed three outliers that were more than 1.5 times the interquartile range above Q3 or below Q1. We proceeded with the analysis of the remaining 482 samples from 270 sites.

We fitted a monotonically increasing and concave frontier line (Parmeter and Racine, 2013) to the relationship between log(MAOC) and clay + silt content of the samples using the smooth, non-parametric frontier line analysis with the R package SNFA (McKenzie, 2022). We calculated the C_Amax and C_def following the approach described in Viscarra Rossel et al. (2024). Each point on the frontier line represents the maximum attainable amount of MAOC that soil could store for a particular clay and silt content.

To estimate uncertainty, we performed 100 non-parametric bootstrap resamples to fit the frontier lines, keeping samples from the same site together during resampling to prevent data leakage. We then averaged all 100 frontier-line fits from the bootstraps. The C_def was calculated as the difference between the estimated mean frontier line and the MAOC content. We also computed the uncertainties of our frontier line estimate by calculating the 95 % confidence limits. All values were then back-transformed to their original units for the spectroscopic modelling.

2.5 Spectroscopic modelling

The mid-IR spectra were interpolated to 32 cm⁻¹ wavenumber intervals to reduce inherent collinearity. Since mid-IR spectra are highly collinear and contain broad absorption features, we interpolated the spectra to 32 cm⁻¹ to reduce the redundant information passed into the machine learning model (Deiss et al., 2020). Visual checks confirmed relevant absorption features remained distinguishable at this resolution. We also checked C_def model performance using spectra interpolated to 8, 16, 24, and 32 cm⁻¹ resolutions and found no significant difference between these resolutions. Preprocessing consisted of an initial offset correction, in which the minimum spectral value minus 0.01 was subtracted from all measurements so that each spectrum was shifted to a common baseline just above zero, followed by a standard normal variate (SNV) transformation, and a final offset correction to address the baseline shift introduced by the SNV transformation. Spectral regions that were either featureless (4000 to 3746 cm⁻¹) or containing distracting features from noise and artefacts from water and CO₂ (2370 to 2082 cm⁻¹) were removed before modelling.

We modelled the MAOC and the estimated C_def with CUBIST. CUBIST is a rule-based regression tree algorithm (Quinlan, 1992; Wang and Witten, 1997). CUBIST creates a tree structure, with branches as a series of “if-then” conditions, then reduced into rules. Each CUBIST rule corresponds to a subset of the data that satisfies the rule's condition. For each rule, a linear regression model is fit to the data using relevant predictors (Kuhn et al., 2012). CUBIST balances accurate predictions and model interpretability through its rule-based structure. CUBIST is tuned by two parameters: committees and neighbours. The number of committees specifies the number of ensembles contributing to the final prediction, with more committees typically improving performance but reducing interpretability, and the number of neighbours specifies how many nearest-neighbours of a sample CUBIST uses to adjust its rule-based predictions. Viscarra Rossel and Webster (2012) described the method for spectroscopic modelling. In our experiments, since our goal was to understand which spectral regions influence predictions and how they relate to soil properties, we prioritised model interpretability by using a single committee to maintain model transparency, avoiding the added complexity of ensemble averaging. We optimised the number of neighbours by testing all values from 0 to 9. Model fitting and validation were carried out using 10-fold cross-validation grouped by site, in which the 270 sampling sites were randomly assigned to 10 folds to ensure that samples from the same site and the three depth layers were kept together within the same fold. We assessed the models based on their coefficient of determination (R²), Lin's concordance correlation coefficient (CCC) (Lin, 1989) and the root mean squared error (RMSE).

We propagated the uncertainty of the frontier line fitting and the CUBIST modelling. From the 100 frontier line fits made with the bootstraps, we derived the upper and lower 95 % confidence intervals (CI) for the frontier line fit and calculated the upper and lower limit of C_def. The upper and lower limits of C_def were also modelled with CUBIST following the same method described above.

2.6 Interpretation

To interpret the models, we extracted each CUBIST rule from the MAOC and C_def models to analyse their rule partitioning. For the MAOC model, we examined the distribution of MAOC values within each rule, while for the C_def model, we analysed the distributions of both MAOC and C_def values within each rule. For the linear models in each CUBIST rule, we examined the wavenumber corresponding to specific absorptions of soil constituents and their coefficients. Furthermore, we calculated the SHAP (SHapley Additive exPlanations) values for each sample for each linear model of the CUBIST rules, focusing the analysis on the C_def model. (SHAP analysis results of MAOC CUBIST model provided in the Supplement) The SHAP values are used to explain the outputs of machine learning models. SHAP is based on game theory (Shapley, 1953) and assigns an importance value to each instance at each feature (in our case, each sample's absorptions at specific wavenumbers) in a model. While the regression coefficients summarise the average effect of a wavenumber within a given rule, SHAP values provide instance-level attributions that quantify each wavenumber's contribution to the prediction for each individual sample. Positive SHAP values indicate a positive impact on the prediction, while those with negative values indicate a negative impact. The magnitude measures the strength of the effect.

All statistical analyses were performed using R (R Core Team, 2024).

3 Results

3.1 The maximum attainable MAOC storage, the MAOC deficit and C sequestration potential

Our samples represent a wide geographical area in Australia (Fig. 1) with large variations in MAOC content and texture (Table 1). The MAOC content ranges from 0.27 to 50.04 g kg soil⁻¹, while silt content ranges from 0.54 % to 31.81 %, and clay content ranges from 2.34 % to 54.25 % (Table 1). The frontier line estimates the maximum C that can be stored in their current environments over their range of clay + silt contents for all 482 samples, with their 95 % confidence intervals shown in Fig. 2. The frontier line increases with increasing clay + silt content to around 20 %–45 %, after which the rate of increase slows. The C_Amax ranges from 5.29 to 45.79 g kg soil⁻¹ with a mean of 32.76 g kg soil⁻¹ (Table 1). The C_def ranges from none to 45.17 g kg soil⁻¹ with a mean of 26.31 g kg soil⁻¹ (Table 1).

Table 1 Summary statistics.

Note: SD = Standard Deviation, Min = Minimum, Q_0.25 = Lower 25 % quartiles, Med = Median, Q_0.75 = Upper 25 % quartiles, Max = Maximum, Skew = Skewness.

Figure 2 Frontier lines and its 95 % confidence interval fitted using all 482 samples.

3.2 Spectroscopic modelling of MAOC content

The CUBIST model predicts MAOC with an RMSE of 2.77 g kg soil⁻¹, is unbiased with R² of 0.86, and CCC of 0.91 (Table 2, Fig. 3 b). The model partitions the data into four rule sets, corresponding to different MAOC content levels, which increase from Rule 1 to Rule 4 (Fig. 3a). Samples in Rule 1 have the least MAOC and are not significantly different from Rule 2 (Fig. 3a). Rule 3 samples have significantly more MAOC than Rule 1 but are not significantly different from Rule 2 (Fig. 3a). Rule 4 samples have significantly more MAOC than all other rules and exhibit the largest spread (Fig. 3a).

Table 2 Tuning parameters and model statistics for MAOC and C_def CUBIST models.

Note: RMSE = Root mean square error, CCC = Lin's concordance correlation coefficient, CI = Confidence interval.

Figure 3 CUBIST model result for MAOC. (a) The distribution of MAOC content for each CUBIST rule and Tukey's HSD between each CUBIST rule. (b) The correlation between observed and predicted MAOC of the CUBIST model, coloured by CUBIST rules. (c) The coefficient of each linear model for each CUBIST rule is plotted over the mean spectra of each CUBIST rule.

The mean mid-IR spectra of the samples of the four rule sets show overall consistent patterns, with differences in absorption intensities at 3700–3500, 2946–2850, 1986–1794, and 1634–1300 cm⁻¹) (Fig. 3c).

Specifically, the mean spectrum of Rule 4 has the highest absorption in the 2946–2850 cm⁻¹ region associated with organic C (C–H vibrations of alkyl CH₂), corresponding to having the highest MAOC content (Fig. 3a, c).

The wavenumbers selected for linear models of the four rules differ, although there is some overlap. All rules use wavenumbers between 2946–2850 cm⁻¹, organic C–H vibrations of alkyl CH₂ groups (Nguyen et al., 1991), though the specific selections vary (Fig. 3c), suggesting that the models rely directly on spectral signals from organic C to predict MAOC content. Rule 1 exhibits densely distributed wavenumbers across both these regions with high coefficient values. Rule 3 shows a similarly dense distribution, concentrated primarily in the 2946–2850 cm⁻¹ region, with the largest coefficient values. Rule 2 displays more sparsely distributed wavenumbers across both regions, while Rule 4 uses only a few select wavenumbers around 2946–2850 cm⁻¹.

Rules 1, 2, and 3 with the smallest MAOC values all use the region between 1986–1794 cm⁻¹, associated with quartz, whereas Rule 4 does not (Fig. 3c). Quartz is chemically inert, carries negligible surface charge, and has a low specific surface area compared to clay minerals, which limits the reactive surface area available for organo-mineral bonding and thus associated with low MAOC content, as found in coarser-textured soils dominated by quartz. Rules 1, 2, and 3 also use region 2515 cm⁻¹ associated with carbonate (Nguyen et al., 1991), as soils with more carbonate commonly form in arid or semi-arid regions with low plant productivity and rainfall, and therefore low organic C.

Rule 4 uniquely includes absorptions at the 3750 cm⁻¹ region, associated with the hydroxyl stretching vibrations of clay minerals (Nguyen et al., 1991). Clay minerals provide a more mineralogically reactive soil matrix that facilitates greater organo-mineral bonding, and the model's use of spectral signals reflecting mineral surface reactivity is particularly informative in soils where mineral-organic associations are most developed. Rule 4 also uses wavenumbers between 1762–1634 cm⁻¹, associated with amide C=O bond (Volkov et al., 2021), as well as wavenumbers around 1154 cm⁻¹, which correspond to the SiO₂ lattice (Spitzer and Kleinman, 1961) and C-OH stretch of aliphatic O–H (Senesi et al., 2003) (Fig. 3c), with Rule 3 also using these latter wavenumbers to a lesser extent. SiO₂ lattice indicating the silicate mineral of the soil matrix, while aliphatic C-OH groups indicate polysaccharide-derived and carbohydrate-like organic matter, and amide indicates protein-derived organic matter. Taken together, this suggests that the model draws on both organic-matter and mineral-related absorptions to estimate MAOC content.

3.3 Spectroscopic modelling of the organic C deficit (C_def)

The model predicts C_def with an RMSE of 3.72 g kg soil⁻¹, R² of 0.89, and CCC of 0.94 while also being unbiased (Table 2, Fig. 4c). The model partitions the data into 3 rule sets, and the linear models of each CUBIST rule also show good precision (Table 3).

Table 3 Model statistics for each linear model of the CUBIST rules in the mean C_def CUBIST model.

Note: RMSE = Root mean square error, CCC = Lin's concordance correlation coefficient.

Figure 4 CUBIST model result for C_def, showing the CUBIST rules separation, including the distribution of (a) C_def and (b) MAOC content for each CUBIST rule and Tukey's HSD between each CUBIST rule. Along with (c) the correlation between observed and predicted C_def of the CUBIST model coloured by CUBIST rules, and (d) the coefficient of each linear model for each CUBIST rule plotted over the mean spectra of each CUBIST rule.

Rule 1 includes samples with the lowest C_def and the highest MAOC content, representing samples that have smaller C sequestration potential, as these samples contain more MAOC (Fig. 4a, b). Rule 2 represents samples with intermediate C_def, and contains little MAOC and clay and silt content, representing coarser-textured soils with more C sequestration potential than samples in Rule 1 because they hold less MAOC (Fig. 4a, b). Rule 3 includes samples with high C_def, low MAOC content and the most clay and silt content. Since these samples contain the finest particles, their capacity is largest and is thus undersaturated with C relative to their potential (Fig. 4a, b).

The three rule sets show similar overall mean spectral patterns but with distinct differences in absorption intensities at key regions, including 2946–2850 cm⁻¹ associated with organic C, 1986–1794 cm⁻¹ associated with SiO₂ overtone and combination bands, and 1538–1218 cm⁻¹) region associated with various organic and mineral absorptions (Fig. 4d). The wavenumbers selected for the models in each CUBIST rule are generally consistent, with the magnitude of the coefficient decreasing from Rule 1 to Rule 3 (Fig. 4d).

In the 2946–2850 cm⁻¹ region, associated with organic C–H vibrations of alkyl CH₂ groups (Nguyen et al., 1991), Rule 1 shows greater average absorption compared to Rule 2 and Rule 3 consistent with Rule 1 having the highest MAOC content (Fig. 4b, d). All three CUBIST rules use wavenumbers within and near this region with relatively large coefficients, but the coefficient magnitude decreases from Rule 1 to Rule 3 (Fig. 4b, d). This pattern reflects that the model predicts C_def by leveraging the spectral signal of existing organic C (MAOC) already occupying reactive mineral surfaces. As C_def represents the remaining sequestration potential, more existing MAOC implies less remaining capacity. Thus, Rule 1 with the highest MAOC and largest organic-region coefficient has the lowest C_def. Rule 2 has intermediate MAOC and coefficient magnitude, and Rule 3 has the lowest MAOC, smallest coefficient, and consequently the highest C_def (Fig. 4a, b, d).

In the region near 1986–1794 cm⁻¹, which is due to the overtones of Si-O vibrations (Volkov et al., 2021), absorption intensity decreases from Rule 2 to Rule 1 to Rule 3, corresponding to decreasing sand content and increasing clay and silt content (Fig. 4d).

All three rules have prominent absorption at and near 1634 cm⁻¹, which are associated with amide, carboxylate and carboxylic acid (Nguyen et al., 1991; Tanykova et al., 2021), aromatic –C=C– stretch (Du et al., 2014), HO–H stretch (Kronenberg, 1994), N–H bend, C=O stretch (Volkov et al., 2021) and absorbed water (Max and Chapados, 2009) (Fig. 4d). This indicates that organic matter types, including polysaccharide-derived, carbohydrate-like, and protein-derived organic C, are used across all three rules similarly to the 2946–2850 cm⁻¹ region, to predict whether the C saturation capacity is filled.

In the fingerprint region (1550–450 cm⁻¹), the band assignments are more challenging due to significant overlaps between mineral and organic absorptions (Soriano-Disla et al., 2014). The region from 1538 to 1218 cm⁻¹, likely associated with quartz minerals as well as organic matter (Volkov et al., 2021), is more prominent in Rule 2 and Rule 1, and lower in Rule 3 (Fig. 4d). Rule 3 exhibits proportionally larger coefficients for wavenumbers in the fingerprint region because of low organic C content and high fine mineral particle content (Fig. 4b, d).

The absorption near 2515 cm⁻¹ due to carbonates shows more prominent absorption in Rule 3. Where C_def is high, and MAOC is low, which matches the tendency of higher carbonate soils from arid or semi-arid regions with low organic C input that leave mineral surfaces unsaturated.

Figure 5 The correlation between observed and predicted C_def of the CUBIST model coloured by CUBIST rules, as well as the observed and predicted C_def estimated from the upper 95 % CI and lower 95 % CI of the frontier line fit. The grey envelopes represent the range of C_def predictions obtained when CUBIST is applied separately to the upper and lower 95 % CI frontier-line estimates, indicating the uncertainty of the frontier-line fit propagating to the C_def predictions.

The model statistics of the CUBIST models of C_def estimated from the upper and lower 95 % CI of C_Amax are shown in Table 2. The model for the C_def estimated with the lower 95 % CI of the frontier line performs better than the model estimated with the upper 95 % CI. This can be attributed to the upper 95 % CI of the frontier line being more uncertain than the lower 95 % CI. Specifically, the upper uncertainty of the frontier line fit is high around 25 % clay + silt content due to the low sample number (Fig. 2). The uncertainty of C_def estimated from CUBIST models of C_def calculated from the upper CI and lower CI of the C_Amax is shown in Fig. 5.

3.4 C_def model interpretation with SHAP

The SHAP contribution of spectral absorption at each wavenumber for the linear model of each CUBIST rule is shown in Fig. 6. The SHAP values coincide with the regression coefficients of the CUBIST rules (Fig. 6). The regression coefficients and SHAP values are generally consistent: large coefficients correspond to strong SHAP model contributions. Rule 1 shows strong contributions primarily from organic C features, and Rule 2 displays a similar pattern but with more contributions from the fingerprint region. For Rule 3, there is a relatively stronger contribution from the absorptions in the double bonds region (including absorption from quartz and the region associated with amide overlapping with other absorptions), and the fingerprint regions have a relatively stronger contribution (Fig. 6).

Figure 6 The mean spectra, key spectral assignment, and the SHAP contribution of the spectral regions used in each linear model of each C_def CUBIST rule. A positive SHAP value indicates a positive contribution to a model with increased absorbance, whereas a negative SHAP value indicates a negative contribution with increased absorbance. The magnitude of SHAP indicates the strength of the contribution. The SHAP values are plotted over the pre-processed spectra of each rule set. The SHAP values are coloured by the normalised absorbance value at each wavenumber, ranging from −1 (lowest absorbance at each wavenumber) to 1 (highest absorbance at each wavenumber). SHAP values of each rule are plotted in different y-axes to accommodate differences in magnitude across rules.

The SHAP values indicate positive and negative contributions from spectral regions associated with characteristic absorption of clay minerals, organic matter, and quartz (Fig. 6). Generally, peaks associated with organic C have a negative model contribution with an increase in absorbance, while the troughs have a positive contribution with increasing absorbance (Fig. 6). Absorbance in these regions indicates existing MAOC that already occupies reactive mineral surface sites. As mineral binding sites fill, the remaining deficit diminishes. Similarly, absorptions associated with clay minerals and silicate have a positive model contribution, while the troughs have a negative contribution (Fig. 6). Absorbance from clay minerals indicates abundant reactive surface area that is available but not yet occupied by organic matter. The positive SHAP contribution reflects unrealised adsorption capacity, whereas quartz has negligible reactive surface area and contributes to coarser texture without contributing to adsorption capacity, which constrains C_Amax.

In comparison to SHAP analysis of the MAOC CUBIST model (Fig. S1 in the Supplement), both models draw heavily on the organic C-H region, but the C_def model shows a progressive shift across rules from organic-C-dominated (Rule 1) to mineral-dominated (Rule 3), indicating the model increasingly relies on mineralogy for available surface in soil with higher C_def. The MAOC model's contributions remain more consistently concentrated in the organic C-H region across all rules, with the fingerprint region playing a lesser role throughout. This reflects that the MAOC model is more driven by the organic C information in the spectra, whereas the C_def model integrates both organic and mineral information.

4 Discussion

Our findings support the hypothesis that mid-IR spectra, combined with machine learning and enhanced by SHAP analysis for interpretability, can accurately estimate soil MAOC content and C_def (Table 2) by elucidating the contribution of specific mid-IR absorptions.

Our results demonstrate that combining soil spectroscopy with machine learning offers a rapid, cost-effective, and robust method for estimating MAOC and C_def. The spectroscopic approach enables many more measurements than conventional methods, thereby enhancing our understanding of how MAOC and C_def vary in soil across space and time (Angers et al., 2011). This approach could also provide essential data for soil biogeochemical and Earth System models, improving their initialisation, validation and ongoing development (Stewart et al., 2007; Georgiou et al., 2022; Abramoff et al., 2022; Vereecken et al., 2016). Given that C storage is a key soil function for maintaining soil health (Lal, 2016; Lehmann et al., 2020), our findings highlight how the current state and potential for C sequestration can be rapidly and cost-effectively measured as part of soil health assessment (Vogel et al., 2019). This aligns with growing evidence that soil spectra, when combined with machine learning, can model soil functions, going beyond predicting individual soil properties (Cohen et al., 2006; Elliott et al., 2007; Cécillon et al., 2009; Viscarra Rossel et al., 2010; Maynard and Johnson, 2018; Deiss et al., 2023).

Two other studies estimated soil C_def using mid-IR spectroscopic modelling Karunaratne et al. (2024); Baldock et al. (2019). Unlike these studies, which used quantile regressions to estimate C_def, our approach avoids under- or over-estimations (of C_def) using bootstrapped frontier lines that more accurately capture the relationship between MAOC and clay + silt content (Viscarra Rossel et al., 2024). Specifically, the frontier-line approach estimates the upper envelope of the MAOC–(clay + silt) relationship under current environmental conditions. Other approaches that aim to avoid under- or overestimation, e.g. Shi et al. (2025), which fit quantile regressions to mineralogically stratified subsets of the data, are effective but still impose a parametric relationship via an internal upper percentile. In contrast, the frontier-line prevents observed values from exceeding the estimated C_Amax and allows them to level off at high clay + silt contents, reflecting diminishing stabilisation due to finite organic inputs. As a result, frontier-line analysis reduces both underestimation (by targeting the upper boundary rather than an internal quantile) and overestimation (by avoiding unconstrained extrapolation at high clay + silt), yielding more realistic estimates of attainable MAOC storage. Additionally, unlike the earlier studies, we characterised two distinct sources of uncertainty: those from the frontier-line fitting and from the cross-validated CUBIST model. The 95 % confidence limits of the frontier-line fit were propagated to C_def predictions by applying CUBIST separately to the upper and lower CI of C_def estimation (Fig. 5).

The MAOC and C_def models relied on spectral regions related to organic functional groups such as the C-H groups near 2900 and 2800 cm⁻¹, the C=O stretch near 1725 cm⁻¹ (Volkov et al., 2021), and 1 : 1 and 2 : 1 clay minerals, which provide surfaces for organic matter adsorption . Absorptions for quartz and other minerals in the fingerprint region were also important in the models, but negatively affected the estimates. The C_def model drew on information on C already present in the soil, which contributed negatively and on soil mineralogy, which indicates what soil minerals could potentially adsorb, contributing positively.

The spectroscopic MAOC and C_def models were developed using CUBIST, which offers good predictability and interpretability, effectively handling non-linearities, and is advantageous compared to linear methods like PLSR. Specifically, CUBIST performs data-driven selection of relevant informative spectral features from the full spectrum input, as well as utilising contextual information across spectral regions from regions without distinct absorption peaks. CUBIST is therefore advantageous, given that mid-IR spectra are high-dimensional and contain regions that vary in information content, and peak positions can shift under varying molecular environments. As a tree-based algorithm, it can be locally interpreted, unlike other algorithms that are limited to global-level interpretation (Viscarra Rossel and Webster, 2012). SHAP values provided additional interpretation, allowing us to not only know how each wavelength contributes to the model and how strongly they contributed to it but also show what direction an increase or decrease in absorbance affects the model, thus identifying which soil constituents (clay minerals, quartz, and organic C) significantly contribute to determining MAOC and C_def. Nevertheless, given the heterogeneity of soil composition, overlapping absorptions make it challenging to distinguish molecular vibrations, particularly in the fingerprint region. Like other regression tree methods, CUBIST can be sensitive to strong collinearity, potentially leading to model instability and overfitting (Kuhn and Johnson, 2013). To minimise the effect of collinearity in our modelling, we interpolated the spectra to a resolution of 32 cm⁻¹ (see Methods section).

This study extends beyond previous research by incorporating samples from various other ecosystems. The samples span Australia's main Köppen-Geiger climate zones, 24 major vegetation groups, and 11 of the 14 Australian soil classification orders (Isbell and the National Committee on Soil and Terrain, 2016). We excluded hydrosols with different C-stabilisation dynamics. Future work will include more samples and a broader representation of soils to develop site-specific C_def estimates. Although MAOC was measured in the ≤ 50 µm fraction and clay + silt content follows the Australian classification (clay <2 µm, silt 2–20 µm), the practical impact of this mismatch is likely modest given that Australian soils tend to have low silt contents and fine fractions dominated by clay-sized particles. Nevertheless, future work should align these operational definitions where possible, for example by directly measuring and modelling ≤50 µm or ≤53 µm silt + clay fractions, to better reflect the true mineral capacity for C stabilisation.

Our method facilitates efficient data acquisition, providing an effective approach to help farmers and land managers gain the insights needed to assess the current and potential for C sequestration on their land. Identifying regions and soil types where increasing organic C storage is feasible enables more targeted resource allocation and informed decision-making.

While our study pertains to Australian soils, the principles of applying laboratory-based mid-IR spectroscopy and machine learning to estimate MAOC and C_def are applicable across various land uses, soil types, and climatic conditions. This approach provides high-throughput MAOC and C_def estimation on sampled soils. Furthermore, such laboratory models can, in the future, underpin and improve the calibration and validation of remote-sensing-based approaches. When used in combinations, these methods provide the rapid assessment capability needed to scale soil C initiatives for monitoring soil organic C and its potential contribution to climate adaptation and mitigation targets under the Paris Agreement and the UN Sustainable Development Goals. The method's ability to support large-scale monitoring of C sequestration potential also makes it relevant to soil C credit systems such as the Australian Carbon Credit Units (ACCU) scheme.

5 Conclusions

We demonstrated that mid-IR spectroscopy combined with machine learning could effectively estimate soil MAOC content (RMSE = 2.77 g kg soil⁻¹, R² = 0.86, CCC = 0.91) and C_def (RMSE = 3.72 g kg soil⁻¹, R² = 0.89, CCC = 0.94). We interpreted CUBIST, confirming the contributions to the models from functional groups related to organic functional groups, clay minerals, and quartz, reflecting existing soil organic C, soil mineralogy, particle size distribution, and surface area available for C adsorption, which are critical for estimating MAOC and C_def. Our approach contributes to the analysis of C sequestration potential using mid-IR spectroscopy and machine learning, supporting the development of rapid and cost-effective soil C sequestration assessment and monitoring.