002-reply-main.Rmd

---
title: 'Comment on Hodgkins et al. (2018)'
author: "Henning Teickner$^{1, *}$, Klaus-Holger Knorr$^1$"
date: ""
output: 
  bookdown::pdf_book:
    keep_tex: true
    toc: false
bibliography: references.bib
header-includes:
   - \usepackage{float}
   - \usepackage{lineno}
   - \linenumbers
   - \usepackage{setspace}
   - \doublespacing
   - \usepackage{titling}
   - \pretitle{\begin{center}
     \bf{This manuscript is a non-peer reviewed preprint submitted to EarthArXiv. It has not been submitted for publication to a journal yet.} \LARGE\\
     \vspace{5ex}}
   - \posttitle{\end{center}}
fontsize: 12pt
---

```{r setup-reply-2, include=FALSE}
knitr::opts_chunk$set(fig.align = "center", 
                      fig.pos = "H",
                      echo = FALSE,
                      cache = TRUE)
```

$^1$ ILÖK, Ecohydrology and Biogeochemistry Group, University of Münster, Heisenbergstr. 2, 48149 Münster, Germany  
$^*$ Corresponding author: Henning Teickner (henning.teickner@uni-muenster.de)  
 
\newpage

# Introduction {.unnumbered}
 
@Hodgkins.2018 used mid infrared spectra (MIRS) to make inferences about the stability of peat against decomposition along a latitudinal gradient from the tropics to northern latitudes. Recently, we have shown that these spectral prediction models are biased and their predictions more uncertain than considered in @Hodgkins.2018 [@Teickner.2022f]. Here, we show what consequences this bias and uncertainty and additional neglected uncertainty sources have for the main conclusions in @Hodgkins.2018.

In particular, we find that:

1. Larger aromatic contents may be necessary to stabilize tropical near-surface peat against decomposition and aromatics may accumulate at a higher rate with depth in tropical peatlands than estimated by @Hodgkins.2018.

2. Larger uncertainties indicate that also larger differences in aromatic contents between (sub)tropical and high latitude peat than estimated by @Hodgkins.2018 are possible, also between deeper peat.

3. More research should focus on how estimates of carbohydrate and aromatic contents from MIRS may be confounded by other organic matter fractions, in particular proteins. As a first step, this requires accurate concepts to name organic matter fractions and variables used in the interpretation of MIRS.

These uncertainties have the potential to change predictions of models on global peat stability as well as what stabilizes peat against decomposition if temperatures rise and thus add in particular to the debate to what extent global warming may increase decomposition of northern deep peat deposits (e.g. @Dorrepaal.2009, @Wilson.2016, @Baysinger.2022).

That said, theoretical considerations of decomposition processes alone support the suggested gradient in peat chemistry and @Hodgkins.2018 made an important contribution in providing the first open access models and estimates for peat carbohydrate and aromatics content on a large scale. Increasing the accuracy of models is an important part of the scientific process and we hope that this comment will focus research on the problems which hamper improving this accuracy.

In the following, we will use the terms holocellulose instead of the more general term carbohydrates and Klason lignin instead of the more general term aromatics. These are the accurate terms for the organic matter fractions quantified in @Hodgkins.2018.


# Uncertainty sources not considered in @Hodgkins.2018 {.unnumbered}

In a recent study, we have shown that the models used to estimate holocellulose and Klason lignin contents in @Hodgkins.2018 (original models) are not valid for peat and this will cause larger uncertainties and biased predictions (uncertainty source 1) --- especially for decomposed peat --- of the estimated holocellulose and Klason lignin contents [@Teickner.2022f].

Besides this, the following two uncertainty sources have not been considered in @Hodgkins.2018: first, the prediction uncertainty of the spectral prediction models (uncertainty source 2), and second, the uncertainty introduced when computing near-surface average holocellulose and Klason lignin contents from estimates for individual peat layers (uncertainty source 3).

We recomputed the models for the latitudinal gradients of near-surface peat holocellulose and Klason lignin contents while considering all three uncertainty sources (supporting information S1). 

Since we do not have yet models to accurately predict peat holocellulose and Klason lignin contents from MIRS, it is of course difficult to quantify the uncertainty and bias introduced by uncertainty source 1. However, the modified models provided in @Teickner.2022f are very likely more accurate and therefore the difference between predictions of the original model and our modified models are a plausible approximation of this additional uncertainty and bias. 

The reanalysis shows that both the estimated slopes for the latitudinal gradient (95% confidence intervals are [`r s1_res_holocellulose %>% dplyr::pull(b_latitude) %>% quantile(probs = c(0.025, 0.975)) %>% round(1) %>% paste0(collapse = ",")`] and [`r s1_res_klason_lignin %>% dplyr::pull(b_latitude) %>% quantile(probs = c(0.025, 0.975)) %>% round(1) %>% paste0(collapse = ",")`] for holocellulose and Klason lignin, respectively; figure \@ref(fig:reply-res-p-sites-latitude) and supporting figure 3<!---\@ref(fig:reply-res-p-latitude-poc-slope)--->) and differences in the depth profiles are highly uncertain due to these uncertainty sources (figure \@ref(fig:reply-res-p-depth-profiles-latitude)). 

Thus, our reanalysis shows that estimates of the latitudinal gradient and depth gradients are much more uncertain than previously assumed.  

(ref:reply-res-p-sites-latitude-caption) Predicted average surface ($\le 50$ cm) peat holocellulose (a) and Klason lignin (b) contents plotted against latitude (compare with Fig. 3 in @Hodgkins.2018). Lines and shaded areas represent average predictions from regression models and 95% confidence intervals. "Modified" is the modified Bayesian hierarchical regression model which simultaneously models individual samples' contents from mid infrared spectra and the latitudinal gradient of average core near-surface peat contents. This model considers prediction uncertainty from the mid infrared spectra and from computing per-core averages. "Original" is the original linear regression model [@Hodgkins.2018] computing only the latitudinal gradient of average core near-surface contents. This model does not consider prediction uncertainty from mid infrared spectra, nor uncertainty from computing per-core averages. Points are average core near-surface contents predicted from the model ("Modified") or computed from the average predictions for individual samples ("Original") with error bars representing 95% confidence intervals. Points for "Modified" are shifted by $+0.1^\circ$.

```{r reply-res-p-sites-latitude, echo=FALSE, out.width="70%", fig.height=3.5, fig.width=6.5, fig.cap='(ref:reply-res-p-sites-latitude-caption)'}
p_poc_latitude6
```

(ref:reply-res-p-depth-profiles-latitude-caption) Predicted median holocellulose (A) and Klason lignin (B) depth profiles of peat core data classified into two latitude categories following @Hodgkins.2018 (compare with Fig. 2 in @Hodgkins.2018). Lines are averages of LOESS smoothers fitted to the predicted values by the models. Shaded regions are corresponding 95% confidence intervals, comprising prediction uncertainty in the holocellulose and Klason lignin estimates, respectively. Vertical dashed lines in each columns represent approximate surface and average mean contents as predicted using the improved models ("Modified").

```{r reply-res-p-depth-profiles-latitude, echo=FALSE, warning=FALSE, out.width="70%", fig.height=5.5, fig.width=6.5, fig.cap='(ref:reply-res-p-depth-profiles-latitude-caption)'}
p_depth_profiles_latitude
```

# Larger aromatic contents may be necessary to stabilize tropical near-surface peat against decomposition {.unnumbered}

Our previous study has shown that the spectral prediction model for Klason lignin is biased, especially for more decomposed peat [@Teickner.2022f]. Our reanalysis using the modified model for Klason lignin from @Teickner.2022f indicates that larger Klason lignin contents may be necessary to stabilize tropical near-surface peat against decomposition as well as that changes in Klason lignin with depth may be more pronounced in the tropics than estimated in @Hodgkins.2018 (supporting information S1).

With our modified model, near-surface peat Klason lignin contents are on average ~10 to 15 mass-% larger across the latitudinal gradient (figure \@ref(fig:reply-res-p-sites-latitude)). Specifically, average (sub)tropical ($<45 ^\circ$N) near-surface peat Klason lignin contents are `r d_peat_pred_best_diff_mean %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude & surface_peat) %>% dplyr::pull(y_hat_diff) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude & surface_peat) %>% dplyr::pull(y_hat_diff) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-% (median, lower and upper 95% prediction interval limit) larger than with the original model. The large uncertainties now made explicit would also allow on average larger differences between deep (sub)tropical ($<45 ^\circ$N) and high-latitude ($\ge45 ^\circ$N) peat (figure \@ref(fig:reply-res-p-depth-profiles-latitude)).

Similarly, residual enrichment of Klason lignin during decomposition may have been underestimated. A rough estimate for the residual enrichment of Klason lignin during decomposition is the difference in Klason lignin content between near-surface peat and deeper peat. With the modified model, this difference is on average for (sub)tropical peatlands `r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude) %>% dplyr::pull(res) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & !high_latitude) %>% dplyr::pull(res) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-% larger (figure \@ref(fig:reply-res-p-depth-profiles-latitude)). For high latitude peatlands the difference is smaller and much more uncertain than previously stated (`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & high_latitude) %>% dplyr::pull(res) %>% median() %>% magrittr::multiply_by(100) %>% round(0)` [`r d_peat_pred_best_diff_mean1 %>% dplyr::mutate(res = y_hat_diff - y_hat_or_diff) %>% dplyr::filter(variable_y == "Klason lignin" & high_latitude) %>% dplyr::pull(res) %>% quantile(probs = c(0.025, 0.975)) %>% magrittr::multiply_by(100) %>% round(0) %>% paste(collapse = ",")`] mass-%).  

Consequently, in general --- and especially in (sub)tropical peatlands --- the average residual enrichment of Klason lignin due to decomposition probably has been underestimated by the original model. A consequence of this is that high latitude peat deposits may experience more decomposition under a warmer climate than suggested in @Hodgkins.2018 because a larger content of Klason lignin is necessary to stabilize peat chemically under warmer conditions. 


# The need to use precise concepts for organic matter fractions {.unnumbered}

We argue that we should differentiate between vague concepts such as carbohydrates and aromatics and precise concepts such as holocellulose and Klason lignin and we should always use the most precise concept possible to describe the organic matter fraction we _plan_ (or intent) to measure, unless more general statements are explicitly warranted.

For example, @Hodgkins.2018 planned to measure Klason lignin contents because their spectral prediction model used Klason lignin data as dependent variable, however this variable is labeled as aromatics. Limitations of the extraction procedure by which Klason lignin are defined are known (e.g., @Hatfield.2005, @Bunzel.2011, @Abu-Omar.2021) and these limitations as well as differences to other procedures get obscured by using vague words such as carbohydrates and aromatics.

However, using precise concepts for organic matter fractions is also important for exactly the opposite reason: to make clear how much accuracy and precision a variable actually has, instead of implying that it would quantify a variable more accurate and precise than is actually the case. This becomes especially important when interpreting peat chemistry based on mid infrared spectra (MIRS) or spectral prediction models because all variables derived from MIRS --- for example the widely used humification indices, as proxy for relative contents of recalcitrant organic matter fractions (e.g. @Broder.2012), or Klason lignin contents predicted from models --- are in fact only indirect estimates of organic matter fractions. These indirect estimates of organic matter fractions can be misleading in abundant ways if they are not sufficiently validated. Stating these limitations requires precise concepts.

Intensities recorded in MIRS really only represent the fraction of the incident infrared radiation which is absorbed by specific molecular structures which happen to absorb at that specific energy level [@Stuart.2004]. Such molecular structures may for example be aromatic C=C bonds, C-N bonds in amides, or C-O bonds in alcohols (e.g. carbohydrates) [@Stuart.2004]. When molecular structures absorb infrared radiation, this causes a change in dipole moment of one or more of their bonds and the larger the change in dipole moment, the more intense is the absorption and hence the larger the peak in a MIRS. The same stretching results in a larger change in dipole moment for more electronegative bonds than less electronegative bonds [@Stuart.2004]. Since C-N bonds are more electronegative than aromatic C=C bonds, this means that already a small amount of proteins in peat can contribute equally large or more to the peaks in MIRS around 1510 and 1630 cm$^{-1}$ than aromatic C=C bonds in the same region of MIRS (see Fig. 3 in @Reuter.2020 for an example). These confounding factors make it necessary to explicitly define the conditions under which a spectral prediction model or other variable derived from MIRS (e.g. a humification index) is a valid proxy for a specific, precisely defined, organic matter fraction.

One major reason why the model for Klason lignin in @Hodgkins.2018 is biased and invalid for peat is that the only predictor variable used in this model, `arom15arom16`, is not a good indicator for Klason lignin because, as explained above, also proteins can absorb in the same energy range [@Stuart.2004] and therefore, estimates are biased depending on the amount of proteins in the samples [@Teickner.2022f]. Therefore, we should interpret `arom15arom16` and similar variables, such as specific humification indices (HI$_{1630/1090}$, HI$_{1510/1090}$) only then as good proxies for the relative abundance of aromatic C=C bonds if we have shown that protein contents do not differ much between the peat samples (thresholds which still need to be established). 

In all other cases, where we cannot validate if a variable is a good proxy for a specific organic matter fraction, we should be precise in our wording by calling the variable `arom15arom16` (instead of Klason lignin or aromatics) to signal that this variable may be no good indicator for aromatics if it is confounded by proteins, i.e. we may agree on `arom15arom16` as precise name for this variable, but it must always be understood that it is first and foremost _only_ defined as sum of the area of two peaks extracted by the specific procedure proposed in @Hodgkins.2018 and nothing more. Obviously, this recommendation equally applies to all other variables derived from MIRS.

@Hodgkins.2018 have made an important first step to actually quantify holocellulose and Klason lignin contents, i.e. specific concepts of carbohydrates and aromatics, from MIRS and this is an important improvement over the qualitative interpretation of peak heights or peak ratios (such as humification indices) used in the past (e.g. @Cocozza.2003, @Broder.2012, @Tfaily.2014). Using precise words for what we have actually measured and what confounding factors we have considered is important to assure that this first step will actually be an improvement over the qualitative interpretation of MIRS.


# Conclusions {.unnumbered}

Our reanalysis of @Hodgkins.2018, taking into account previously unconsidered sources of uncertainties, shows that:

1. Larger aromatic contents may be necessary to stabilize tropical near-surface peat against decomposition and aromatics may accumulate at a higher rate with depth in tropical peatlands than estimated by @Hodgkins.2018.

2. Larger uncertainties indicate that also larger differences in aromatic contents between (sub)tropical and high latitude peat than estimated by @Hodgkins.2018 are possible, also between deeper peat.

3. More research should focus on how estimates of carbohydrate and aromatic contents from MIRS may be confounded by other organic matter fractions, in particular proteins. As a first step, this requires accurate concepts to name organic matter fractions and variables used in the interpretation of MIRS.

The results of this reanalysis also apply to a more recent extension of @Hodgkins.2018 with a larger dataset [@Verbeke.2022].


# Supporting information {.unnumbered}

Supporting information S1 is available as appendix to this manuscript.

# Acknowledgements {.unnumbered}

We like to thank Jeffrey Chanton and Suzanne Hodgkins for a constructive discussion which greatly helped to improve a first version of this manuscript. This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) grant no. KN 929/23-1 to Klaus-Holger Knorr and grant no. PE 1632/18-1 to Edzer Pebesma.

# Competing interests {.unnumbered} 

The authors declare no competing interests.

# Author contributions {.unnumbered}

HT performed the calculations, prepared the figures, and wrote the original text. Both HT and KHK revised and edited the text. KHK provided funding for this study.

# Data and Code availability {.unnumbered}

Data and code to reproduce our analyses are available via <https://doi.org/10.5281/zenodo.10276230> [@Teickner.2023a].


# References {.unnumbered}