the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 01 Jun 2026
Compositional spatial modelling of soil organic and inorganic carbon fractions with calibrated joint uncertainty propagation
Abstract. Farm-scale soil-carbon assessments require more than a map on total carbon. They need the organic fractions, the inorganic pool, and calibrated uncertainty around each estimate. We developed a probabilistic compositional framework that propagates uncertainty jointly from mid-infrared (mid-IR) spectroscopic predictions through probabilistic trend and Bayesian spatial modelling. The framework preserves closure among particulate organic carbon (POC), mineral-associated organic carbon (MAOC) and an instrument-defined residual organic carbon (ROC), and preserves mass balance among total organic carbon (TOC), total inorganic carbon (TIC) and total carbon (TC). We applied the framework at a Mediterranean-type semi-arid farm to map POC, MAOC, ROC, TOC, TIC and TC at 0–10, 10–30 and 0–30 cm. Spectroscopic uncertainty was represented by bootstrap prediction distributions, propagated through Natural Gradient Boosting (NGBoost) trend models and Bayesian spatial models based on stochastic partial differential equations (SPDE), estimated using the Integrated Nested Laplace approximation (INLA). Predictive calibration was strong: 95 % probability-integral-transform (PIT) coverage was 0.94–0.95 across all response-depth combinations. Posterior intervals also bracketed bulk laboratory measurements (Kling–Gupta efficiency, KGE 0.64–0.79) and independent measurements (KGE 0.12 for ROC to 0.66 for MAOC, and up to 0.84 in the managed-pasture cohort). The maps showed consistent land-use effects on organic carbon. Cropping, managed pasture and natural vegetation formed the ordering crop < managed < natural for every organic-C pool and depth, with the largest deficits at the surface. Cropping also shifted composition toward the protected pools, with a lower labile-to-protected ratio (POC/[MAOC+ROC]) than pasture. TIC and ROC showed little land-use contrast. Spatial controls differed among pools: gamma-radiometric ratios dominated MAOC, electromagnetic induction conductivity dominated POC at depth, and topographic redistribution organised pools integrating multiple mechanisms. The calibrated posterior, rather than the point estimate, is the appropriate basis for soil-C management, monitoring and accounting.
Competing interests: R.A.V.R is an executive editor or SOIL
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.- Preprint
(45394 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
- RC1: 'Comment on egusphere-2026-2905', Anonymous Referee #1, 19 Aug 2026
-
RC2: 'Comment on egusphere-2026-2905', Anonymous Referee #2, 28 Sep 2026
This is a well-written manuscript that addresses a relevant problem, using advanced methods and a solid case study. The methodology is highly complex. Perhaps too complex. I lacked the necessary background to do a thorough review of the methodology. I strongly advise that an expert in Bayesian hierarchical modelling reviews the paper as well.
I did not review the entire paper but do have many comments. These are mostly minor to moderate but there are also some that in my view require a major revision.
L47 I don’t think you can claim that methods are now capable to “produce predictions accurate enough to inform management at sub-hectare resolution”. It all depends on the case at hand (number of observations, predictive power of covariates, degree of residual spatial autocorrelation). We often still find that prediction uncertainty is large.
L49-51 I do not understand this sentence, how can measurement error be a single point of truth and how can this be used to estimate the trend? I guess the word “error” should not have been there? I also do not understand that you write that a predictive distribution is collapsed into a single residual variance. If uncertainty is ignored there is no predictive distribution and no variance, we only have a point prediction. Note also that DSM often reports maps of quantiles of the predictive distribution, so more thatn just a variance. And is this paper not also guilty of reducing a predictive distribution to a variance (see L250-251)?
L51-52 This sentence also needs better explanation to be understandable. What do you mean by “aggregated estimates”? Please be also more consistent throughout the paper in the use of “estimate” and “prediction”, these are not synonyms.
L59 I guess this is only one order of magnitude given that the cost reduction is factor 10?
L100 Replace “by” by “in” (you refer to a publication, not to people)?
L105 Measurements were composite samples from three cores. This is the only time in the entire paper that mention is made of the spatial “support” of observations and predictions. Given the strong focus on prediction uncertainty and uncertainty propagation, I would have expected that more attention was paid to this. It is well-known that uncertainty is highly dependent on the spatial support of the observations and predictions, with much lower uncertainty at larger spatial supports (e.g. compare point and block kriging variances). Perhaps the case study should not address uncertainty at ‘point’ support but at ‘block’ support? I am not asking that such analysis is done, but the Discussion could make clear that results strongly depend on the support and that users are likely more interested at predictions and prediction uncertainties at block support than point support.
L125 & L150 Be careful with the term ‘representative’. Better specify what you mean by this.
L150 Is there any overlap between this subset of n=200 samples and the subset of n=100 samples mentioned in L125? And are these n=200 samples evenly distributed over topsoil and subsoil?
L166-174 It was not clear to me that data leakage was prevented. Usually one uses a nested CV-approach for that (splitting the dataset in three: calibration, validation and test), but here it seems the dataset was only split in two?
L175-181 Spectroscopic uncertainty was quantified using bootstrapping. The text is very condensed and it was not clear to me how it all worked in detail, but I fear that the procedure used severely underestimates the spectroscopic prediction uncertainty. Bootstrapping means creating datasets that have the same size as the original dataset, using sampling with replacement from the original dataset. If uncertainty is quantified from the variability between bootstrap results, then this uncertainty is only representing sampling uncertainty (i.e., the fact that only a finite dataset is available for model training and prediction). It will not capture other, more important sources of uncertainty, for example the limited capability of spectra to predict a soil property, even with an infinitely large training dataset. Note that the variance estimated by bootstrap sampling converges to zero as the sample size tends to infinity. In summary, a major concern I have is that the bootstrap method only quantifies the sampling uncertainty, which is only a (small) part of the overall prediction uncertainty. I would expect that this comes out from calculating the PICP. Authors use PIT, which I guess is quite similar, but they only do so for the final results (see L371-374). I ask that this is also done for just the spectroscopic uncertainty.
L175-181 One other key problem is that the prediction errors associated with spectroscopic predictions are likely correlated (i.e., not statistically independent) between samples. It seems this study ignores this completely and assumes that they are all independent. The consequence of this is that spectroscopic uncertainty hardly propagates to DSM, because uncertainties largely cancel out. But is this realistic?
L178 Why call this a ”discrete predictive distribution”? The variable of interest is continuous-numerical, so we need a continuous probability distribution. I think a better formulation is to write that you consider the 30 predictions a random sample from the continuous conditional distribution of the variable of interest. And is not B=30 a much too small number?
L196-198 I realise that reduction of number of covariates using PCA is commonly done but I am sceptic about the validity of that. Variables that are not correlated with others will not contribute much to the first PCs and so will not be retained if you only keep PCs that have at least 10% of the variance. But such variables might well have strong explanatory power, so why delete them without checking? Many machine learning methods have in-built protections against overfitting, so I would delete this step.
L204-211 Again considerable tuning, how did you prevent data/information leakage? This requires that test data were removed before tuning, but I did not see that mentioned yet.
Section 2.5 I am not an expert in compositional data so could not verify the soundness of this section. I thought the problem was easier and that log-ratio transformation already does the trick. For instance, when applied to soil texture (sand, silt, clay) we can use the ALR transform to reduce the three fractions to two variables, do the modelling and prediction on those two, and later back-transform to the original scale (such as done in Poggio et al. 2021)?
L232 Is it realistic to assume that log(TOC) and log(TIC) are independent? And even if this assumption is realistic, I fear that this does not automatically imply that the prediction errors of these two variables will also be independent. Please check.
L236 It was not clear to me what “multiple-imputation” means. It is only explained much later, in Section 2.7, but perhaps this could be clarified earlier on. As I understand it, it is no more than a Monte Carlo loop over the B=30 bootstrap samples, is this correct?
L241-242 Here you assume that spectroscopic prediction errors are statistically independent, which is likely unrealistic. See also an earlier comment. At least this assumption should be made explicit. If I understood correctly, spectroscopic uncertainty is not accounted for in trend modelling with NGBoost, but only in the kriging part through Eq. 3. Perhaps it should indeed not be included in trend modelling (part 2) because you take an approach where you sample from the conditional distribution of the ‘true’ response variable, given spectroscopic uncertainty (i.e., stochastic simulation), but if that is true we also should not account for spectroscopic uncertainty in the spatial modelling (part 3) (i.e., there should be no variance inflation of the epsilon in Eq. 2)? I also did not understand the rationale behind Eq. 3. Is this equation theoretically derived or a heuristic approach?
L241 It was not clear to me why INLA-SPDE was used instead of kriging. What is the added value? Maybe this is needed to take a Bayesian approach, but where does the Bayesian come in and is this justified (is the added complexity worth the effort)? Which parameters of the spatial model were treated as stochastic and how large were the variances of their posteriors? If small it is not worth the effort. I also noticed that not much effort was placed in defining priors. Please note I am not an expert in Bayesian hierarchical modelling and could not check the validity of these parts (i.e. Section 2.6.3).
L248-252 It was not clear to me whether you used a stochastic or deterministic trend model. L249 mentions that each fit returns a predictive distribution but is this because the input to NGBoost is stochastic (i.e., the 30 bootstrap predictions from the spectroscopic model) or is this because NGBoost accounts for its own prediction uncertainty? If the latter, how is this done? And why is only the standard deviation of the full predictive distribution used? And why is NGBoost’s per-point predictive variance only used to down-weigh the fit of the Bayesian spatial model through an increase of the residual variance in Eq. 3 (L282-284)? Can this approach be mathematically derived?
L261-262 I think it is bad science to remove outliers only because they are influential and distort the trend. In my view it is only justified to remove outliers if there is sufficient evidence that something went wrong during measurement or processing (e.g. negative SOC) or that the observation location does not belong to the population of interest (e.g. the location is on a road or manure heap).
L263-265 I do not understand why this prevents overfitting, but this might be because I am not an expert in these methods. Explanation is needed. Maybe it is explained in Section 2.7?
L217 Much of Section 2.6.3 is beyond my knowledge and so I could not check it. I am still not convinced why INLA-SPDE was needed and could not be replaced by the much easier residual kriging. The added value must come from making parameters of the random field stochastic, but is not parameter uncertainty negligibly small in this study? Could you compare your results with those obtained using a simple residual kriging approach?
L278 Why not force beta_0 to 0 and beta_1 to 1? I hope you did allow w(s) to include a nugget (short-distance variability)? Authors will know that short-distance spatial variation is often large for soil properties, so nugget is more than random measurement error. Regrettably, the equation in L296 suggests that this was not included.
L331 Why replace B by M?
L331-332 Here I guess independence between specimens was assumed, which might not be realistic.
L342-347 Text was not entirely clear to me but am I right that you used a stochastic simulation (instead of prediction) approach to derive quantiles of the back-transformed variables? It seems you did because L334 mentions that you took 500 draws from the posterior distribution. I think this makes sense, and it would be interesting to separate the uncertainty resulting from spectroscopic uncertainty from that of other sources, and I think you do so in Section 2.8. I can see that separating out the spectroscopic uncertainty from the M x 500 draws is easy, but I do not see how you can separate out the other four sources uncertainty (as mentioned in L352) from these draws. Is this clearly explained somewhere in the paper?
L364-365 Not clear what CV you refer to here and that data leakage is prevented. As mentioned above, I would assume a nested cross-validation where outer held-out folds intended for performance evaluation are created right at the start, not used in any way during modelling, also not for model selection and hyperparameter tuning.
L389 Were these 59 sites a subset of the sites previously sampled or were these a new set of sites? I guess that for this to be a meaningful check these should be newly sampled locations, ideally a probability sample from the study area. L392 suggests this was not the case.
L424-706 I did not review this part of the paper (I had already spent much time on the review and made many comments).
Citation: https://doi.org/10.5194/egusphere-2026-2905-RC2
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 194 | 75 | 34 | 303 | 28 | 21 |
- HTML: 194
- PDF: 75
- XML: 34
- Total: 303
- BibTeX: 28
- EndNote: 21
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
The paper is well-written and the analyses are rigorous. However, the method section is very long and difficult to follow for the non-specialist. I would recommend transferring parts to an annex.
Is there any reason why the equations are not numbered anymore after equation 3?
Line 392 I there any reason why you were not able to re-sample 4 out of 59 locations within an acceptable distance?
Figs 3, 4, 10, 11 What does the red dot refer to? Please explain in the caption.
Table 2 and all tables and figures. Abbreviations should be defined in the caption
Section 4.1 I miss a consideration on the practical application. To what extend is this rather complex analysis realistic for farm management? Is there any key indicator that the authors can recommend for practical purposes?
Lines 674 and 692 What does ‘bulk soliTOC’ ?