the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 16 Feb 2026
Quantifying soil organic carbon stocks above the treeline in the Swiss Alps
Abstract. More than 90 % of the total carbon (C) in alpine ecosystems is stored belowground, yet spatial estimates of soil organic carbon (SOC) stocks remain scarce due to limited accessibility and the demanding nature of SOC stock estimates in rocky alpine terrain. By combining new measurements at 144 sites across the Swiss Alps with data from existing inventories, we compiled a comprehensive dataset on SOC stocks totalling 307 sites from treeline to the permafrost region (1750 m – 3100 m a.s.l.). We predicted the spatial distribution of SOC by linking stock measurements to environmental covariates using Quantile Regression Forests (QRF) and produced a SOC stock map at 25 m resolution illustrating the spatial SOC variability in alpine terrain. Our results show that SOC stocks average 7.3 ± 3.3 kg m⁻² in alpine grasslands and 1.8 ± 1.7 kg m⁻² in partly vegetated areas around and above the vegetation line. Overall, the alpine region of Switzerland, which covers one-third of the total country area, stores an estimated amount of 47.6 Mt SOC, representing a non-negligible share of the Swiss greenhouse gas inventory. Vegetation productivity, represented by the Normalized Difference Vegetation Index (NDVI) and topo-climatic covariables, together with vegetation-derived indicators of humus content and soil pH, were highly informative for spatial predictions. This study identifies hotspot regions of SOC storage and influential spatial predictors of its distribution, providing a quantitative baseline for assessing the status-quo and future changes in alpine SOC stocks under continued climate and land-use change. The observed increase in SOC stocks with increasing NDVI suggests that climate change-driven greening at high elevations, where vegetation cover is currently sparse, may enhance SOC storage, although the rates and magnitude of these changes require further investigation.
Competing interests: One of the coauthors, Frank Hagedorn, is a member of the editorial board of Biogeosciences.
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
(4276 KB) - Metadata XML
-
Supplement
(106999 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2026-835', Anonymous Referee #1, 12 Mar 2026
-
AC1: 'Reply on RC1', Michael Zehnder, 29 May 2026
Dear anonymous reviewer,
We sincerely appreciate the constructive feedback on our manuscript. Should we be granted the opportunity to submit a revised version of our work, we are confident in our ability to address all the concerns raised and to substantially improve the quality of the submitted manuscript. Please find our responses to each comment and suggested changes. Comments of the referee are included in italic font.
Kind regards on behalf of all listed co-authors,
Michael Zehnder
This study quantifies soil organic carbon stocks above the treeline in the Swiss Alps by combining 307 field sites with pedotransfer-based subsoil extrapolation and spatial modelling using Quantile Regression Forests. It estimates 47.6 Mt SOC stored across alpine terrain and identifies NDVI, vegetation-derived humus indicators, temperature, and topography as key predictors of SOC spatial variability.Overall, the manuscript presents a valuable and timely contribution to SOC analysis in alpine ecosystems. It is well-written, provides an important baseline for SOC stocks and fits very well the scope of the journal. However, clarifications are needed before the manuscript could be considered for publication, especially regarding the robustness of PTF for subsoil extrapolation, the predictor circularity, the coarse fragment integration, the pseudo-zero weighting, and the full uncertainty propagation.
We thank Referee 1 for the positive, constructive, and thoughtful comments on our manuscript. Below, we address the points raised and hope that our responses help clarify the remaining uncertainties. Most importantly, we clarify the concerns of the reviewer regarding PTF by moving detailed description from the Appendix to the main text, and by providing a SOC map from 0-20 cm depth, which does not rely on PTFs and provides additional information on surface soils.
Main concerns
1- First, a major methodological concern relates to the use of pedotransfer functions (PTF) to estimate subsoil SOC stocks. Subsoil values are extrapolated from only 17 soil profiles stratified by bedrock type and elevation bands, yet these relationships are applied to 307 sampling sites across the alpine region. Given the well-known heterogeneity of alpine soils, including high skeletal content, cryoturbation, discontinuous soil depth, and highly variable pedogenic pathways, the robustness and transferability of these PTF require more rigorous justification and uncertainty assessment. Another option is to move all this analysis into supplementary information. In particular, the procedure of selecting the “best-fitting” profile for each site may introduce overfitting and underestimate the structural uncertainty associated with subsoil extrapolation. Related to this, the manuscript does not clearly explain how soil depth is estimated or constrained across the study area. Without explicit constraints on soil depth, extrapolation of SOC stocks below the measured 0–20 cm layer becomes difficult to evaluate. Moreover, soil types such as Umbrisols, Podzols, Rankers, Rendzinas, and Cambisols can occur over short spatial distances and exhibit very different vertical SOC distributions and bulk density patterns.
We appreciate the points raised and agree that the referee has identified important aspects of our decision to use PTFs to approximate SOC stocks below the measured depth of 20 cm. To address these concerns, we will move relevant information from Appendix B to the main text for greater transparency.
To clarify our PTF procedure: we did not simply select the best-fitting function and applied it to a given site, which would limit the fitted depth-trajectories of FE density, carbon content and the coarse rock fragment fraction to the exact values measured at the 17 profile sites. Instead, after identifying the environmentally matching profile sites for a given topsoil site, the depth-function from environmentally matching profiles (in most cases 3-5), is fitted to the measured values of the topsoil sites at different depths (see Figure 2, but also Appendix B; Figure B5). We will clarify our PTF procedure in the revision of the manuscript as explained further below.
Before responding to the individual points in detail, however, we would like to clarify what appear to be some misunderstandings regarding the PTFs, as reflected most clearly in two comments raised further below. We address these two points first to provide context for our subsequent responses.
-“If I have well understood your study, you first spatialized SOC stocks in the upper 0-20cm horizon using QRF model, then you used PTF to assess SOC stocks for the whole profile, right?”
No, this is not right, and this is a critical point for understanding how we used PTF in our study. Very importantly, the PTF were implemented before the spatial modelling on site-level, and not on pixel-level for specific elevational or environmental ranges as apparently assumed by referee 1. To avoid this misunderstanding, we will address the PTF step more explicitly in the Method section and add it to the workflow illustration (Figure 1) for improved visual clarity, as also kindly suggested by Referee 1 in a later comment.
- "L177: Why did you choose linear regression for FE density while you used exponential decline for SOC content?"
The type of regression was chosen based on the type of regression that best fit the measured data in the profile sites. In the case of FE density, linear regression fit the best. For SOC content, exponential functions fitted better. In paragraph “2.5 Pedotransfer functions (PTF)” (line 170) we repeatedly referred to Appendix B), which can be found at the end of the manuscript, for detailed descriptions of the PTF approach, where we described the linear relationships for FE density and exponential relationships with SOC content (Figure B2, B3 and B4). As the figures B2 - B5 shown in Appendix B) required a lot of space and are not critical for the main storyline, which is the spatial modelling of SOC stocks rather than PTF, we moved the more detailed description of PTF to the appendix and cross-referenced this in the main text.
Suggested changes related to PTF:
Given the apparent unclarity regarding the implementation of PTF, we think that it will significantly benefit the manuscript if we will describe the PTF approach better and more detailed in the main text, and not just in Appendix B. First of all, we suggest increasing the general level of detail in the description of the procedure in the main text by elaborating on the fitting procedure. At the same time, we believe that most of the unclarities occurred because referee 1 understood that PTF were applied on a pixel-basis – a point that we will improve in the manuscript by explicitly including it at the right place in the overview figure, and explaining the order of steps (i.e. first step: PTF implementation on site-level based on measured topsoil stocks, second step: training the spatial model with PTF-modelled subsoil stocks below 20cm added to measured SOC stocks in the uppermost 20 cm, third step: spatial extrapolation, whereby the spatial model already learnt the relationships between all the variables and total SOC stocks down to bedrock) better in the paragraph on PTF functions.
In direct relation to unclarities with PTF implementation, we will elaborate on how soil depth is estimated in the method section 2.5. We will add the following few sentences to the PTF approach description in the main text:
“If soil depth exceeded the sampled depth of 20 cm, soil depth was estimated through the linear increase in the PTF-modelled coarse-rock fraction. Once this fraction reaches 100 %, the parent material is reached. The exact maximum soil depth thus varies from site to site and is dependent on the measured values of coarse rock fractions and how these change within the 0-5, 5-10 and 10-20 cm intervals, but also on which of the environmentally matching subset of the total 17 profile sites best fitted the measured values used to model the increase of the coarse rock fraction.”
As mentioned at the beginning of our response, we will provide a map of the measured 20 cm topsoil, which is independent of the PTF approach.Lastly, we would like to emphasize that soil sampling below 20 cm in alpine regions is a very strenuous task (particularly the excavation of soil profiles), which is possibly one of the reasons for scarcity of soil data at the profile level in alpine regions. The number of high elevation soil profiles (going up to 3000 m a.s.l.) used in our study to constrain PTF and the overall number of more than 300 samples is unprecedented and forms a strong basis for the spatial product provided.
2- Second, Landolt-H and Landolt-R are partly circular predictors: both are vegetation-derived proxies likely co-varying with SOC. Their strong importance may reflect shared response rather than mechanistic control, potentially inflating predictive performance while limiting process inference. Moreover, the uncertainties of these variables are not integrated in the modelling approach, and I assume that these uncertainties are spatially distributed.
We thank the referee for raising this important point. We fully agree with this concern and intended to make clear throughout the manuscript that the relationship involving the vegetation-derived humus indicator reflects co-variation rather than a mechanistic link. For example, in line 519 we explicitly stated that “the high rank in the importance metrics reflects shared variation with SOC rather than an independent mechanistic driver or a controlling factor of SOC.” See also our phrasing in line 515: “Landolt-H was frequently used to partion variation”.
We also emphasized that vegetation-derived spatial information (available as point data at much broader spatial coverage, i.e. millions of georeferenced humus-indicator values, and translated to a robust country-wide map by Descombes et al. 2020; line 213 and 520) is used as a proxy to support the spatial model, rather than as a mechanistic predictor of SOC dynamics. We will further clarify the distinction and the use of Landolt values for mapping purposes in both the discussion and the description of the spatial covariates used by adding the following to statements after line 519:
- “Landolt-H and is used here purely as spatial mapping tool, i.e. to improve the spatial prediction of SOC stocks, and are not intended as mechanistic predictors of SOC.”
- “Variable importance metrics reported here reflect each covariate's contribution to model predictive performance and should not be interpreted as evidence of mechanistic processes.”
Concerning Landolt-R, our pH proxy; This variable arises from the fact that there are simply no existing pH maps or soil mineralogical maps of the surface (quaternary) geology available for our study area above the treeline (see, line 212). The use of plant-derived pH proxies is therefore the best available option and significantly improved model performance (Appendix A, Figure A4). This also means that the data required to include a mechanistic discussion, as suggested by the referee in a comment further down, is not available on a spatial basis, as exemplified by the use of Landolt-R. These limits (or even prohibits) the use of terms “drivers” or mechanistic link as we fully agree. We have discussed this limitation explicitly in line 505 and consistently refer to Landolt-R as a pH proxy, while avoiding terminologies such as “predictive power” or “mechanistic link” for these variables, where no such mechanistic link exists. We are happy to make this limitation even more explicit in the manuscript. On the other hand, we think that our innovative approach of using spatially available vegetation-derived data has many benefits for spatial mapping of soil properties, which we have only partially exploited in this manuscript.
We are not entirely sure what referee 1 is referring to regarding the uncertainties of Landolt-H and Landolt-R not being integrated. Our modelling framework does not differentiate between uncertainties associated with plant-derived variables and those of other spatial predictors. Instead, the uncertainty analysis (based on the standard deviation derived from quantiles estimated by the quantile regression forest) treats all predictor variables consistently, which is a standard practice in this type of modelling approach. If uncertainty distributions were available for each individual input layer (e.g., CHELSA temperature, CHELSA precipitation, NDVI, etc.), it would in principle be possible to propagate these uncertainties through computationally intensive Monte Carlo simulations. However, this is rarely implemented in practice and is not standard in digital soil mapping workflows. The focus of our paper is to quantify SOC stocks above the treeline using existing spatial tools, not to develop or refine geostatistical methodology and spatial uncertainty propagation (see also our response to major comment nr. 4). In the vast majority of studies using quantile regression forests, uncertainty in spatial input layers is implicitly assumed to be comparable across predictors, and overall uncertainty is instead approximated using the model-internal quantile distributions, which is the approach we applied here. We will, however, add a sentence acknowledging that the provided uncertainty map reflects model prediction uncertainty only and does not account for uncertainties inherent to the input covariate layers.
Lastly, we implemented a spatial blocking in the cross-validation of the QRF model to reduce possible spatial autocorrelation. We will add a remark on this in line 251: “The data were partitioned into five spatially contiguous folds based on coordinate proximity, such that, in each iteration, 80 % of the data were used for model training and 20 % were held out for spatially independent validation using the blockCV package (Valavi et al., 2019) (Figure A2), reducing the risk of performance overestimation due to spatial autocorrelation between training and validation samples."
3- Third, the treatment of coarse fragments and bare-rock areas raises methodological concerns that may bias SOC stock estimates. The sampling protocol did not allow the quantification of coarse fragments larger than 5 cm, which are common in alpine soils. In addition, the treatment of bare-rock areas through the introduction of “pseudo-zero” SOC values based on expert judgement lacks sufficient methodological justification. This assumption directly influences SOC predictions in sparsely vegetated, high-elevation zones and may significantly affect regional SOC stock estimates. It is unclear why bare-rock areas were not explicitly excluded from SOC calculations, unless the study explicitly considers petrogenic organic carbon or organic matter stored in lithic substrates.
We thank the reviewer for pointing out this lack of clarity. While this may not have been described explicitly enough, our sampling protocol did indeed allow us to account for rock fragments larger than the corer diameter. In fact, this is one of the distinctive features of our approach at IMIS sites, and we agree that it should be emphasized more clearly in the manuscript. We therefore suggest adding the following explanation to the methods section:
“Four spatial soil sampling replicates were taken, whereby each replicate consisted of three pooled cores, collected within approximately 1m distance in every quarter, at depth intervals of 0–5, 5–10, and 10–20 cm (Figure C3). Notably, if one of the individual soil cores hit rock, the unextracted proportion and the volume below this depth was treated as rock, therefore allowing the quantification of the volume of coarse rock fragments larger than the corer diameter of 2 cm.”
We hope this clarifies the issue and highlights the value of this approach in alpine terrain, where coarse rock fragments indeed play a major role in determining SOC stocks.
Regarding the pseudo-zeros, we believe there may be some misunderstanding. This is also reflected in a later comment where Referee 1 asks whether pseudo-zeros contain the value “0”: Yes, they do. We suspect the current manuscript may not sufficiently explain how pseudo-zeros influence predictions in bare rock areas or, more generally, at higher elevations, and what we define as “bare rock” in this context. We believe these points can be clarified by revising a few sentences in the manuscript (see below), and we thank the referee for highlighting these ambiguities.
First of all, the classification of "bare rock" is defined by solid rock cover of >= 80% (Appendix D; Table D1) and is spatially defined within the TLM3D layer provided by SwissTopo – the standard reference for Switzerland’s land cover classification. This definition indicates that even within the cover category of "bare rock", there can be pockets with substantially vegetated cover, with non-negligible amounts of SOC, especially in the lower alpine belt! We will add a more exact definition of bare rock in line 295. (“Within the TLM3D land cover classification, "bare rock" is defined as areas with solid rock cover of at least 80%, meaning that up to 20% of such areas may still support vegetation and non-negligible SOC stocks (Table D1)”). Moreover, we have noticed that Figure 3, containing the representative images of the cover categories, may falsely imply that "bare rock" is exclusively found at very high elevation and dominated by grey, solid-rock landscapes, as depicted in Figure 3c. We will add a better picture that is less misleading.
Furthermore, the reason for the relatively high total SOC stocks in this land-cover class are relatively high (and this correctly so) is that many areas classified as “bare rock” occur in the lower alpine belt, where rocky terrain often contains vegetation pockets, heathlands, or other patches with substantial SOC stocks. The comparatively high model predictions in these “bare rock” areas are therefore not erroneous, as the model is trained on continuous predictors such as NDVI, elevation, and related environmental variables and not the categorical land-cover classes. The land cover classes mainly serve to calculate the total SOC stocks for a given land cover category, for instance to delineate the treeline, and to crop out ice-covered areas (line 303). Strictly removing all pixels classified as “bare rock” would therefore exclude a substantial proportion of SOC-rich areas at lower alpine elevations from the final map. While excluding bare rock in nival regions could arguably be justified, we did not identify a straightforward and objective way to distinguish these areas from rocky terrain in the lower alpine belt. We therefore chose to keep the land-cover classifications provided by SwissTopo. That said, we agree that this point should be explained more clearly in the discussion and also in “2.9 Spatial mask and land cover classes”, particularly why SOC stocks in the “bare rock” category appear relatively high. At the same time, we recommend keeping these areas included in the final map and not masking them out. Users inspecting the underlying `.tif` output will see that predictions for bare rock areas above approximately 3200 m are close to zero (or zero), whereas rocky areas around 2300 m often show substantially higher values, sometimes comparable to grasslands, which we consider to be well within realistic conditions given the definition of “bare rock”. In the revised manuscript, we will add the following sentences to explain the relatively high SOC stocks in the “bare rock category” discussed in line 439:
“The comparatively high total SOC stocks estimated for the bare rock category reflect the fact that a substantial proportion of areas classified as bare rock occurs in the lower alpine belt, where vegetation pockets and organic-rich soils within the permissible ≤20% non-rock cover (Table D1) contribute meaningfully to SOC stocks in this class.”
The referee raised the concern that the “treatment” of bare rock areas using pseudo-zeros may introduce bias and lacks sufficient justification. We would like to clarify that pseudo-zeros were not included in our initial modelling workflow (to answer the question of the reviewer further below). However, models trained without them produced notable overestimations of SOC stocks at high elevations because our field sampling did not include truly bare, exposed rock surfaces where no soil was present – for practical and logical reasons, as discussed in line 429. As a result, the quantile regression forest model learned that some SOC is always present, which is clearly unrealistic for large parts of the nival zone dominated by exposed rock, where organic matter stocks are minimal or absent (excluding petrogenic organic carbon, which was beyond the scope of this study). Accurately representing these nival bare rock environments would have required direct sampling of such locations. Since this would ultimately serve the same purpose as manually introducing pseudo-zero observations during model development, we opted for the latter approach. The inclusion of pseudo-zeros substantially improved model performance and better calibrated predictions in nival regions, where the model now frequently predicts values close to zero or exactly zero. Again; without pseudo-zeros, predictions were constrained by the lowest observed SOC values in the training dataset and therefore systematically overestimated bare rock areas. We therefore agree with referee 1 that pseudo-zeros meaningfully influence regional SOC stock estimates, but we emphasize that they do so in a direction that reduces, rather than introduces, bias from the absence of truly exposed rock samples in the field data.
We agree that the exact number of eleven pseudo-zero sites included appears arbitrary. The number of eleven was simply chosen to match the number of sites categorized (by the landcover class layer of SwissTopo) as on “bare rock” (n=11, line 300), with the goal to train the model with truly non-SOC containing pixels in an approximately equal number, while giving the pseudo-zeros three times higher weights. A detailed assessment of the optimal numbers of pseudo-zeros required for model optimization, and the exact weighting of these pseudo-zeros in the training of the model exceeds the scope of the manuscript. To clarify that we did not sample exposed bare rock with no soil in the field in line 300, we will add the following sentences:
“While permanently snow- and ice-covered areas and completely exposed bedrock contain no soil and were therefore not sampled in the field, we added a set of expert-led pseudo-zero points (n = 11) located above 3150 m a.s.l. in permanently snow- and ice-covered areas and on bare rock to inform model behaviour in these high-elevation environments as suggested by Wadoux et al., (2020). Pseudo-zeros can be used to constrain model behavior in areas where SOC stocks can be assumed to be zero with high confidence, as applied for example in SoilGrids 250m (Hengl et al., 2017), where hundreds of pseudo-zero observations were used in desert and ice-covered regions.”
Hengl et al (2017): https://doi.org/10.1371/journal.pone.0169748
4- Finally, although the geostatistical framework appears robust and overall well implemented, the treatment of uncertainty remains incomplete. The reported total SOC stock estimate does not appear to fully propagate the different sources of uncertainty inherent to the methodology. In particular, uncertainties associated with pedotransfer functions, covariate uncertainties, and model structural uncertainty are not explicitly integrated into the final confidence intervals. In addition, the assessment of spatial autocorrelation in the model outputs appears to rely primarily on visual inspection and not on formal statistical evaluation (variogram analysis or Moran’s I statistics on residuals).
We thank the referee for the positive comment on the general implementation of our spatial model and we confirm that the analysis of residuals was done visually (but see also Appendix A, Figure A5). We will implement a Moran’s I statistics analysis on residuals.
We agree that each covariate in our model carries some sort of uncertainty. A comprehensive, joint propagation of uncertainties in heterogeneous covariates has so far hardly been established in practical DSM/SOC studies and is rarely implemented; existing approaches are mostly case-study- or component-specific. Even topographic data derived from high-accuracy digital elevation models, carry an inherent uncertainty, yet these uncertainties are not propagated into final mapping products in the vast majority of DSM studies to the best of our knowledge, e.g. none of the practical studies cited in our manuscript included uncertainty propagation on the level of individual covariables. In fact, this is an ongoing challenge of DSM, which is happening, however, rather on a theoretical level:
E.g. see at the end of the editorial: Mulder, V. L., Roudier, P., & Arrouays, D. (2023). Editorial: Digital soil mapping – advancing the knowledge frontiers. Frontiers in Soil Science, 3, 1225672. https://doi.org/10.3389/fsoil.2023.1225672
Or point 5 in “4 Summary and future work”: Xia, Y., McSweeney, K., & Wander, M. M. (2022). Digital mapping of agricultural soil organic carbon using soil forming factors: A review of current efforts at the regional and national scales. Frontiers in Soil Science, 2, 890437. https://doi.org/10.3389/fsoil.2022.890437
Or: Takoutsing, B., Heuvelink, G. B. M., & Stoorvogel, J. J. (2022). Accounting for analytical and proximal soil sensing errors in digital soil mapping. European Journal of Soil Science, 73(2), e13226. https://doi.org/10.1111/ejss.13226
Moreover, for the majority of data layers included in our analysis, there exists no spatially available error quantification in the corresponding source publications. As opposed to the input data used, we actually provided an error quantification on pixel-level for our spatial product by using quantile regression forests that keep the full spread of model predictions (as opposed to traditional random forests).
We agree that the use of PTF to approximate subsoil stocks carries inherent uncertainty. However, we acknowledged this in the section starting in 549 (“…the accuracy of PTF-based extrapolation depends on the similarity of local conditions to those represented in the reference dataset of 17 deep soil profiles by Udke et al (submitted). Sites to model subsoil stocks were chosen based on best-possible environmental match between each topsoil site and a corresponding soil profile site (Figure B5). While the set of profiles used for PTF cover a large environmental gradient, they do not necessarily reflect the exact in-situ conditions at all topsoil sites, thus introducing an additional model uncertainty, which we did not additionally validate through excavating additional soil pits at randomly picked topsoil sites, as this exceeded the extent and labour of this study.”). Again, we believe that a large part of unclarity related to PTF can be avoided by the more detailed description of PTFs in the revised Method section (see our answers related to PTF).
Specific comments
L11: “… of the total [organic] carbon…”, because in carbonated context, some of carbon is include in carbonates.
We agree and will implement accordingly.
L21: Using soil pH as explanatory variable is contestable, as SOC content (organic acid) decreases soil pH.
This is an important point raised by the referee. We suggest discussing this as a limitation but although organic acid decrease soil pH, geogenic parent material will still exert a dominant influence . Nonetheless, we would like to list a few studies that use spatial information of soil pH in order to map SOC distributions – some of which were co-authored by leading scientists in the field of digital soil mapping.
• https://doi.org/10.1016/j.geoderma.2024.116873. - British Columbia (BC), Canada
• https://doi.org/10.1016/j.geodrs.2021.e00437 - Switzerland
• https://doi.org/10.1371/journal.pone.0125814 - Africa, continent-wide
• https://doi.org/10.5194/bg-8-1053-2011 - France; note that the authors write: “The best candidate among soil properties [for spatial SOC prediction] would be the soil pH”Spatial approaches sometimes rely on predictors which would not be chosen in studies on mechanistic relationships.
L23-25: I agree with this statement for the upper alpine grasslands, but what about alpine grasslands near the treeline, containing the highest SOC stocks?
L23-25 specifically refers to higher elevations because the spatial patterns that we found only allow for a trend forecast for higher elevations, while the trajectory for already SOC rich lower alpine belt is more complex. Based on the data we have, we refrain from suggestive statements regarding future changes in the lower alpine belt, as this is not a key takeaway from our study. The distinction between lower alpine belt and higher nival regions is discussed in-depth in the paragraph “4.4. Outlook”.
L35: “due to the difficulty of accessing as well the high heterogeneity” → “due to the difficulty of access as well as the high heterogeneity”.
Thanks for catching this typo!
L42: Greening of alpine area is only one consequence of climate change on carbon dynamics, influencing quantity. Other consequences could be the change in land cover and thus change in SOC chemical composition, influencing also its stability.
We thank referee 1 for pointing this out. We suggest adding that land use and land cover change in general (e.g. shrubification) may additionally affect carbon dynamics in the alpine region.
L58: “suppress” à “slow, reduce,…”
We will change “suppress microbial decomposition” with “slow down microbial decomposition”.
L75-76: The end of the sentence does not follow on from the preceding sentences.
We thank referee 1 for catching this. We will try to improve the flow of this text passage.
L90: Why you mention the “time” factor above in the introduction (L53), and you didn’t include it in the analysis?
We mention time for completeness within Jenny's framework and explicitly acknowledge in Section 4.2 that no adequate spatial data layer exists for this factor. This is a limitation shared by virtually all regional-scale digital soil mapping studies.
L97: Given the elevation of study sites, could some areas be recently deforested? e.g. at 1750 m, potentially leading to legacy effect on SOC stocks not integrated by your model.
We thank the referee for raising this point. This is indeed possibly as livestock cultivation has lowered the treeline in many parts of the Swiss Alps. Accessing this historic data is, however, very difficult and we suggest discussing this as a limitation in the discussion rather than quantitatively including it in our models.
L105: The mosaic of parent material leads to wider variations in soil conditions than soil pH only, with some of them that could influence SOC stocks (Fe, Al, Ca, clay, sand contents, etc.). Your study didn’t really integrate the “parent material” factor. Is the integration of geologic properties could increase model performance?
In digital soil mapping, the variables included in any model are constrained to what data is available spatially for predictions. As previously mentioned, there is no map for soil texture or geological properties. Therefore, the Landolt-R, the pH proxy, is the best approximation of parent material available on a spatial basis. The vegetation accurately reflects bedrock material on the coarse level (dolomitic, calcarious, amphibolitic) and this variable indeed improved model performance, as would a map of parent rock classes – if it were available. However, the inclusion of a factorial map of soil properties (if it existed) comes with its own limitations for QRF, as factorial variables are notoriously difficult to integrate with in combination with continous variables in QRF.
L107: Can you give here an idea of the study area surface?
Absolutely, we can add this here.
L122: It is not clear why there is 20 sites from GLORIA. 6 summits x 4 cardinal directions = 24 sites?
We thank the referee for catching this. As one Gloria summit was located outisde the Swiss borders, it was automatically excluded as most covariable data is constrained to the country area. This reduces the number to 20 sites. We will adjust this accordingly.
L129: mean [annual] precipitations I supposed?
As described in line 196, the precipitation variable represents mean daily precipitation averaged over 1981–2010 (Zilker & Karger, 2025; https://www.doi.org/10.16904/envidat.689). We agree this is counterintuitive and have contacted the data provider while working on the analyses, who is currently developing a mean annual precipitation product at this resolution. In the meantime, this represents the best available high-resolution gridded precipitation data for Switzerland and represents the best available data – even though the unit may be not optimal.
L140: It is always problematic to assess the bulk density. Given the available data of the three datasets you have, your approach seems to be the best, although the difference in soil volume could introduce bias on the coarse fragment content estimation. Because in IMIS dataset, coarse fragments >2cm were not integrated in bulk density assessment, while in BDM only coarse fragments >4,8cm (and in GLORIA for coarse fragments >10cm) were not integrated. And given the importance of such coarse fragments in alpine soils (as we can see on the soil profile pictures in appendix B), this omission could have great importance on SOC stock calculation. Did you try to compare soil coarse fragment content between the three datasets? And the bulk density?
We fully agree that this is a limitation of soil studies, particularly those involving sampling in alpine terrain. As explained in our comment to referee 1’s 2. major comment, the IMIS sampling quite possibly is the best way of the three campaigns to get good estimated of the bulk density. There are still differences that can partly be explained by different elevational distributions of the three campaigns, about which we wrote several pages in Appendix C. However, a bias will always exist given that more than one dataset is combined (which is not uncommon in digital soil mapping, see e.g. Gupta et al. (2024) https://doi.org/10.1016/j.geodrs.2023.e00747 or Kasraei et al (2024) https://doi.org/10.1016/j.geoderma.2024.116873. Reducing the dataset to only IMIS or only BDM is (in our eyes) not a solution and we discussed these limitations accordingly. We are open to discussing these even more in-depth.
Furthermore, we suggest recommending the IMIS sampling approach for future studies explicitly in the discussion.
L145: “replica” → “replicate”?
We can adjust this.
L170: In this paragraph, it is not clear how do you predict spatially the soil depth.
As explained in our answer to the major comment, we can elaborate on this.
L171: It is not clear for me which PTF you assign for each pixel of your study area. If I have well understood your study, you first spatialized SOC stocks in the upper 0-20cm horizon using QRF model, then you used PTF to assess SOC stocks for the whole profile, right? First, I suggest adding the use of PTF in the workflow of figure 1. Second, the stratification by bedrock type and elevation is debatable. Are these two parameters the only ones that control soil properties and stratification? Coarse content is more dependent on geomorphic processes rather than bedrock and elevation.
This is addressed in our answer to major comment nr 1.
L177: Why did you choose linear regression for FE density while you used exponential decline for SOC content?
This is addressed in our answer to major comment nr 1.
L218: “humus-indictor map” → “humus-indicator map”.
Thank you for catching this typo!
L260: What does the activation of this option mean?
thank you for pointing out this unclarity. With quantreg = TRUE, the ranger QRF no longer just predicts a single mean value per pixel, but instead it retains the full distribution of individual tree predictions across the forest in the training of the model. This allows to extract any quantile of the prediction distribution at each location. In practice this means that we can predict not just the median SOC estimate but also the 5th and 95th percentiles when the model is applied to each pixel. This again allowed us to estimate pixel-wise uncertainty and compute the uncertainty map.
We will add a brief explanation following the mentioning of the “quantreg = TRUE option” statement with a reference to paragraph 2.8 Model uncertainty:“We trained the final QRF used for spatial extrapolation of SOC stocks on the entire dataset with the quantreg = TRUE option in the ranger package (Wright & Ziegler, 2017), enabling the prediction of quantile-based uncertainty intervals during spatial extrapolation (see 2.8 Model uncertainty)."
L276: Did you exclude the pixels out of AOA from your analysis?
We did not exclude the pixel as we consider this as a decision made by the user of the map reader of the map.
L294: Did you predict SOC in forests?
We did not. This is already explicitly mentioned in line 302: “Sites located within forested areas (n = 21) were excluded from the analysis, and areas assigned to “forest” were removed from the final SOC stock map.”
L298: I am not familiar with such studies, but did you see: https://doi.org/10.21425/F5FBG61746?
We were not aware of this study. Thank you for pointing this out. If such a product existed for the Swiss Alps, we would certainly integrate it in the separation of land cover areas and the calculation of stocks per-area. We have checked with experts, and, unfortunately, such products do not yet exist here.
We will change line 298 to “Heath vegetation above the treeline was included within the alpine grassland class, as no land-cover product is currently available that reliably distinguishes ericaceous shrub belts in Swiss alpine grasslands (but see Bayle et al., 2024).”
We will also address this knowledge gap in the outlook, i.e. we will add the following in line 614:
“Shrub encroachment poses an ongoing challenge in the management of alpine ecosystems with widespread implications for above-belowground processes (Loarden et al; 2025; Marchal et al. 2025; Broadbent et al., 2022). Litter from ericaceous shrubs, rich in tannins and phenolic compounds, decomposes slowly and is associated with elevated SOC stocks (Gavazov et al., 2022), yet the spatial separation of heath- and grassland remains challenging (Bayle et al., 2024). Explicitly mapping alpine heathland communities and integrating them into regional SOC assessments would allow a better quantification of their contribution to total carbon stocks and their role in the context of ongoing greening and shrubification."
-Gavazov, K., Canarini, A., Jassey, V. E. J., Mills, R., Richter, A., Sundqvist, M. K., Väisänen, M., Walker, T. W. N., Wardle, D. A., & Dorrepaal, E. (2022). Plant-microbial linkages underpin carbon sequestration in contrasting mountain tundra vegetation types. Soil Biology and Biochemistry, 165. https://doi.org/10.1016/j.soilbio.2021.108530
-Laorden, L., Karl, C., Elena, G., García, T., Lyonnard, B., Pascale, M., Christiane, C., Tappeiner, U., Leitinger, G., & Lavorel, S. (2025). Shrub encroachment modifies soil properties through plant resource economics traits. Plant and Soil, 514(2), 2083–2104. https://doi.org/10.1007/s11104-025-07506-3
-Marchal, C., Tello-García, E., Laorden-Camacho, L., Binet, M.-N., Grigulis, K., Colace, M.-P., Périgon, S., Arnoldi, C., Rioux, D., Miquel, C., Laporte, F., Gallet, C., Gonindard-Melodelima, C., Leitinger, G., Mouhamadou, B., & Lavorel, S. (2025). Ericaceous shrub encroachment influences the structure of soil fungal communities in subalpine grasslands. Applied Soil Ecology, 208, 105985. https://doi.org/https://doi.org/10.1016/j.apsoil.2025.105985
-Broadbent, A. A. D., Bahn, M., Pritchard, W. J., Newbold, L. K., Goodall, T., Guinta, A., Snell, H. S. K., Cordero, I., Michas, A., Grant, H. K., Soto, D. X., Kaufmann, R., Schloter, M., Griffiths, R. I., & Bardgett, R. D. (2022). Shrub expansion modulates belowground impacts of changing snow conditions in alpine grasslands. Ecology Letters, May 2021, 52–64. https://doi.org/10.1111/ele.13903L300: This pseudo-zero method needs further clarification. How did you choose the site? Did you attribute a stock of 0?
As written in line 301, we chose sites on ice-covered ares and bare, exposed rock areas where we are certain that there is 0 SOC based on aerial images.
L303: Did the forest soil data be used for model calibration?
We are not sure if we understand this comment of referee 1. However, our approach focuses on areas above the treeline. as written in line 97 (study area description), line 302 (removal of sites in forests), and line 303 (removal of forested areas in the SOC map), as well as in the title.
L324: No comma or conjunction after 0.49
We will adjust this.
L329: “negligable” → “negligible”)
We will adjust this.
L329: “among [the covariates]”?
Thanks for catching this. We will adjust this!
L377: Can you give in supplementary the elevational density of your samples? Some elevations may be underrepresented in your dataset.
We suggest adding a density histogram of the elevations covered by our samples. However, note that appendix figure A8 shows the elevational distribution of the sites already.
L399: The discussion needs more in-depth mechanisms explanation and not just a repetition of the results with a comparison of other studies.
We thank the reviewer for this suggestion. However, digital soil mapping studies typically refrain from in-depth mechanistic interpretation, as the underlying process-level data are rarely available at spatial scales. We believe our current discussion already approaches the speculative limit of what the data can support. For a detailed treatment of soil chemical controls on SOC, we refer the reviewer to Udke et al. (in review, Global Change Biology), which addresses these mechanisms explicitly with additional soil-specific data collected at profile sites used for PTF development in our manuscript. However, our manuscript here is focused on delivering a spatial product, which inherently constrains the depth of mechanistic discussion.
L403: Can you give an estimated value of SOC stocks at national scale for comparison?
We thank the referee for the suggestions. We suggest moving some of the numbers that we have provided later on in line 442 comparing the total stock in alpine regions with agricultural and forested regions. However, since no study has quantified the SOC stock in alpine regions (line X), nor mapped all regions simultaneously (due to data scarcity in alpine regions), there is no number that quantifies the total SOC stock of Switzerland.
L415: What do you [they] mean by “managed” grassland?
We suggest clarifying this, by adding that managed grassland refers to extensive livestock grazing by cows or sheep.
L419: Yes, illustrating the crucial role to integrate the coarse fragments content, especially rocks and stone.
L425: And what about petrogenic organic matter?
We suggest adding a specific note on the fact that our study does not include petrogenic organic matter in the methods section and a note on that in the discussion.
L432-436: Did you try to compare the model results in high elevation with and without integrating pseudo-zero samples?
See our explanation to major comment nr. 3.
L451: Landolt-R is not really an estimation of the parent material properties.
It is the best spatially available data for the Swiss Alps that approximates pH conditions in the topsoil - which we discussed in the manuscript. See our explanation to major comment nr. 2.
L451-460: What are the uncertainties of Landolt-R and Landolt-H variables? Did you integrate these uncertainties into your model?
See our explanation to major comment nr. 2. We would like to highlight that the spatial layers of Descombes et al (2020) are based on 6 Mio georeferenced plant observation (on the relatively small area of Switzerland) and can thus be considered very robust.
L471: Are there any bog and mire samples in your training dataset?
Yes, there are indeed a few samples that happened to be fens or located in depressions, where at least historically there were mires or fens. If this is preferred, we can add this in the sampling description. We don’t think this is needed as this would implicitly require a more detailed site description or classification of all 307 sites, raising further issues.
L497: pH conditions on carbonate bedrock could also limit microbial decomposition if they are too alkaline. See “optimal windows” from https://doi.org/10.1007/s10533-017-0410-1
We agree that this is an interesting point. In fact, we referred to this so called “window of opportunity” described Rowley et al (2019) in a previous version of our manuscript. We can surely add it again.
L509: yes, but plant communities also depend of SOC content…
We agree, and this is not in disagreement with our manuscript. Notably, L509 refers to the reaction value (Landolt-R) and not the humus value (Landolt-H).
L515: Beware of circular reasoning: you predicted SOC using a variable that indicates high SOC content…
We agree that the use of a humus indicator is not a mechanistic predictor of SOC (see also our response to the major comment nr. 2). The line 515 mentioned by referee 1 is followed by the sentence: “Accordingly, the high rank in the importance metrics reflects shared variation with SOC rather than an independent mechanistic driver or a controlling factor of SOC” in line 518. We do not see any circular reasoning here. It is explicitly labelled and referred to as “covariance”, aligning with the concerns of referee 1.
L530: And also the uncertainties due to the covariates themselves.
This is a shared limitation of digital soil mapping studies. We suggest explicitly mentioning this additionally in the section on uncertainties. See also our answer to the major comment nr 4.
L549: Yes but this estimation remains highly uncertain, you didn’t estimate the soil depth.
We are not sure what this comment points to. However, see our response to the major comment nr. 1 where we address the point about how soil depth was estimated.
L560: The spatial distribution of residuals was only assessed through visual representation. Did you use Moran’s I test to assess spatial autocorrelation?
We suggest implementing Moran’s I test and thank the referee for the suggestion.
L568: Did you sample bogs and fens in your study?
Not specifically, but yes, some Sampling points coincidentally fell into areas that were at least historically a mire / fen. We would classify at least five sites as natural fens based on the vegetation composition. As described in our sampling description, the choice of sites was not based on presence / absence of specific habitats or soil types.
L595: This is only true for upper grasslands. What about grasslands just above the treeline?
The statement refers to Appendix Figure A6, referenced in the preceding sentence. We suggest adding in L595 that for grasslands around the treeline, SD may reach up to 5 kg/m2 and refer to Appendix Figure A6.
L599: I agree but the time scale between soil formation and SOC accumulation could be very different.
We suggest acknowledging this explicitly by adding the additional sentence in italic:
“These changes may facilitate soil formation and SOC accumulation through increased biomass production and C input in high alpine regions that further promotes colonization by upward shifting vegetation – analogous to SOC formation along retreating glaciers or debris flow fans (Egli et al., 2010; Y. Zhang et al., 2023). However, the timescales over which soil formation and SOC accumulation respond to ongoing greening remain highly uncertain and are difficult to project from current evidence (but see Juselius et al, 2022). "
Juselius, T., Ravolainen, V., Zhang, H. et al. Newly initiated carbon stock, organic soil accumulation patterns and main driving factors in the High Arctic Svalbard, Norway. Sci Rep 12, 4679 (2022). https://doi.org/10.1038/s41598-022-08652-9
L609: You mean that MAOC could be saturated in alpine soils? This seems rather unlikely, given that such results are rarely observed at low elevations.
We agree with the reviewer that MAOC is unlikely saturated in alpine soils and rather refer to a higher degree of C saturation relative to iron and aluminum oxides and hydroxides in the soil. We will clarify this in the section.
L622: The conclusion is concise, clear and summarizes well the main findings of the study.
Figure 1: Why are “forest” and “water/ice” include as cover class although they are not analyzed in the study?
This is described in paragraph “2.9 Spatial mask and landcover classes”, but we suggest further specifying in the figure caption that the forest coverclass is necessary to crop out forested areas (mentioned line 304), and similarly, water/ice covered sites are excluded from the final SOC map, as there is arguably no soil and thus no SOC stored in these areas.
Figure 2: It is debatable and not common to predict the rock fraction using pedotransfer function. The relationship between coarse content and depth is not obvious and it is dependent on the coarse fragment considered (stone, blocks, pebbles, …).
We agree that the approach is uncommon. However, we found relatively strong relationships between depth and the coarse rock fragment fraction for the 17 soil profiles and consider it as the only reasonable approach based on the data on alpine soils available.
Figure 6: SOC decreases the soil pH, so the more the SOC content, the lower the soil pH. Thus, SOC could also explain the soil pH and the vegetation preference, and not the opposite.
The partial dependence plots in figure 6 shows the marginal statistical association between Landolt-R and predicted SOC within the random forest model, not a causal direction. The key point is that Landolt-R is used as a spatial mapping tool capturing the calcareous vs. siliceous bedrock contrast at the broad landscape scale, where the dominant pH gradient is dependent on bedrock geology (and well reflected in Landolt-R). We will add a clarifying sentence to the Figure 6 caption noting that the partial dependence reflects a statistical association:
“Partial dependence plots reflect statistical associations within the QRF model and do not imply causal direction”
Citation: https://doi.org/10.5194/egusphere-2026-835-AC1
-
AC1: 'Reply on RC1', Michael Zehnder, 29 May 2026
-
RC2: 'Comment on egusphere-2026-835', Anonymous Referee #2, 07 Apr 2026
The manuscript “Quantifying soil organic carbon stocks above the treeline in the Swiss Alps” is a very important contribution for carbon stock estimations as it is the first comprehensive study that quantifies SOC stocks in the Swiss Alps above the treeline. The paper is well written and the methodology is scientifically sound using established methods and quantifies model uncertainty as well as uses the Area of Applicability to assess where the model is reliable based on the training data.
I agree with comments raised by R1 especially on the uncertainty of the usage of PTF for lower soil depths and also the availability of only 17 sites for PTF. While at the same time I would like to acknowledge that sampling in this terrain is challenging and labor intensive, I am just curious how the amount of 17 sites is justified?
In addition, I have the following concerns:
- I was wondering what the modelled depth for SOC stocks was. From the figures I see that it was mapped up to 80 cm depth. Did you build separate models for each depth increment (increments shown f.ex. in Figure 2 and 4) or calculated the total stock and built one model for the entire SOC stock? I assume accuracy for the upper soil (0-20 cm) would be higher than for the lower soil for which PTF are used. A map for the upper soil will likely have a higher accuracy than for the whole depth (0-80 cm) combined in one model.
- Was the number of sampling points for each landcover class (grasslands, scree, and bare rock) proportionate to their respective area fractions? Bare rock has only 11 sampling sites, is this representative?
Specific comments:
Line 33: “ with the majority being stored in the northern circumpolar regions” – Could you give a number?
Line 16: You mention “SOC stock map at 25 m resolution”. This is important information, but it would also be helpful to specify the soil depth considered
Line 21: You describe the variables as “highly informative” Since this assessment is based on variable importance, it would be clearer to state this explicitly
Line 64/65: “Therefore, distinguishing these three major land-cover types (alpine grasslands, partly vegetated scree, and bare rock) is essential for quantifying current SOC stocks.”
Out of curiosity, could separate models for each class have improved prediction results, given that these classes have fundamentally different SOC stocks? I realize that the sample size may make this infeasible, but it could be interesting to discuss briefly.Line 90: “addressing the questionS: …” The research questions are included in the text, but consider presenting them as numbered or bulleted points. This can improve clarity and make it easier to refer back to them in the Results and Conclusion sections
Line 134-135: This refers to the additionally sampled sites for IMIS. Consider moving this text after lines 14–17 or clarifying it further, as it was a bit confusing on first reading. It might also be helpful to include a table summarizing the total number of samples used in the study, distinguishing between those from existing datasets and those collected additionally
Line 138: I am not sure if I did understand “subsoil” in the context of this study correctly. I guess you define the soil below 20 cm as subsoil?
Line 194-195: Consider including McBratney et al. (2003) and the “SCORPAN” model. It would be helpful to indicate in Table A1 which variables correspond to which soil-forming factors and clarify why the older concept by Jenny was chosen/ mentioned instead of the newer one by McBratney.
Line 300: Could you clarify why 11 pseudo-zero points were selected? Are there any guidelines for this, and how does the number of pseudo-zeros relate to the total number of sampling points?
Line 319: since you used PTF for estimations below 20cm- how high is the uncertainty associated to the estimations < 20 cm depth?
Line 450 - 451: How is the factor “time” considered? I see it is discussed as not being adequately represented, but it may also be worth mentioning it in the Methods section already.
Line 624: at 25 m spatial resolution
Line 626 and elsewhere: Consider replacing “informative” with “important” or referring to the “variable importance” to clearly indicate that the importance is derived from the model
Conclusion: Ensure that the conclusion refers back to the research questions posed in the Introduction and addresses them fully. For example, consider including results for all landcover classes.
Line 1000: “..with an_ 8-cell smoothing”
Overall, the authors present a robust and well-executed study given the inherent data constraints; however, the manuscript would benefit from a more explicit discussion of methodological limitations and associated uncertainties, especially propagation of errors.
Citation: https://doi.org/10.5194/egusphere-2026-835-RC2 -
AC2: 'Reply on RC2', Michael Zehnder, 29 May 2026
Dear anonymous reviewer,
We sincerely appreciate the constructive feedback on our manuscript. Should we be given the opportunity to revise, we are confident we can address all concerns raised and substantially improve the manuscript. Below, we outline our proposed responses to each comment.
Kind regards on behalf of all listed authors,Michael Zehnder
The manuscript “Quantifying soil organic carbon stocks above the treeline in the Swiss Alps” is a very important contribution for carbon stock estimations as it is the first comprehensive study that quantifies SOC stocks in the Swiss Alps above the treeline. The paper is well written and the methodology is scientifically sound using established methods and quantifies model uncertainty as well as uses the Area of Applicability to assess where the model is reliable based on the training data.
We thank the referee for the overall positive feedback.
General comments
I agree with comments raised by R1 especially on the uncertainty of the usage of PTF for lower soil depths and also the availability of only 17 sites for PTF. While at the same time I would like to acknowledge that sampling in this terrain is challenging and labor intensive, I am just curious how the amount of 17 sites is justified?
We thank Referee 2 for the positive and constructive feedback. While the choice of 17 sites may initially appear arbitrary, there is a logical explanation, and we will expand on this in the main text. The soil profiles used to develop the PTFs were excavated along three elevational transects: two on siliceous bedrock and one on calcareous bedrock (Udke et al., in review, Global Change Biology; additional publications with these profiles are in prep.). Profiles were sampled at approximately 175 m elevational intervals from the treeline to the summit, with the aim of obtaining ~six profiles per transect. This sampling design was largely achieved. However, at the summit of the calcareous transect, it was virtually impossible to excavate a soil profile beyond 5 cm depth (Appendix B, line 1126), as the site consisted almost entirely of rock. As a result, the final dataset comprised 17 rather than 18 actual soil profiles. To account for this gap, we created a pseudo-profile for this specific site using information from the high-elevation site on the siliceous transect, with the rationale explained in the lines following line 1126. We continue to refer to the true count of 17 sampled sites, as we considered this the more transparent and accurate representation of the field effort. However, we would like to emphasize that the additional pseudo-profile improved the practical implementation of the workflow and was needed to avoid bias in the PTF workflow (line 1130).
In addition, I have the following concerns:
• I was wondering what the modelled depth for SOC stocks was. From the figures I see that it was mapped up to 80 cm depth. Did you build separate models for each depth increment (increments shown f.ex. in Figure 2 and 4) or calculated the total stock and built one model for the entire SOC stock? I assume accuracy for the upper soil (0-20 cm) would be higher than for the lower soil for which PTF are used. A map for the upper soil will likely have a higher accuracy than for the whole depth (0-80 cm) combined in one model.We built one spatial model, based on the total SOC stocks of the sites, including PTF-modelled subsoil stocks down to bedrock. We will expand on this in the method section, see below (and see also our comments to referee nr. 1 for concrete suggestions).
We calculated stocks for the depth compartments 0-5, 5-10, and 10-20 cm based on the measured values of the variables needed for calculation of SOC stocks. Below 20 cm, we calculated them for 1cm increments, i.e. we modelled FE density, SOC% and the coarse rock fragment contents for 20.5, 21.5, 22.5, … 79.5 cm midpoints (based on the PTF for each variable) and used 1cm depth intervals. This is more appropriate given the exponential decline of SOC% with increasing depth. This procedure is described in line 180 (“SOC stocks for the subsoil layers were then calculated at 1 cm resolution according to Eq. 1, vertically aggregated within 10 cm depth intervals, and added to the measured topsoil stocks.”) in the methods section, but we understand that we may have to elaborate on this a bit more. The aggregation of 1 cm increments to 10 cm intervals mentioned in line 180 was primarily done for visual purpose (e.g. Figure 2 and Figure 4b). For the total SOC stocks per site we then added up the measured 0 – 20cm stocks to the modelled 1 cm SOC stocks from 20-80 cm as the referee has correctly understood (although very often PTF-modelled SOC stocks below ~40cm are approximating 0 since either SOC% approaches 0 or coarse rock fragment fraction reaches 100 % i.e. parent material). We suggest clarifying this further in the methods section.
In that sense, accuracy for the measured 20 cm topsoil is more accurate. As we suggested in response to the first referee’s comments, we add an additional map based on only measured topsoil stocks. However, we would like to emphasize that the use of PTFs is particularly beneficial for estimating total SOC stocks in the lower alpine belt, where restricting measurements to the upper 20 cm of soil would often substantially underestimate total SOC storage due to the much greater soil depths > 20 cm commonly found in these areas.
Was the number of sampling points for each landcover class (grasslands, scree, and bare rock) proportionate to their respective area fractions? Bare rock has only 11 sampling sites, is this representative?
We thank the referee for this thoughtful comment. We reported the area fractions in line 362, which are approximately 2:1:1 for grassland:scree:bare rock. We agree that these proportions are not well reflected in our sampling dataset, where the 307 observations were distributed as 259:37:11 across these land-cover classes. However, random forests do not require predictors to be sampled in proportion to their natural prevalence, though they rely on sufficient representation of the environmental space to learn reliable predictor–response relationships. For this reason, we introduced pseudo-zero observations with SOC values of “0” in underrepresented SOC-free areas (see also our response to major comment 3 from Referee 1). While this still does not perfectly match the actual land-cover proportions, this does not inherently bias the random forest model. Instead, the pseudo-zeros helped to ensure that truly SOC-free environments were represented in the training space and prevented systematic overprediction in bare rock areas (based on the minimum of what was measured on sites that fell in the “bare rock” cover class). We suggest elaborating on this in the manuscript.
Specific comments:
Line 33: “ with the majority being stored in the northern circumpolar regions” – Could you give a number?
We thank the referee for this hint. We agree that this calls for an explicit number, which we will happily provide from the cited literature (e.g. Hugelius et al., 2014; Mishra et al., 2022), which estimate SOC stocks at ~1000 Pg C (0-3 m).
Line 16: You mention “SOC stock map at 25 m resolution”. This is important information, but it would also be helpful to specify the soil depth considered
We will clarify this by writing: “SOC stocks map at 25 m resolution down to the bedrock“
Line 21: You describe the variables as “highly informative” Since this assessment is based on variable importance, it would be clearer to state this explicitly
We thank the referee for pointing this out. We had deliberately removed the explicit reference to “variable importance” because the term can be somewhat vague for readers who are not familiar with random forest modelling terminology (variable importance – “important for what”?). Our intention was to make the statement more accessible by referring more generally to variables as “highly informative.” However, we agree that explicitly linking this assessment to the variable importance metric improves clarity (at least for the random forest-literate audience), and we are happy to reintroduce this wording – possibly with an additional statement that introduces the term “node-impurity”
Line 64/65: “Therefore, distinguishing these three major land-cover types (alpine grasslands, partly vegetated scree, and bare rock) is essential for quantifying current SOC stocks.”
We thank the referee for catching this and we will implement this accordingly.
Out of curiosity, could separate models for each class have improved prediction results, given that these classes have fundamentally different SOC stocks? I realize that the sample size may make this infeasible, but it could be interesting to discuss briefly.This is an interesting suggestion. In principle, separate models for each land-cover class could be beneficial if predictor–response relationships differed fundamentally between classes and sufficient training data were available for each subset. However, quantile regression forests already partition the predictor space recursively and can implicitly learn class-specific relationships through variables such as elevation, NDVI, and vegetation indicators that implicitely divide the data into broad branches of the decision trees of the model. Given the relatively small sample sizes for scree and especially bare rock sites, splitting the dataset would likely reduce model stability and increase overfitting risk. We therefore opted for a single model, although we thought about this during the course of the project. The suggestion is very good and such an approach may be beneficial when modelling country-wide SOC stocks across even more habitats. We suggest discussing our choice briefly in the discussion.
Line 90: “addressing the questionS: …” The research questions are included in the text, but consider presenting them as numbered or bulleted points. This can improve clarity and make it easier to refer back to them in the Results and Conclusion sections
We thank the referee for the suggestion. We happily provide the research questions as bullet points if this is preferred. We don’t have a strong preference here.
Line 134-135: This refers to the additionally sampled sites for IMIS. Consider moving this text after lines 14–17 or clarifying it further, as it was a bit confusing on first reading. It might also be helpful to include a table summarizing the total number of samples used in the study, distinguishing between those from existing datasets and those collected additionally
We thank the referee for raising this point. We are open to moving this paragraph further up to line 314-317 (which we believe is meant by the referee). We opted for a separate paragraph as the profiles used for PTF development should not be confused with the topsoil data used for modelling. We will consider a table that summarizes the sampling design better and we thank the referee for the constructive suggestion.Line 138: I am not sure if I did understand “subsoil” in the context of this study correctly. I guess you define the soil below 20 cm as subsoil?
Yes, this is exactly what is meant. We suggest introducing the term explicitly at the point where it is first mentioned (line 138).
Line 194-195: Consider including McBratney et al. (2003) and the “SCORPAN” model. It would be helpful to indicate in Table A1 which variables correspond to which soil-forming factors and clarify why the older concept by Jenny was chosen/ mentioned instead of the newer one by McBratney.
We thank the referee for suggesting the use of McBratney et al. (2003). There is no strong conceptual reason for our choice, besides the fact that Jenny developed his soil forming factors in Switzerland (before going to the US). Actually, he developed them in Grisons where a lot of our sites are located, and MacBratney is based in Australia. Therefore, we suggest sticking with Jenny and indicating in Table A1 which soil forming factors each variable represents.
We thank the referee for the constructive suggestion in the Appendix Table A1. We will implement the suggestion accordingly.
Line 300: Could you clarify why 11 pseudo-zero points were selected? Are there any guidelines for this, and how does the number of pseudo-zeros relate to the total number of sampling points?
We thank the referee for raising this point. It was chosen to match the number of sites that fell in the category of “bare rock”.
The concept of pseudo-zeros is not new to digital soil mapping, e.g. Siewert (2018) https://doi.org/10.5194/bg-15-1663-2018 chose 10 pseudo-zeros (with 0 kg/m2 SOC) in combination with 47 pedons, while SoilGrids 250m from Hengl et al (2017) https://doi.org/10.1371/journal.pone.0169748 was trained with hundreds of pseudo-zeros in the Sahara or on ice-covered areas. The rationale was the same as ours; sites with 0 kg/m2 are normally not sampled. Hengl et al (2017) suggests 1–2% of the total number of training points, which we match quite well. We suggest adding this in the methods section in line 300 and citing Siewert (2018) and Hengl et al (2017) as an example as described in our comment to the first referee’s major comment nr. 3
"“While permanently snow- and ice-covered areas and completely exposed bedrock contain no soil and were therefore not sampled in the field, we added a set of expert-led pseudo-zero points (n = 11) located above 3150 m a.s.l. in permanently snow- and ice-covered areas and on bare rock to inform model behaviour in these high-elevation environments as suggested by Wadoux et al., (2020). Pseudo-zeros can be used to constrain model behavior in areas where SOC stocks can be assumed to be zero with high confidence, as applied for example in SoilGrids 250m (Hengl et al., 2017), where hundreds of pseudo-zero observations were used in desert and ice-covered regions.”
Line 319: since you used PTF for estimations below 20cm- how high is the uncertainty associated to the estimations < 20 cm depth?
We assume the referee refers to SOC estimates below (>) 20 cm depth. We agree that uncertainty associated with PTF-derived subsoil estimates is an important consideration.
The most robust way to quantify this uncertainty would be an independent validation using additional full soil profiles excavated adjacent to randomly selected topsoil sampling locations. SOC stocks could then be quantified directly down to bedrock (in the traditional way) and compared with the stocks predicted from topsoil measurements using our PTF approach. In fact, not the SOC stocks should be compared but the in-situ measured variables, i.e. fine earth densities, the SOC % and the coarse rock fraction, for depth intervals below 20 cm. However, such validation was beyond the scope of this study due to the logistical constraints of excavating additional full profiles in remote alpine terrain (a maximum of one site can be visited per day).In principle, we could also have performed a cross-validation exercise using the two elevational transects on comparable silicious bedrock by predicting SOC stocks below 20 cm in transect A from environmentally similar profiles in transect B (and vice versa). However, this would represent a highly conservative test because only half of the available deep soil profiles would be used for parameterization at any given location, likely leading to an overestimation of prediction uncertainty.
Instead, we provided a pseudo-validation of the PTF-derived subsoil estimates using the same sites for which full profile information was available (Appendix C, Figure C1a). We acknowledge that this validation is limited because it does not represent an independent test and may underestimate true uncertainty. This is mentioned in line 556 in the discussion of uncertainties.
As described above, we will add an additional map for only the top 20cm of soil with measured SOC stocks independent of PTF
Line 450 - 451: How is the factor “time” considered? I see it is discussed as not being adequately represented, but it may also be worth mentioning it in the Methods section already.We thank the referee for this thoughtful comment. We suggest adding a note on the fact that no spatial layers capturing the soil forming factor of time was available in the methods sections (where we introduce Jenni/SCORPAN), and adding a brief note on the potential usefulness of such a layer in the discussion, where we encourage efforts to produce a spatial pH layer for the alpine landscape (line 508).
Line 624: at 25 m spatial resolution
We thank the referee for catching this missing word.
Line 626 and elsewhere: Consider replacing “informative” with “important” or referring to the “variable importance” to clearly indicate that the importance is derived from the model
We thank the referee for pointing this out. We will consider the replacement but see also our comment above on the term “variable importance”.
Conclusion: Ensure that the conclusion refers back to the research questions posed in the Introduction and addresses them fully. For example, consider including results for all landcover classes.
We thank the referee for raising this point. We agree that all research questions should be covered in the conclusions and will refine this section accordingly.
Line 1000: “..with an_ 8-cell smoothing”
We thank the referee for catching this typo in the Appendix.
Overall, the authors present a robust and well-executed study given the inherent data constraints; however, the manuscript would benefit from a more explicit discussion of methodological limitations and associated uncertainties, especially propagation of errors.
We thank referee 2 for the constructive feedback and the overall positive reception of our manuscript. We will happily incorporate the suggestions.
Citation: https://doi.org/10.5194/egusphere-2026-835-AC2
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 956 | 762 | 75 | 1,793 | 159 | 109 | 146 |
- HTML: 956
- PDF: 762
- XML: 75
- Total: 1,793
- Supplement: 159
- BibTeX: 109
- EndNote: 146
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Overview
This study quantifies soil organic carbon stocks above the treeline in the Swiss Alps by combining 307 field sites with pedotransfer-based subsoil extrapolation and spatial modelling using Quantile Regression Forests. It estimates 47.6 Mt SOC stored across alpine terrain and identifies NDVI, vegetation-derived humus indicators, temperature, and topography as key predictors of SOC spatial variability.
Overall, the manuscript presents a valuable and timely contribution to SOC analysis in alpine ecosystems. It is well-written, provides an important baseline for SOC stocks and fits very well the scope of the journal. However, clarifications are needed before the manuscript could be considered for publication, especially regarding the robustness of PTF for subsoil extrapolation, the predictor circularity, the coarse fragment integration, the pseudo-zero weighting, and the full uncertainty propagation.
Main concerns
1- First, a major methodological concern relates to the use of pedotransfer functions (PTF) to estimate subsoil SOC stocks. Subsoil values are extrapolated from only 17 soil profiles stratified by bedrock type and elevation bands, yet these relationships are applied to 307 sampling sites across the alpine region. Given the well-known heterogeneity of alpine soils, including high skeletal content, cryoturbation, discontinuous soil depth, and highly variable pedogenic pathways, the robustness and transferability of these PTF require more rigorous justification and uncertainty assessment. Another option is to move all this analysis into supplementary information. In particular, the procedure of selecting the “best-fitting” profile for each site may introduce overfitting and underestimate the structural uncertainty associated with subsoil extrapolation. Related to this, the manuscript does not clearly explain how soil depth is estimated or constrained across the study area. Without explicit constraints on soil depth, extrapolation of SOC stocks below the measured 0–20 cm layer becomes difficult to evaluate. Moreover, soil types such as Umbrisols, Podzols, Rankers, Rendzinas, and Cambisols can occur over short spatial distances and exhibit very different vertical SOC distributions and bulk density patterns.
2- Second, Landolt-H and Landolt-R are partly circular predictors: both are vegetation-derived proxies likely co-varying with SOC. Their strong importance may reflect shared response rather than mechanistic control, potentially inflating predictive performance while limiting process inference. Moreover, the uncertainties of these variables are not integrated in the modelling approach, and I assume that these uncertainties are spatially distributed.
3- Third, the treatment of coarse fragments and bare-rock areas raises methodological concerns that may bias SOC stock estimates. The sampling protocol did not allow the quantification of coarse fragments larger than 5 cm, which are common in alpine soils. In addition, the treatment of bare-rock areas through the introduction of “pseudo-zero” SOC values based on expert judgement lacks sufficient methodological justification. This assumption directly influences SOC predictions in sparsely vegetated, high-elevation zones and may significantly affect regional SOC stock estimates. It is unclear why bare-rock areas were not explicitly excluded from SOC calculations, unless the study explicitly considers petrogenic organic carbon or organic matter stored in lithic substrates.
4- Finally, although the geostatistical framework appears robust and overall well implemented, the treatment of uncertainty remains incomplete. The reported total SOC stock estimate does not appear to fully propagate the different sources of uncertainty inherent to the methodology. In particular, uncertainties associated with pedotransfer functions, covariate uncertainties, and model structural uncertainty are not explicitly integrated into the final confidence intervals. In addition, the assessment of spatial autocorrelation in the model outputs appears to rely primarily on visual inspection and not on formal statistical evaluation (variogram analysis or Moran’s I statistics on residuals).
Specific comments
L11: “… of the total [organic] carbon…”, because in carbonated context, some of carbon is include in carbonates.
L21: Using soil pH as explanatory variable is contestable, as SOC content (organic acid) decreases soil pH.
L23-25: I agree with this statement for the upper alpine grasslands, but what about alpine grasslands near the treeline, containing the highest SOC stocks?
L35: “due to the difficulty of accessing as well the high heterogeneity” → “due to the difficulty of access as well as the high heterogeneity”.
L42: Greening of alpine area is only one consequence of climate change on carbon dynamics, influencing quantity. Other consequences could be the change in land cover and thus change in SOC chemical composition, influencing also its stability.
L58: “suppress” à “slow, reduce,…”
L75-76: The end of the sentence does not follow on from the preceding sentences.
L90: Why you mention the “time” factor above in the introduction (L53), and you didn’t include it in the analysis?
L97: Given the elevation of study sites, could some areas be recently deforested? e.g. at 1750 m, potentially leading to legacy effect on SOC stocks not integrated by your model.
L105: The mosaic of parent material leads to wider variations in soil conditions than soil pH only, with some of them that could influence SOC stocks (Fe, Al, Ca, clay, sand contents, etc.). Your study didn’t really integrate the “parent material” factor. Is the integration of geologic properties could increase model performance?
L107: Can you give here an idea of the study area surface?
L122: It is not clear why there is 20 sites from GLORIA. 6 summits x 4 cardinal directions = 24 sites?
L129: mean [annual] precipitations I supposed?
L140: It is always problematic to assess the bulk density. Given the available data of the three datasets you have, your approach seems to be the best, although the difference in soil volume could introduce bias on the coarse fragment content estimation. Because in IMIS dataset, coarse fragments >2cm were not integrated in bulk density assessment, while in BDM only coarse fragments >4,8cm (and in GLORIA for coarse fragments >10cm) were not integrated. And given the importance of such coarse fragments in alpine soils (as we can see on the soil profile pictures in appendix B), this omission could have great importance on SOC stock calculation. Did you try to compare soil coarse fragment content between the three datasets? And the bulk density?
L145: “replica” → “replicate”?
L170: In this paragraph, it is not clear how do you predict spatially the soil depth.
L171: It is not clear for me which PTF you assign for each pixel of your study area. If I have well understood your study, you first spatialized SOC stocks in the upper 0-20cm horizon using QRF model, then you used PTF to assess SOC stocks for the whole profile, right? First, I suggest adding the use of PTF in the workflow of figure 1. Second, the stratification by bedrock type and elevation is debatable. Are these two parameters the only ones that control soil properties and stratification? Coarse content is more dependent on geomorphic processes rather than bedrock and elevation.
L177: Why did you choose linear regression for FE density while you used exponential decline for SOC content?
L218: “humus-indictor map” → “humus-indicator map”.
L260: What does the activation of this option mean?
L276: Did you exclude the pixels out of AOA from your analysis?
L294: Did you predict SOC in forests?
L298: I am not familiar with such studies, but did you see: https://doi.org/10.21425/F5FBG61746 ?
L300: This pseudo-zero method needs further clarification. How did you choose the site? Did you attribute a stock of 0?
L303: Did the forest soil data be used for model calibration?
L324: No comma or conjunction after 0.49
L329: “negligable” → “negligible”)
L329: “among [the covariates]”?
L377: Can you give in supplementary the elevational density of your samples? Some elevations may be underrepresented in your dataset.
L399: The discussion needs more in-depth mechanisms explanation and not just a repetition of the results with a comparison of other studies.
L403: Can you give an estimated value of SOC stocks at national scale for comparison?
L415: What do you [they] mean by “managed” grassland?
L419: Yes, illustrating the crucial role to integrate the coarse fragments content, especially rocks and stone.
L425: And what about petrogenic organic matter?
L432-436: Did you try to compare the model results in high elevation with and without integrating pseudo-zero samples?
L451: Landolt-R is not really an estimation of the parent material properties.
L451-460: What are the uncertainties of Landolt-R and Landolt-H variables? Did you integrate these uncertainties into your model?
L471: Are there any bog and mire samples in your training dataset?
L497: pH conditions on carbonate bedrock could also limit microbial decomposition if they are too alkaline. See “optimal windows” from https://doi.org/10.1007/s10533-017-0410-1
L509: yes, but plant communities also depend of SOC content…
L515: Beware of circular reasoning: you predicted SOC using a variable that indicates high SOC content…
L530: And also the uncertainties due to the covariates themselves.
L549: Yes but this estimation remains highly uncertain, you didn’t estimate the soil depth.
L560: The spatial distribution of residuals was only assessed through visual representation. Did you use Moran’s I test to assess spatial autocorrelation?
L568: Did you sample bogs and fens in your study?
L595: This is only true for upper grasslands. What about grasslands just above the treeline?
L599: I agree but the time scale between soil formation and SOC accumulation could be very different.
L609: You mean that MAOC could be saturated in alpine soils? This seems rather unlikely, given that such results are rarely observed at low elevations.
L622: The conclusion is concise, clear and summarizes well the main findings of the study.
Figure 1: Why are “forest” and “water/ice” include as cover class although they are not analyzed in the study?
Figure 2: It is debatable and not common to predict the rock fraction using pedotransfer function. The relationship between coarse content and depth is not obvious and it is dependent on the coarse fragment considered (stone, blocks, pebbles, …).
Figure 6: SOC decreases the soil pH, so the more the SOC content, the lower the soil pH. Thus, SOC could also explain the soil pH and the vegetation preference, and not the opposite.