dd-paper.Rmd

---
title: "Estimation of the Degree of Decomposition of Peat and Past Net Primary Production from Mid-Infrared Spectra"
date: "`r format(Sys.time(), '%d %B, %Y')`"
journal: "`r rticles::copernicus_journal_abbreviations(journal_name = 'SOIL')`"
author:
  - given_name: Henning
    surname: Teickner
    affiliation: "1, 2"
    email: henning.teickner@uni-muenster.de
    corresponding: true
  - given_name: Julien
    surname: Arsenault
    affiliation: "3"
    email:
  - given_name: Mariusz
    surname: Gałka
    affiliation: "4"
    email:      
  - given_name: Klaus-Holger
    surname: Knorr
    affiliation: "1"
    email: klaus-holger.knorr@uni-muenster.de
#
affiliation:
  - code: 1
    address: Ecohydrology & Biogeochemistry Group, Institute of Landscape Ecology, University of Münster, 48149, Germany
  - code: 2
    address: Spatiotemporal Modelling Lab, Institute for Geoinformatics, University of Münster, 48149, Germany
  - code: 3
    address: "Département des Sciences biologiques, Université du Québec à Montréal, Montréal, H2X 1Y4, Canada"
  - code: 4
    address: "University of Lodz, Faculty of Biology and Environmental Protection, Department of Biogeography, Paleoecology and Nature Conservation, Banacha 1/3, 90-237 Łodz, Poland"
bibliography: '`r targets::tar_read(dd_references_file)`'
running:
  title: Estimation of the Degree of Decomposition of Peat and Past Net Primary Production
  author: Teickner et al.
# This section is mandatory even if you declare that no competing interests are present.
competinginterests: |
  The authors declare no competing interests.
# See https://publications.copernicus.org/for_authors/licence_and_copyright.html, normally used for transferring the copyright, if needed. 
# Note: additional copyright statements for affiliated software or data need to be placed in the data availability section. 
copyrightstatement: |
#  The author's copyright for this publication is transferred to institution/company. 
### The following commands are for the statements about the availability of data sets and/or software code corresponding to the manuscript.
### It is strongly recommended to make use of these sections in case data sets and/or software code have been part of your research the article is based on.
### Note: unless stated otherwise, software and data affiliated with the manuscript are assumed to be published under the same licence as the article (currently Creative Commons 4.0)
availability:
  #code: |
  #  use this to add a statement when having only software code available
  #data: |
  #  use this to add a statement when having only data sets available
  codedata: |
    Peat data are available from the pmird database [@Teickner.2025c]. Litter data are available from the pmird database (undecomposed litter, data from @Reuter.2019 and @Reuter.2019a), and @Arsenault.2024b (all data except spectra, which are available from @Teickner.2025n). The code to reproduce this manuscript and the mid-infrared spectra for samples in @Arsenault.2024b are available from @Teickner.2025n. An R-package to estimate the mixing model for peat samples is available from @Teickner.2025i.
#  sample: |
#    use this section when having geoscientific samples available
# videosupplement: |
#  use this section when having video supplements available
authorcontribution: |
  HT: Conceptualization, methodology, software, validation, formal analysis, investigation, visualization, writing - original draft. KHK: supervision, funding acquisition. JA: Planned and performed the litterbag experiment in @Arsenault.2024a. MG: Performed the macrofossil analysis in @Diaconu.2020, @Galka.2022a, and @Galka.2022. All authors: writing - review & editing.
disclaimer: |
#  We like Copernicus.
acknowledgements: |
    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.
appendix:
language: "english" # set the manuscript language. The copernicus.cls supports: english, german, and french
output: 
  bookdown::pdf_book:
    base_format: rticles::copernicus_article
    keep_tex: true
    #latex_engine: xelatex
#fontsize: 12pt
#link-citations: true
header-includes:
  - \usepackage{float}
  - \usepackage{booktabs}
  - \usepackage{multirow}
  - \usepackage{xr} \externaldocument[si-]{dd-supporting-info}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, eval = TRUE, fig.align='center', fig.pos = 'H')
```

```{r packages, include=FALSE}
library(knitr)
library(kableExtra)
library(pmird)
library(targets)
library(magrittr)
library(bib2df)
library(RMariaDB)
```


```{r}
targets::tar_source()
tar_load(dd_supporting_info)
targets::tar_load(dd_model_info)
targets::tar_load(dd_data_model)
targets::tar_load(dd_data_pmird_peat_cores_yhat)
dd_data_pmird_peat_cores_yhat <- purrr::map(dd_data_pmird_peat_cores_yhat, readRDS_rvars)
tar_load(dd_citations_data)
tar_load(dd_stan_1_mcmc_settings)
tar_load(dd_stan_1_kfold_comparison)
tar_load(dd_stan_1_kfold_2_comparison)
tar_load(dd_data_model_evaluation_1)
tar_load(dd_stan_1_fit_mcse)
tar_load(dd_stan_1_fit_rhat)
tar_load(dd_data_pmird_mineral)
targets::tar_load(dd_data_pmird_peat_cores_irpeat)
# dd_reconstruction_initial_mass_rate_1 <- readRDS_rvars(tar_read(dd_reconstruction_initial_mass_rate_1))
tar_load(dd_plot_19_1)
tar_load(dd_data_bona2018_moss_npp)
tar_load(dd_data_bona2018_moss_npp_summary_1)
tar_load(dd_data_bengtsson2021_moss_npp_summary_1)
targets::tar_load(dd_simulation_1)
targets::tar_load(dd_simulation_2)
targets::tar_load(dd_simulation_3)
targets::tar_load(dd_simulation_3_species)
dd_stan_1_kfold_2_rmse <- readRDS_rvars(targets::tar_read(dd_stan_1_kfold_2_rmse))
targets::tar_load(dd_simulation_6)
```



```{r, results='asis'}
cat(dd_make_author_list(rmarkdown::metadata$authors))
cat("  \n")
cat("  \n")
cat(dd_make_affiliation_list(rmarkdown::metadata$affiliations))
```


\introduction[Introduction]

```{r}
f_ash <- function(c0, ct) {
  1 - c0/ct
}

f_ash_ct <- function(c0, gamma) {
  c0 / (1 - gamma)
}

d_ash <-
  tidyr::expand_grid(
    c0 = c(0.25, 1.03)/100,
    c0_assumed = c(0.25, 1.03)/100,
    gamma = 
      c(0, 0.1, 0.5, seq(0, 0.999, length.out = 20)) |>
      unique()
  ) |>
  dplyr::mutate(
    ct = f_ash_ct(c0 = c0, gamma = gamma),
    gamma_ash = f_ash(c0 = c0_assumed, ct = ct),
    gamma_ash =
      dplyr::case_when(
        gamma_ash < 0 ~ 0,
        TRUE ~ gamma_ash
      ),
    gamma_error = gamma - gamma_ash
  )
```

Peatland mass and carbon balance are the difference of net primary production (NPP) as mass input and of decomposition as mass output. Theory and modeling studies suggest a complex feedback between decomposition and net peat accumulation [@Bauer.2004; @Belyea.2006; @Morris.2011b; @Waddington.2015; @Mahdiyasa.2022]. As a direct effect, decomposition decreases net peat accumulation because mass is lost from the system, but as an indirect effect, this mass loss decreases peat height, porosity, and hydraulic conductivity, which can increase surface wetness and therefore may increase litter production and decrease decomposition over longer time periods. Under some conditions, more decomposition can therefore increase future net peat accumulation [e.g., @Quillet.2013]. Therefore, if we could accurately describe how decomposition is controlled by environmental conditions, we could improve both reconstructions of past and predictions of future peatland dynamics.  
When analyzing peat accumulation and decomposition processes, an important quantity is the degree of decomposition ($\gamma$), the fraction of initial mass that was lost due to decomposition. It is directly measured in litterbag experiments and allows to estimate decomposition rates and how they are controlled by environmental conditions. Peatland models use these estimates to simulate decomposition in peatlands for decades to millenia [e.g., @Frolking.2010]. However, litterbag experiments only cover short time periods which leads to imprecise estimates for parameters that take effect only over longer time periods, such as the slow-down of decomposition due to increasing recalcitrance of the material left behind [@Frolking.2001]. The degree of decomposition is also used to estimate mechanical and hydrological properties that control the ecohydrological feedback, such as the Young's modulus [@Mahdiyasa.2022; @Mahdiyasa.2023], bulk density [@Frolking.2010], and saturated hydraulic conductivity [@Frolking.2010; @Morris.2011b; @Morris.2015a], but these relations have only been tested against qualitative decomposition indicators [@Dykes.2008; @Morris.2015a], if at all. Estimates for $\gamma$ of peat are also useful to disentangle peat accumulation scenarios: an unfortunate consequence of the ecohydrological feedback is that both shallow and deep water table levels may lead to very similar peat heights and masses, two variables often used to test peatland models, but litter production, decomposition losses, the net mass balance, and the degree of decomposition differ between these contrasting scenarios [@Ramirez.2023]. For this reason, estimates of $\gamma$ for peat samples would be useful to distinguish between these scenarios and to obtain more accurate estimates for processes rates [@Ramirez.2023]. Finally, if $\gamma$ is known for a peat sample, one can directly estimate the initial mass of the sample as $m/(1 - \gamma)$, where $m$ is the final mass of the sample. The NPP can be computed from the initial mass if the peat core is dated and roots are not dominant in peat formation. Thus, $\gamma$ is a central quantity when estimating parameters and reconstructing process rates.  
It is currently not possible to measure $\gamma$ for peat samples. Direct measurements of $\gamma$ are currently only possible in litterbag experiments. For peat, several measurable properties (decomposition indicators) have been suggested to estimate differences in $\gamma$. Currently, only the ash residue method [e.g., @Leifeld.2011a; @Leifeld.2011] allows to estimate $\gamma$ on an absolute scale, but this method has large errors [@Kruger.2016] and indirectly presupposes that one already knows the initial mass of a sample (we abbreviate $\gamma$ estimated with the ash residue method as $\gamma_\text{ARM}$). All other decomposition indicators we are aware of have only been used to estimate qualitative differences in the degree of decomposition of peat [@Biester.2014; @Zaccone.2018]. Problems in the application of decomposition indicators are that the exact relation to $\gamma$ is unknown, that confounding variables are known to weaken or to reverse the relation to $\gamma$, or worse, potential confounders are unknown, and that the values of decomposition indicators may be highly variable already for undecomposed peat because of differences in litter chemistry [e.g., @Biester.2014; @Scheffer.2001; @Limpens.2003; @Mahdiyasa.2022; @Broder.2012; @Stuart.2004; @Teickner.2022f].<!-- Such error sources are difficult to take into account in peatland models and it would therefore be more useful to develop a decomposition indicator that estimates $\gamma$ on an absolute scale. To better understand how available decomposition indicators relate to $\gamma$ and how differences in litter chemistry confound this relation, measurements of $\gamma$ and decomposition indicators from litterbag experiments could be analyzed for different litter types.-->    
To address these problems with measuring $\gamma$, we suggest to define properties of an ideal decomposition indicator and based on this develop a decomposition indicator that has these properties. An ideal decomposition indicator is linearly related to $\gamma$, well tested against litterbag data, does not confound differences in litter chemistry of undecomposed litter with decomposition, and is robust against other potential confounding factors, in particular mineral contents and mixing different litter types.  
To develop such a decomposition indicator, we suggest to develop a model that predicts $\gamma$ measured in litterbag experiments from mid-infrared spectra (MIRS). Predictions by these models are the suggested decomposition indicator, $\gamma_\text{MIRS}$. Such a model may fulfill the requirements of an ideal decomposition indicator because MIRS from litterbag experiments contain information about the relative abundance and interactions of many molecular structures [@Stuart.2004] which change in abundance during decomposition and differ between plant species [e.g., @Arsenault.2024a; @Chapman.2001; @Cocozza.2003; @Tfaily.2014]. In contrast to other decomposition indicators that use simple formulas and no training data, such as C/N or humification indices, using a statistical model and undecomposed litter from various species as reference should avoid confounding differences in the chemistry of undecomposed litter with differences in chemistry due to decomposition. In contrast to $\gamma_\text{ARM}$, no knowledge of initial states is required because $\gamma$ is predicted only with information from the decomposed sample. Previous studies used MIRS to predict various peat properties [e.g., @Artz.2008; @Chapman.2001], but, to our knowledge, not $\gamma$.  
Peat samples can be much more complex than samples used in typical litterbag experiments because peat can be a mixture of multiple components: litter of different plant species and organs, possibly with different degree of decomposition and different changes of chemical components (e.g., cellulose, lignin, lipids) during decomposition, and minerals. Therefore, it will be important not only to evaluate $\gamma_\text{MIRS}$ with homogeneous litter data, but also to evaluate $\gamma_\text{MIRS}$ for mixtures of litter types similar to peat samples. To generate a diverse set of mixtures, one can add spectra of litterbag samples multiplied by scale factors. According to the Beer-Lambert law and with constant path length, these scale factors are the mass fractions of the components [@Stuart.2004] because MIRS intensities are proportional to the relative abundance of molecular structures and the relative abundance of a molecular structure is a weighted average of the relative abundances of the molecular structure of all components, where weights are the components' mass fractions. Also $\gamma$ can be calculated for these mixtures from $\gamma$ of the individual components and the initial litter masses, which allows us to test $\gamma_\text{MIRS}$ for peat samples.  
Theoretical considerations suggest that $\gamma_\text{MIRS}$ will underestimate $\gamma$ for mixtures of litter types: whereas both MIRS intensities and $\gamma$ are weighted averages of the components' values, the weights are not the same. For MIRS intensities, the weights are the mass fractions of the components in the sample, as mentioned above, whereas for $\gamma$, the weights are the mass fractions of the components in the undecomposed sample. This difference also propagates to $\gamma_\text{MIRS}$ predictions: 

\begin{equation}
\begin{aligned}
\gamma_\text{MIRS}(t) &=& f^{-1}\left(\beta_0 + \sum_{k=1}^K\sum_{w=1}^W \beta_w h\left(I_{k,w}(t)\frac{m_k(t)}{m(t)}\right)\right)\\
\gamma(t) &=& 1 - \frac{m(t)}{m(t_0)} = \sum_{k=1}^K \gamma_k(t) \frac{m_k(t_0)}{m(t_0)},
\end{aligned}
(\#eq:eq-dd-gamma-bias-1)
\end{equation}

where $t_0$ is the time when the sample was undecomposed, $t$ is some time point later than $t_0$, $k$ and $w$ are indices for the component and wavenumber, $I_{k,w}(t)$ is the MIRS intensity at wavenumber $w$ in component $k$, $m_k(t)$ and $m_k(t_0)$ are the mass and initial mass of component $k$, $m(t)$ and $m(t_0)$ are the mass and initial mass of the mixture, $f^{-1}$ is the inverse link function of the prediction model for $\gamma_\text{MIRS}$, where we assume that the model is linear on the link scale with intercept $\beta_0$ and coefficients $\beta_w$ for each wavenumber, and $h$ is a transformation of raw intensity values due to spectral preprocessing (for example standard normal variate (SNV) normalization and $z$-transformation of predictors). Therefore, more decomposed components have a smaller weight for $\gamma_\text{MIRS}$ than for $\gamma$, suggesting that $\gamma_\text{MIRS} < \gamma$, unless all components have the same $\gamma$ or the sample consists only of one component.  
If the bias is not negligible, alternative strategies to accurately estimate $\gamma$ also for peat samples to the direct measurement of $\gamma$ with $\gamma_\text{MIRS}$ need to be developed. When $\gamma_\text{MIRS}$ accurately estimates $\gamma$ for individual litter types and $h$ is a function such that $h\left(I_{k,w}(t) \frac{m_k(t)}{m(t)}\right) = h(I_{k,w}(t)) \frac{m_k(t)}{m(t)}$, then it should be possible to estimate $\gamma_\text{MIRS}$ for mixtures from $\gamma_k(t)$ for each component $k$ as:

\begin{equation}
\gamma_\text{MIRS} = f^{-1}\left(\beta_0 + \sum_{k=1}^K (f(\gamma_k(t)) - \beta_0)\frac{m_k(t)}{m(t)} \right)
(\#eq:eq-dd-gamma-bias-2)
\end{equation}

At first glance, it does not seem helpful to estimate $\gamma_\text{MIRS}(t)$ from $\gamma_k(t)$ to estimate $\gamma(t)$, but many peat samples consist of one or few dominant components which then also dominate the value of $\gamma(t)$. Using equation \@ref(eq:eq-dd-gamma-bias-2) in a mixing model [e.g., @Cotrufo.2023] then allows us to estimate $\gamma(t)$ from $\gamma_\text{MIRS}(t)$ and we even get estimates for $\gamma_k(t)$ for free. Such a mixing model works when we know the mass fraction of each component ($m_k(t) / m(t)$). To analyze existing data, we suggest to approximate $m_k(t) / m(t)$ from the volume fractions in macrofossil analysis, as is also assumed in current peatland models [e.g., @Frolking.2010; @Tuittila.2013; @Quillet.2015]. The same strategy also works when integrating $\gamma_\text{MIRS}$ into peatland models. We note that theory also suggests similar biases for other decomposition indicators (that measure a property of the decomposing organic material), were these indicators used to estimate $\gamma$. Therefore, the bias is a general limitation when we attempt to estimate $\gamma$.  
<!--To our knowledge, this bias has not been analyzed before and it will therefore be interesting to see whether the theoretical model proposed above can describe predictions of $\gamma_\text{MIRS}$ and how much they differ from $\gamma$. --><!--If the bias is not negligible, it is of course necessary to develop alternative strategies to accurately estimate $\gamma$ also for peat samples than directly estimating it with $\gamma_\text{MIRS}$. We note that theory also suggests similar biases for other decomposition indicators (that measure a property of the decomposing organic material), were these indicators used to estimate $\gamma$. Therefore, the bias is a general limitation when we attempt to estimate $\gamma$.-->
Here, we address the problems of estimating $\gamma$ of individual litter types and peat, and the question whether it is possible to estimate $\gamma$ of individual litter samples and of peat with $\gamma_\text{MIRS}$. In this context, our aims are: 

<!--To summarize, open questions are how available decomposition indicators relate to $\gamma$,  how much they are confounded by differences in litter chemistry, and whether $\gamma$ of litter and peat samples can be predicted from MIRS. To address these knowledge gaps, the aims of our study are:-->

1. To develop a novel decomposition indicator, $\gamma_\text{MIRS}$, which we define as $\gamma$, as measured in litterbag experiments, predicted from MIRS with a spectral prediction model.

2. To evaluate how well the developed models predict $\gamma$ for peat samples as mixtures of litter types from (1) different species, (2) with different degree of decomposition, and (3) with admixtures of silicate minerals.  

3. To develop a simple mixing model that allows to estimate $\gamma$ of peat samples from bulk measurements and macrofossil volume fractions. 

4. To illustrate and discuss how $\gamma_\text{MIRS}$ can be used to address open questions in peatland research (including the relation between saturated hydraulic conductivity and $\gamma$), to improve process models, to reconstruct past NPP, and to define site-specific references for $\gamma$ and NPP for peatland restoration. 

<!--3. To analyze how well other commonly used decomposition indicators relate to $\gamma$ as measured in litterbag experiments.-->

To this end, we collected available data from undecomposed litter and litterbag experiments for which MIRS and measurements of some commonly used decomposition indicators are available. This enabled us to compute beta regression models that predict $\gamma$ from MIRS. Several models are computed because MIRS can be preprocessed in different ways and this may lead to differences in prediction errors and robustness. We evaluated the models on individual litter types using cross-validation and residual plots versus N contents, a known confounder of predictions from MIRS models [e.g., @Teickner.2022f]. To evaluate applicability to peat samples, we apply the models to three sets of mixed spectra: mixtures of undecomposed material from different species, mixtures of material from different species and with different degree of decomposition, and mixtures of litter with a silicate-rich peat sample. Here, we test whether $\gamma_\text{MIRS}$ underestimates $\gamma$ for mixtures of litter types as predicted by the theoretical model in equations \@ref(eq:eq-dd-gamma-bias-1) and \@ref(eq:eq-dd-gamma-bias-2). Based on this analysis, we develop strategies to avoid this bias and show how $\gamma_\text{MIRS}$ can be used to estimate $\gamma$ and reconstruct past NPP for peat cores.   
<!--To evaluate the suitability of other decomposition indicators, we selected a subset of suggested decomposition indicators (C/N, H/C, O/C, $\delta^{13}$C, $\delta^{15}$N, HI$_{1630/1090}$, $\Delta$G$_\text{f}^0$, and NOSC) to $\gamma$ for the same litter samples used to develop the prediction models. Element ratios and isotope values have been discussed as decomposition indicators in recent reviews and previous studies [e.g., @Biester.2014; @Zaccone.2018; @Drollinger.2019]. HI$_{1630/1090}$ is the humification index (HI) computed from MIRS as the intensity at 1630 cm$^{-1}$ divided by the intensity at 1090 cm$^{-1}$ [@Broder.2012]. Other HI have also been defined, but they are often closely related to HI$_{1630/1090}$ and we therefore only analyze HI$_{1630/1090}$. $\Delta$G$_\text{f}^0$ is the standard Gibbs free energy of formation and has been suggested as decomposition indicator in @Worrall.2018. NOSC is the nominal oxidation state of carbon [@Masiello.2008] and is expected to decrease with more decomposition because of the preferential decomposition of O- and H-rich carbohydrates.  -->
By developing and evaluating models that predict $\gamma$ of litter from MIRS, we contribute to the development of a method that allows to accurately estimate $\gamma$ of peat, therefore to better understand decomposition processes, reconstruct past NPP, and improve process models. This will contribute to a better understanding of peat accumulation processes, the definition of ecological baselines, and prediction of future peatland dynamics. <!--By analyzing how other decomposition indicators relate to $\gamma$ in litterbag experiments, we improve understanding of the confounding factors and the comparability of decomposition indicators among peat samples.-->




# Methods

## Training data

For our analyses, we used the following litterbag data: *Sphagnum capillifolium* (capitulum and first 5 cm of the stem) and *Typha latifolia* (leaves) samples incubated in bog peat and pools (oxic and anoxic conditions) in Canada [@Arsenault.2024a; @Arsenault.2024b], and *Phragmites australis* (leaves and rhizomes) grown under different nutrient availability and incubated under anoxic conditions in peat with different nutrient availability [@Reuter.2019; @Reuter.2019a; @Reuter.2020]. In addition, we used undecomposed litter samples (some collected under natural conditions, some grown in laboratory experiments with different nutrient availability) of `r dd_data_model |> dplyr::filter(dd_dataset_label == "pmird_undecomposed_litter") |> dplyr::mutate(taxon_rank_value = paste0("*", taxon_rank_value, "*")) |> dplyr::pull(taxon_rank_value) |> unique() |> sort() |> knitr::combine_words()` [`r dd_citations_data$pmird_undecomposed_litter_citations |> dplyr::mutate(BIBTEXKEY = paste0("@", BIBTEXKEY)) |> dplyr::arrange(YEAR) |> dplyr::pull(BIBTEXKEY) |> unique() |> paste(collapse = "; ")`] from the pmird database [@Teickner.2025c]. For undecomposed litter samples, we assumed a degree of decomposition of 0%. In total, `r nrow(dd_data_model)` litter samples were available for model training. All MIRS used in this study are measured in transmission mode on Fourier-transform infrared spectrometers (FTIR) as milled samples mixed with potassium bromide into transparent pellets. MIRS were measured on different devices (`r dd_data_model |> dplyr::mutate(measurement_device_first_30char_duplciated = measurement_device |> stringr::str_sub(start = 1L, end = 30L) |> duplicated()) |> dplyr::filter(! measurement_device_first_30char_duplciated & ! is.na(measurement_device)) |> dplyr::pull(measurement_device) |> knitr::combine_words()`).    

## Peat data {#dd-methods-peat-data}

To evaluate the models and to illustrate the application of $\gamma_\text{MIRS}$, we used bog peat core measurements derived from the pmird database: a core from the peatland Odersprungmoor (OD2) [@Galka.2022a], a core from the peatland Martinskapelle (MK1) [@Galka.2022], and a core from the peatland Mohoş (MH1) [@Diaconu.2020]. Macrofossil volume fractions from these cores were used to estimate the mass fractions of individual litter components in the peat samples and average water table depth reconstructed from testate amoebae (TA-WTD) were used to relate estimated $\gamma$ and NPP to changes in WTD. Age-depth models were estimated with $^{14}$C dates using rbacon [@Blaauw.2017] (supporting section \@ref(si-dd-sup-age-depth-model)). In addition, we also used data from three ombrotrophic Patagonian peat cores from @Broder.2012 and two bog peat cores with highly degraded drainage layers from the Venner Moor to evaluate prediction domains of the models (see subsection \@ref(dd-methods-prediction-domain-coverage)).

<!--All cores have MIRS, all cores except from @Broder.2012 have elemental contents, $\delta^{13}$C and $\delta^{15}$N values, cores from the Venner Moor have bulk density measurements, and cores from @Galka.2022a, @Galka.2022, and @Diaconu.2020 have macrofossils [---todo: update these information once the pmird database is updated].-->

## Prediction model development

We used beta regression models to predict $\gamma$ from intensities in preprocessed MIRS. To harmonize the spectra, we interpolated them to unit wavenumber resolution and clipped them to the range 600 to 4000 cm$^{-1}$. Depending on the variable and properties of the MIRS, different preprocessing steps may maximize predictive performance. Therefore, we computed models with different preprocessing steps in addition to the harmonization: (1) Baseline correction, signal normal variate (SNV) (model `r dd_model_names(1, dd_model_info)`), (2) First derivative spectra, SNV (model `r dd_model_names(2, dd_model_info)`), (3) Second derivative spectra, SNV (model `r dd_model_names(3, dd_model_info)`). After that, spectra were binned (bin width 10 cm$^{-1}$) to reduce the number of predictors and reduce correlation between predictors, and all predictors were $z$-transformed across all spectra. All preprocessing steps were performed with the ir package [@Teickner.2022i]. Baseline correction was done with an automated convex hull procedure [@Beleites.2021]. Spectra of some undecomposed litter samples have artefacts from CO$_2$, as indicated by peaks around 2360 cm$^{-1}$ [@Wallace.1997a]. To avoid that these peaks have an effect on the normalization or models, we linearly interpolated intensities between 2300 and 2380 cm$^{-1}$ and excluded variables from this range from the models. CO$_2$ also produces peaks around 670 cm$^{-1}$ and 3580 to 3780 cm$^{-1}$, but we did not correct these because they overlap with peaks from organic matter and are smaller than the peaks between 2300 and 2380 cm$^{-1}$ [@Wallace.1997a].   
Beta regression models were computed with brms [@Burkner.2018], using a logit link function, assuming a constant shape parameter, using a normal prior for the intercept, and regularized horseshoe priors [@Piironen.2017a; @Piironen.2017b] for the slopes (for each predictor variable). The regularized horseshoe prior shrinks coefficients to zero except where they are strongly related to the response variable, conditional on other predictors. To reduce overfitting, we defined a large amount of shrinkage, by assuming that `r unique(dd_model_info$par_ratio) * nrow(dd_model_info$x[[1]]$spectra[[1]] )` of the `r dd_model_info$x[[1]]$spectra[[1]] |> nrow()` predictor variables have non-zero coefficients [@Piironen.2017b]. Bayesian inference was done using Markov Chain Monte Carlo (MCMC) sampling with Stan [@StanDevelopmentTeam.2021a], using `r dd_stan_1_mcmc_settings$warmup` warmup iterations and `r dd_stan_1_mcmc_settings$iter - dd_stan_1_mcmc_settings$warmup` sampling iterations. Across all models, the maximum Monte Carlo standard error [@Vehtari.2021] for the degree of decomposition predicted for the training data and all analyzed peat samples was `r dd_stan_1_fit_mcse |> dplyr::filter(variable == "dd_mcse_mean") |> dplyr::pull(max) |> max() |> magrittr::multiply_by(100) |> round(1)`% for the mean, `r dd_stan_1_fit_mcse |> dplyr::filter(variable == "dd_mcse_sd") |> dplyr::pull(max) |> max() |> magrittr::multiply_by(100) |> round(1)`% (for standard deviations), and `r dd_stan_1_fit_mcse |> dplyr::filter(variable %in% c("dd_mcse_lower", "dd_mcse_upper")) |> dplyr::pull(max) |> max() |> magrittr::multiply_by(100) |> round(1)`% (for lower and upper 95% prediction interval boundaries). No model had divergent transitions and the largest rank-normalized $\hat{R}$ for model parameters was `r dd_stan_1_fit_rhat |> dplyr::pull(rhat) |> max() |> round(2)`, indicating convergence of the chains [@Vehtari.2021].

## Model evaluation for litter samples

### Prediction errors for individual litter types

We estimated the predictive accuracy of the models with 10-fold cross-validation (CV). CV folds were defined in two ways to evaluate predictive performance for two scenarios. In the first scenario, we used stratified CV, where each observation is assigned to a categorical variable and observations are split into the 10 folds while approximately preserving relative category frequencies. Each observation was assigned to a category based on the study and species, or, for samples from @Reuter.2019, based on the study, species, and the site where the litter was grown (to balance occurrence of rhizome litter with different nutrient contents). This procedure approximates a test of the models with samples similar to the training data, as far as this is possible with our small and heterogeneous dataset.  
In the second scenario, we used grouped CV, where each observation is assigned to a categorical variable and observations for the same categorical variable all are assigned to the same fold. The categorical variable was the same as for the stratified CV. This estimates predictive performance for novel litter types that are not part of the training data, for example the only litterbag data with *Sphagnum* litter is from @Arsenault.2024a, and therefore in this scenario a model not trained on decomposed *Sphagnum* litter predicts the degree of decomposition for decomposed *Sphagnum* litter. The procedure also implies that some folds only contain undecomposed litter because there are fewer categories for litterbag data than CV-folds. The purpose of this second scenario is to test model robustness for novel litter types.   
To compare models, we used the expected log predictive density (ELPD) computed on observations held out during CV. Model evaluation was performed with the loo package [@Vehtari.2019]. Following rules of thumb [@Sivula.2022], we assumed models to have equivalent predictive performance (according to the capability of our evaluation) when the difference of their ELPD ($\Delta$ELPD) is smaller than 4, and otherwise when $\Delta$ELPD is larger than two times its standard error (using normal approximation for $\Delta$ELPD). To give an easier to interpret performance metric, we also computed the root mean square error (RMSE).  

### Prediction domain coverage {#dd-methods-prediction-domain-coverage}

We analyzed what fraction of the analyzed peat samples is within the prediction domain of the models. The prediction domain of a model (based on @Wadoux.2021, defined in @Teickner.2023) is the range of the predictor variables covered by the training data. A sample is within the prediction domain if its preprocessed spectrum is within the range for each predictor variable, otherwise it is outside the prediction domain. If a sample is outside the prediction domain, it has spectral properties that cannot be interpolated from the training data and therefore the model extrapolates when making predictions for the spectra. In particular with models that overfit, this may lead to larger prediction errors than fitting errors for the training data.  

### Relations of model residuals to N content

We analyzed whether residuals are related to the nitrogen (N) content as proxy for proteins because it is known that protein content differs between species and increases during decomposition [@Reuter.2020; @Scheffer.2001; @Limpens.2003; @Biester.2014], and that proteins, despite their comparatively small amount in peat, are a main control of the height and shape of peaks that are also caused by aromatics [@Stuart.2004]. Thus, if one considers that the training data contains decomposed litter only from three species, the fit of the models may be confounded by proteins and this may bias predictions, especially under extrapolation. For some of the litter samples, no N measurements were available. For these samples, we predicted N contents from MIRS using prediction models from the irpeatmodels package [@Teickner.2025g; @Teickner.2023]. MIRS-predicted N contents are labelled N$_\text{MIRS}$ and we indicate whether predictions are within the training prediction domain of the model [@Teickner.2023].


## Model evaluation for mixtures of peat components

### Admixtures of minerals

The impacts of minerals on the estimated degree of decomposition was evaluated by adding a scaled version of a spectrum of a peat sample from the pmird database [@Drollinger.2019; @Drollinger.2020] with large mineral content (`r dd_data_pmird_mineral |> dplyr::mutate(res = (1 - loss_on_ignition)) |> dplyr::pull(res) |> round(2)` g g$^{-1}$, computed from loss on ignition; based on the mid-infrared spectrum (Fig. \@ref(si-fig:dd-plot-16)), many of them are silicates [@Stuart.2004; @Parikh.2014]) to selected spectra from the training data from different taxa and with different degree of decomposition. A range of scaling factors from absence of silicates to a clear dominance of silicates was chosen and we predicted the degree of decomposition for each spectrum created this way. For the same samples from the training data, we evaluated the impact of adding increasing amounts of a very decomposed sample (based on HI$_{1630/1090}$, peat structure, C content, and degree of decomposition predicted by all models).  

### Mixtures of undecomposed and decomposed litter

To evaluate how robust the models are for mixtures of litter from different taxa and with different degree of decomposition, we created pairwise mixtures from two sets of spectra. The first set consisted of randomly selected undecomposed litter samples from 15 taxa and the two most decomposed samples per taxon (where available). The spectra in this set were scaled in the range `r range(dd_simulation_3$components$scale_factor) |> knitr::combine_words(and = " to ")`. The second set consisted of the two most decomposed samples and the least decomposed sample per taxon, which are not already in the first set. All pairwise mixtures of the spectra from the first and second set were created and $\gamma_\text{MIRS}$ predicted with all three models. Values of $\gamma$ were computed with equation \@ref(eq:eq-dd-gamma-bias-1). The value of $\gamma_k(t)$ is measured in the litterbag experiments, the value of $m_k(t_0)$ is computed by assuming that the scale factors correspond to masses (as outlined in the Introduction) as: $m_k(t_0) = m_k(t) / (1 - \gamma_k(t))$. Values of $m(t)$ are computed as sum of the scale factors of the components. We then used $\gamma(t)$ measured for each component and equation \@ref(eq:eq-dd-gamma-bias-2) to estimate $\gamma_\text{MIRS}$ to test whether this theoretical model can reproduce $\gamma_\text{MIRS}$ for the mixtures, irrespective of the litter type.  
To evaluate in more detail how robust the models are for mixtures of undecomposed litter, we selected one sample per taxon from the undecomposed litter samples and created all pairwise mixtures with equal amounts of of both components. It was tested whether predictions of the models for these mixtures match the expected value for $\gamma$ of 0 g g$^{-1}$.    
<!--### Prediction of $\gamma_\text{MIRS}$ from $\gamma$ in process models

Our analysis described in the previous subsection showed that mixing two or more components with different $\gamma$ tends to underestimate $\gamma_\text{MIRS}$ when applying our models. As outlined in the Discussion, this is not surprising if one recalls that MIRS signals are proportional to the relative abundance of molecular structures and therefore each component contributes to the MIRS signal relative to its mass fraction in the mixture. Therefore, if two components started with the same mass at a degree of decomposition of 0% and reach different degrees of decomposition until measurement, then they contribute equally to $\gamma$ of the mixture (equation \@ref(eq:eq-dd-gamma-bias-1)), but the more decomposed component contributes less to the MIRS signal because of its smaller mass fraction in the decomposed sample. This implies that bulk MIRS measurements cannot be used to estimate $\gamma$ from $\gamma_\text{MIRS}$. [---todo: maybe move to discussion]-->



<!--### Qualitative behavior for peat cores

Finally, we evaluated qualitative behavior of the models for two peat cores from the peatland Venner Moor [@Teickner.2025c] and three peat cores from Patagonian peatlands [@Broder.2012]. The Venner Moor is a former bog that was used for small-scale peat extraction, drained in 1895 for channel constructions and partly rewetted since the 1970's [@Wilkens.1955; @NaturschutzzentrumKreisCoesfeldeV.2017]. The first core was taken in the part of the peatland that was used for peat extraction and rewetted; here, predictions should indicate a small degree of decomposition for peat formed after rewetting and a sharp transition to the remaining peat that was decomposed during drainage. The second core was taken up to the mineral base in an area without peat extraction that was not rewetted and is overgrown by trees; here, predictions should indicate a small degree of decomposition for the top 5 cm where undecomposed leaves and *Sphagnum* grow, a large degree of decomposition for the part above the water table level, and a smaller degree of decomposition for the, based on C/N, HI$_{1630/1090}$, and photographs, undecomposed water saturated lower part (except for the bottom sections). The Patagonian cores feature ash layers and therefore allow to explore how the models extrapolate for peat samples with larger mineral contents. We emphasize that for peat samples, $\gamma_\text{MIRS}$ underestimates $\gamma$ and can be interpreted as the largest possible degree of decomposition of the least decomposed component and the smallest possible degree of decomposition of the most decomposed component. Therefore, qualitative evaluations of $\gamma_\text{MIRS}$ with the cores mentioned above is possible.-->

<!--For the qualitative evaluation against the peat cores mentioned in the Introduction, we plot the suggested decomposition indicator as depth profiles. Macrofossil abundances and TA-WTD are plotted alongside depth profiles where they are available [---todo: consider doing cross-correlation analysis with the TA-WTD?].-->

<!--## Relation of other decomposition indicators to the degree of decomposition

First, we analyzed the relation between other decomposition indicators (C/N, H/C, O/C, $\delta^{13}$C, $\delta^{15}$N, HI$_{1630/1090}$, $\Delta$G$_\text{f}^0$, and NOSC) and the degree of decomposition measured in the litterbag samples by computing Pearson correlations and creating scatterplots that differentiate between litter types (plant species and organ).  
Second, we compared the variation of HI$_{1630/1090}$ for the undecomposed litter samples with the variation in the degree of decomposition predicted by the models for the same samples and with the variation of HI$_{1630/1090}$ for all analyzed peat samples with a C content larger than 0.3 g g$^{-1}$, in order to analyze in more detail whether HI$_{1630/1090}$ confounds differences in the chemistry of undecomposed litter with decomposition, compared to the suggested decomposition indicator.-->  
<!--Third, we analyzed the relation between the other decomposition indicators (now also including bulk density) and the degree of decomposition predicted by the models for the peat samples by computing Pearson correlations and creating scatterplots that differentiate between peat cores.-->
<!--To test hypotheses regarding EAC and EDC, we predict EAC and EDC from MIRS for two sets of peat samples: peat samples with *Sphagnum* abundance $>80%$, peat samples with sedge abundance $>60%$, and peat samples with shrub abundance $>60%$ from cores with macrofossil information; and for the litterbag samples, and plot the EAC and EDC versus the suggested degree of decomposition (or the measured degree of decomposition for litter samples).-->
<!--We used models from the irpeat package to predict O/C, H/C, $\Delta$G$_\text{f}^0$, and NOSC, missing values for bulk density, C and N contents, and C/N, and for all samples humification index HI$_{1630/1090}$ (the ratio of the intensity at 1630 cm$^{-1}$ to the intensity at $1090$ cm$^{-1}$) as defined in @Broder.2012. In many cases, the litter samples are outside the prediction domains of the corresponding models from the irpeatmodels package and therefore may be biased [@Teickner.2023]. For $\delta^{13}$C and $\delta^{15}$N values, the prediction models have large prediction errors and we therefore omitted missing values in our analyses. For these reasons, we consider our analysis of the relation of the degree of decomposition to O/C, H/C, $\Delta$G$_\text{f}^0$, and NOSC as exploratory analyses. To make clear whether values are measured or predicted from MIRS, we label all MIRS predicted variables with a subscript "MIRS", for example C$_\text{MIRS}$, or use graphical means to differentiate measurements from predictions.-->


## Estimation of $\gamma$ for peat cores from bulk measurements of $\gamma_\text{MIRS}(t)$ and reconstruction of past NPP

We developed a simple Bayesian mixing model based on equation \@ref(eq:eq-dd-gamma-bias-2) that estimates $\gamma_k(t)$ for all components in a peat sample from bulk measurements of $\gamma_\text{MIRS}(t)$ (predicted by any of the three models) and the mass fractions of the components:

\begin{equation}
\begin{aligned}
\gamma_\text{MIRS} &\sim& \text{beta}(\mu_{\gamma_\text{MIRS}} \phi_{\gamma_\text{MIRS}}, (1 - \mu_{\gamma_\text{MIRS}}) \phi_{\gamma_\text{MIRS}})\\
\mu_{\gamma_\text{MIRS}} &\sim& \text{beta}(\mu \phi, (1 - \mu) \phi)\\
\mu &=& \text{logit}^{-1}\left(\beta_0 + \sum_k^K(\text{logit}(\gamma_k) - \beta_0)\frac{m_k}{m} \right)\\
\phi &\sim& \text{gamma}(\alpha_{\phi}, \beta_{\phi})\\
\gamma_k &\sim& \text{beta}(\alpha_{\gamma_k}, \beta_{\gamma_k})\\
\beta_0 &\sim& \text{normal}(\beta_{0_\mu}, \beta_{0_\sigma})\\
\end{aligned}
(\#eq:eq-dd-gamma-mixing-model-1)
\end{equation}

where the third line is the same as equation \@ref(eq:eq-dd-gamma-bias-2), the first line is a measurement error model to consider prediction errors in $\gamma_\text{MIRS}$ (here, we approximate errors of the prediction model with a beta distribution with scale parameter $\phi_{\gamma_\text{MIRS}}$ computed from the average prediction error of the respective model), $\mu$ is the estimated true value of $\gamma_\text{MIRS}$ (without measurement errors), $\phi$ is the scale parameter of the beta distribution (with a gamma prior with shape and rate parameters $\alpha_{\phi}, \beta_{\phi}$), and $\gamma_k$ of component $k$ is modeled with a beta distribution with shape and rate parameters $\alpha_{\gamma_k}$ and $\beta_{\gamma_k}$. As outlined in the Introduction, linking this model to peat core data can be done with bulk measurements for $\gamma_\text{MIRS}$ and by assuming that $m_k / m$ equals the macrofossil volume fractions of the litter components. For individual litter types, the mixing model will only yield correct estimates for $\gamma_k$ when they have a relatively large mass fraction because $\gamma_\text{MIRS}$ is not sensitive to $\gamma_k$ if component $k$ has only a small mass fraction. To consider this error source, we used a $\text{beta}(1,1)$ prior for $\gamma_k$, which gives uniform weights to all possible values. Unless the sample consists of many components, the inability to correctly estimate $\gamma_k$ for components with small mass fraction in the sample has only a small effect on $\gamma$ and an even smaller effect on the reconstructed NPP. We expect large errors not accounted for by the mixing model only when some initially dominant component is completely or nearly completely decomposed and therefore not detectable.  
We applied this model to the three cores, OD2, MK1, and MH1 described in section \@ref(dd-methods-peat-data) to estimate $\gamma_k$ of all $K$ litter components and the $\gamma$ of the peat samples. With these estimates for $\gamma_k$, bulk densities estimated from MIRS [@Teickner.2023; @Teickner.2025g], and age-depth models for each core, we could then estimate the initial mass of each layer, which divided by the time range of aboveground litter formation in each layer without roots equals the aboveground NPP. Reconstruction of belowground NPP is only possible with peatland models. Our analysis considers errors in estimates for $\gamma_\text{MIRS}$, bulk density, peat ages, and $\beta_0$. We also developed an R package that allows to estimate the mixing model for own data [@Teickner.2025i].  


<!--## Reconstruction of past *Sphagnum* NPP {#dd-methods-reonstruction-npp}

Besides testing peatland models, a potential application of a decomposition indicator that directly estimates the degree of decomposition is to reconstruct initial layer masses by dividing the mass of a peat sample by the fraction of initial mass remaining ($1 - \gamma(t)$). Since $\gamma_\text{MIRS}(t)$ underestimates $\gamma(t)$, this approach also underestimates the initial mass or NPP and we therefore call an approach that reconstructs initial layer masses with an uncorrected $\gamma_\text{MIRS}(t)$ a naive reconstruction of initial masses (or NPP in the absence of roots). Based on the results of our simulations described in the previous subsection, we suggest that under some conditions the naive approach can be used to approximately reconstruct past aboveground NPP of individual litter types (see the Discussion, subsection \@ref(dd-discussion-gamma-mirs-bias), and supporting section \@ref(si-dd-naive-reconstruction-applied)). For the dated peat cores (MH1, MK1, OD2), we used this approach to estimate the average NPP for some *Sphagnum* species and compare these estimates to measured *Sphagnum* NPP from previous studies. The details of this analysis are described in supporting section \@ref(si-dd-naive-reconstruction-applied).
-->

<!--NPP was reconstructed for layers that met the following conditions: the entire layer had a macrofossil volume fraction of *Sphagnum* $\ge ` dd_reconstruction_initial_mass_rate_1 |> dplyr::pull(macrofossil_volume_fraction_threshold) |> unique()`$%, model `r dd_model_names(2, dd_model_info)` predicts a median degree of decomposition less than -->

## Relation between saturated hydraulic conductivity and $\gamma$

Using the estimated $\gamma$ for the peat samples and saturated hydraulic conductivity ($K_\text{sat}$) estimated from MIRS [@Teickner.2025g; @Teickner.2023], we tested whether $K_\text{sat}$-$\gamma$ relations hypothesized in previous studies agree with $K_\text{sat}$ and $\gamma$ for our peat cores. The hypothesized relations are:

1. An assumed sigmoidal relation as used in the Holocene Peatland Model (corrected versions of equations (10) and (15) in @Frolking.2010; corrected versions are available from the supporting info to @Treat.2022, $K_\text{sat}$ in cm min$^{-1}$, according to @Morris.2015a):

    \begin{equation}
    K_\text{sat} = \frac{150}{70} - \frac{3}{70} \left( \rho_\text{min} +     \Delta\rho\left(1 - 0.5 \left(1 - \text{erf}\left( \frac{(1 - \gamma) - c_3}{\sqrt{2} c_4}\right) \right) \right) \right),
    (\#eq:eq-dd-ksat-gamma-frolking2010-1)
    \end{equation}
    
    where $\rho_\text{min}$, $\Delta\rho$, $c_3$, and $c_4$ are parameters defined in @Frolking.2010 and $\text{erf}$ is the error function.

2. An assumed exponential relation [@Morris.2011b] as used in the DigiBog model [@Morris.2012]:

    \begin{equation}
    K_\text{sat} = 0.001 \exp(8(1 - \gamma))
    (\#eq:eq-dd-ksat-gamma-morris2011b-1)
    \end{equation}

3. An empirical relation that was developed with $\gamma$ estimates derived from C/N ratios [@Morris.2015a]:

    \begin{equation}
    K_\text{sat} \approx 1.2 (1 - \gamma) + 2.6147~\text{depth} - 0.877~\text{hollow}  - 0.020~\text{center} - 2.87,
    (\#eq:eq-dd-ksat-gamma-morris2015a-1)
    \end{equation}

    where $\text{depth}$ is the mid depth of a layer, $\text{hollow}$ is a dummy variable which is 1 when the layer is from a hollow, and $\text{center}$ is a dummy variable which is 1 when the layer is from the center of the peatland (in our calculations, we use $\text{depth} \in \{0.3, 0.5\}$ m, $\text{hollow} = 0$, and $\text{center} = 1$). 

# Results

## Predictive performance

```{r}
dd_stan_1_kfold_2_comparison_df <- 
  dd_make_kfold_comparison_dataframe(
    kfold_comparison = dd_stan_1_kfold_2_comparison, 
    kfold_target_name = "dd_stan_1_kfold_2", 
    dd_model_info = dd_model_info
  ) |>
  dplyr::mutate(
    elpd_diff_is_significant = elpd_diff != max(elpd_diff) & elpd_diff + 2 * se_diff <= 0.0
  )
dd_stan_1_kfold_comparison_df <- 
  dd_make_kfold_comparison_dataframe(
    kfold_comparison = dd_stan_1_kfold_comparison, 
    kfold_target_name = "dd_stan_1_kfold", 
    dd_model_info = dd_model_info
  ) |>
  dplyr::mutate(
    elpd_diff_is_significant = elpd_diff != max(elpd_diff) & elpd_diff + 2 * se_diff <= 0.0
  )
```


All models have a similar fit to the training data (Fig. \@ref(fig:dd-plot-15)) and similar predictive performance as measured by their ELPD during CV (Tab. \@ref(tab:dd-tab-loo-comparison-1)). Model `r dd_stan_1_kfold_2_comparison_df |> dplyr::arrange(elpd_diff) |> dplyr::slice(length(elpd_diff)) |> dplyr::pull(model_name)` had the best average predictive accuracy (= maximum average ELPD) for the stratified CV and model `r dd_stan_1_kfold_comparison_df |> dplyr::arrange(elpd_diff) |> dplyr::slice(length(elpd_diff)) |> dplyr::pull(model_name)` for the grouped CV. Using a normal approximation for $\Delta$ELPD and a two-sided significance level of 5%, only model `r dd_stan_1_kfold_2_comparison_df |> dplyr::filter(elpd_diff_is_significant) |> dplyr::pull(model_name)` in the stratified CV `r if(any(dd_stan_1_kfold_comparison_df$elpd_diff_is_significant) || ! all(dd_stan_1_kfold_2_comparison_df$model_name[dd_stan_1_kfold_2_comparison_df$elpd_diff_is_significant] == dd_model_names(3, dd_model_info))) {stop("[---todo: check if(any(dd_stan_1_kfold_comparison_df$elpd_diff_is_significant ||...]")}` has a detectably worse predictive accuracy than the model with best average ELPD (Tab. \@ref(tab:dd-tab-loo-comparison-1)). Large errors of $\Delta$ELPD in the grouped CV indicate a large variability in predictive performance of the model for the different folds, as expected when testing the model with litter samples from species not included in the training data.  
For all models, the fit to the training data is better on average (smaller average RMSE) than for predictions in the two CV scenarios (Fig. \@ref(fig:dd-plot-15)), indicating that all models overfit. The two models using derivative spectra had, on average, the best fit to the training data, but worse or not much better predictive accuracy in the CV scenarios, indicating an, on average, larger overfitting risk than for model `r dd_model_names(1, dd_model_info)`, which does not use derivative spectra.  
Predictive accuracy of all models was worse for the grouped CV than for the stratified CV (Fig. \@ref(fig:dd-plot-15)), as expected, because in the stratified CV models were tested with samples from similar studies and species as the data they were trained with, whereas in the grouped CV some folds contained samples from species the model was not trained with. In the grouped CV, $\gamma$ is underestimated by all models for decomposed *T. latifolia* and *S. capillifolium* samples from @Arsenault.2024a. Conversely, $\gamma$ is overestimated for undecomposed litter samples, with prediction errors larger than 10% for model `r dd_model_names(1, dd_model_info)` for `r dd_data_model_evaluation_1 |> dplyr::filter(id_model == 1 & cv_type == "Grouped CV" & abs(y - yhat_mean) > 0.1 & y <= 0.01) |> dplyr::mutate(taxon_rank_value = paste0("*", stringr::str_replace(taxon_rank_value, pattern = ".+ ", replacement = paste0(stringr::str_sub(taxon_rank_value, 1, 1), ". ")), "*")) |> dplyr::pull(taxon_rank_value) |> unique() |> knitr::combine_words()`, and for the other two models for `r dd_data_model_evaluation_1 |> dplyr::filter(id_model != 1 & cv_type == "Grouped CV" & abs(y - yhat_mean) > 0.1 & y <= 0.01) |> dplyr::mutate(taxon_rank_value = paste0("*", stringr::str_replace(taxon_rank_value, pattern = ".+ ", replacement = paste0(stringr::str_sub(taxon_rank_value, 1, 1), ". ")), "*")) |> dplyr::pull(taxon_rank_value) |> unique() |> knitr::combine_words()`. Model `r dd_model_names(3, dd_model_info)` has, both in the grouped and stratified CV, the largest average prediction errors for undecomposed litter (Fig. \@ref(fig:dd-plot-15)).

(ref:dd-plot-15-caption) Fitted values (first column) or predictions for folds held out during CV for the two CV scenarios (second and third columns) versus the measured $\gamma$ for the training data, for each model (rows). For litterbag data, different colors represent different taxa, for undecomposed litter all samples are labelled as "Other" and have the same color. For *P. australis* rhizomes from @Reuter.2020, we additional differentiate between initial N contents (low, medium, and high N content). Error bars are 95% prediction intervals. For each case, we give the average RMSE and lower and upper 95% confidence intervals computed from MCMC draws.

```{r dd-plot-15, fig.cap='(ref:dd-plot-15-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_15))
```

(ref:tab-dd-tab-loo-comparison-1-caption) Difference of expected log predictive density of all models to the best model ($\Delta$ELPD) and corresponding standard error (SE($\Delta$ELPD)) for stratified or grouped cross-validation. In each case, the best model is at the top.

```{r dd-tab-loo-comparison-1, tab.cap='(ref:tab-dd-tab-loo-comparison-1-caption)', message=FALSE, warning=FALSE}
dplyr::bind_cols(
  dd_stan_1_kfold_2_comparison_df |>
    dplyr::select(-elpd_diff_is_significant),
  dd_stan_1_kfold_comparison_df |>
    dplyr::select(-elpd_diff_is_significant)
) |>
  dplyr::mutate(
    dplyr::across(
      dplyr::where(is.numeric),
      function(.x) round(.x, digits = 2L)
      )
  ) |>
  kableExtra::kable(
    format = "latex",
    booktabs = TRUE, 
    longtable = FALSE,
    escape = FALSE, 
    row.names = FALSE,
    col.names = rep(c("Model", "$\\Delta$ELPD", "SE($\\Delta$ELPD)"), 2L)
  ) |>
  kableExtra::add_header_above(c("Stratified CV" = 3, "Grouped CV" = 3))
```


## Prediction domain coverage

A comparison of the prediction domains of the models to the spectral range covered by the analyzed peat samples indicates that many of the peat samples have MIRS different to the training data. HI$_{1630/1090}$ as alternative decomposition indicator, C or C$_\text{MIRS}$, and photographs of the peat samples, suggest that samples outside prediction domains are either more decomposed or have large mineral contents (Fig. \@ref(si-fig:dd-plot-17) and \@ref(si-fig:dd-plot-18)), which is not surprising since there are comparatively few decomposed samples and no mineral-rich samples in the training data. Data coverage by prediction domains differs between the models: `r round(sum(dd_data_pmird_peat_cores_yhat[[1]]$degree_of_decomposition_in_pd) / nrow(dd_data_pmird_peat_cores_yhat[[1]]) * 100, digits = 1)`% of the peat samples are within the prediction domain for model `r dd_model_names(1, dd_model_info)` and `r round(sum(dd_data_pmird_peat_cores_yhat[[2]]$degree_of_decomposition_in_pd) / nrow(dd_data_pmird_peat_cores_yhat[[2]]) * 100, digits = 1)`% and `r round(sum(dd_data_pmird_peat_cores_yhat[[3]]$degree_of_decomposition_in_pd) / nrow(dd_data_pmird_peat_cores_yhat[[3]]) * 100, digits = 1)`% for the other two models, which indicates that derivative spectra differ more between the peat and the training data data than underived spectra. Thus, the models extrapolate for more decomposed and mineral-rich peat and the prediction domain of model `r dd_model_names(1, dd_model_info)` covers most of the peat samples. 

## Relations of model residuals to N contents

<!--We analyzed whether residuals are related to N$_\text{MIRS}$ as proxy for proteins because it is known that protein content differs between species and increases during decomposition, and that proteins, despite their comparatively small amount in peat, are a main control of the height and shape of peaks that are also caused by aromatics. Thus, if one considers that the training data contains decomposed litter only from three species, the fit of the models may be confounded by proteins and this may bias predictions, especially under extrapolation.--> A plot of residuals versus N for different levels of the $\gamma$ suggests that model `r dd_model_names(1, dd_model_info)` underestimates $\gamma$ for samples in the training data with a $\gamma$ between ca. 20 and 40% and N content smaller than ca. 0.02 g g$^{-1}$ (Fig. \@ref(fig:dd-plot-8)). For the models computed with derivative spectra, this bias is smaller (Fig. \@ref(fig:dd-plot-8)). This may indicate that model `r dd_model_names(1, dd_model_info)` makes less accurate predictions for the peat samples analyzed here (`r dd_data_pmird_peat_cores_irpeat |> dplyr::mutate(nitrogen_content_1 = sum(as.numeric(nitrogen_content_1) < 0.02) / length(nitrogen_content_1)) |> dplyr::slice(1) |> dplyr::pull(nitrogen_content_1) |> magrittr::multiply_by(100) |> round(0)`% have a N content smaller than 0.02 g g$^{-1}$) and in general [@Loisel.2014].<!--, which could explain why model ... predicts a smaller average degree of decomposition for many peat samples from the Venner Moor and the Patagonian peatlands (Fig. \@ref(fig:dd-plot-10-2)). -->

(ref:dd-plot-8-caption) Model residuals (measured $\gamma$ - $\gamma_\text{MIRS}$) versus N or N$_\text{MIRS}$ for different levels of $\gamma$ (columns) and the three models (rows). Points are averages and lines and shaded areas averages and 95% confidence intervals of regression models fitted to the average residuals. "Reliability" refers to how reliable the N or N$_\text{MIRS}$ values are. Values of N$_\text{MIRS}$ are less reliable when they are not in the prediction domain, which is the case only for some of the undecomposed litter samples.

```{r dd-plot-8, fig.cap='(ref:dd-plot-8-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_8))
```


## Predictions with additions of a silicate-rich sample

<!--We analyzed the robustness of the models under extrapolation for decomposed and mineral-rich samples. As we did not have access to independent validation data, we conducted qualitative sensitivity analyses by selecting some litter samples from the training data, adding multiples of spectra of a silicate-dominated or a very decomposed (based on HI$_{1630/1090}$, peat structure, C content, and degree of decomposition predicted by all models), and then predicting the degree of decomposition for these modified spectra.   -->
For silicates, prediction errors of all models increase the more silicates are mixed into the samples such that 95% prediction intervals cover nearly the whole range of possible values for $\gamma_\text{MIRS}$ (Fig. \@ref(si-fig:dd-plot-1-1) to \@ref(si-fig:dd-plot-1-3)). The behavior of predicted medians differs between the models. For model `r dd_model_names(1, dd_model_info)`, predictions overall increase for undecomposed samples and decrease for more decomposed samples, but the differences are never larger than ca. 30%. Model `r dd_model_names(3, dd_model_info)` behaves similar, but predictions always increase with more silicate influence and differences to the true value may exceed 50%. For model `r dd_model_names(2, dd_model_info)`, predictions decrease with more silicate influence for all but the undecomposed samples and compared to the other models more posterior probability is placed on smaller degrees of decomposition.  

## Predictions with additions of strongly decomposed samples 

For decomposed peat, the median predicted $\gamma_\text{MIRS}$ increases for all samples and models, as expected when adding a highly decomposed peat sample to less decomposed litter samples (Fig. \@ref(si-fig:dd-plot-5-1) to \@ref(si-fig:dd-plot-5-3)). Median predictions for the maximum addition of the decomposed peat are smallest for model `r dd_model_names(1, dd_model_info)` and largest for model `r dd_model_names(3, dd_model_info)`. As for the silicate experiment, also prediction errors of all models increase the more of the decomposed sample is mixed into the samples, but prediction intervals are narrower than for the silicate experiment and, especially for already decomposed litters, their width also differs between models: Model `r dd_model_names(1, dd_model_info)` estimates the largest prediction errors and prediction errors increase the more of the decomposed peat is added to a sample, even for decomposed litter. Model `r dd_model_names(2, dd_model_info)` and in particular model `r dd_model_names(3, dd_model_info)` have narrower prediction intervals for already decomposed litter, but in most of the cases predictions do not differ significantly from those of model `r dd_model_names(1, dd_model_info)`. Thus, model `r dd_model_names(1, dd_model_info)` estimates a smaller median $\gamma_\text{MIRS}$ with larger prediction errors than the models using derivative spectra, but all models make extrapolations that are qualitatively reasonable.   


## Predictions for mixtures of different litter components

Predictions of all models for samples consisting of two litter types with different degree of decomposition underestimate $\gamma$. As shown in Fig. \@ref(fig:dd-plot-27), the bias is well described by our theoretical expectation (equation \@ref(eq:eq-dd-gamma-bias-2)). Thus, the magnitude of the bias is related to the difference in $\gamma$ of the mixed litter types and their relative mass fractions in the mixture. This relation is non-linear; when one of the two litter types is clearly dominant (mass ratio of 1:1000), the bias vanishes, such that for mass ratios in-between these extremes, the bias reaches a maximum. This maximum value is larger the larger the difference in $\gamma$ between the two litter types is. The difference in $\gamma$ also controls the position of the maximum. The more decomposed one of the components is compared to the other component, the more is the maximum shifted towards a smaller mass fraction of the more decomposed component.  
There is some deviation from this pattern due to prediction errors in $\gamma_\text{MIRS}$ for individual litter samples. Discrepancies not related to model prediction errors for individual litter types only occur when undecomposed litter of one species is mixed with strongly decomposed litter of some different species, where the mass ratio of the more decomposed litter type is between one third and two thirds (upper right and lower left panels in Fig. \@ref(fig:dd-plot-27), where diverse litter types are mixed with *P. australis* and where *S. capillifolium* litter is mixed with *P. australis* litter). If only undecomposed litter from different species and organs is mixed, $\gamma_\text{MIRS}$ correctly estimates a small $\gamma$ (values range between `r c(dd_simulation_6$mixture$degree_of_decomposition_1, dd_simulation_6$mixture$degree_of_decomposition_2, dd_simulation_6$mixture$degree_of_decomposition_3) |> mean() |> range() |> round(2) |> knitr::combine_words(" and ")` g g$^{-1}$ across all mixtures and models), which is also in line with our expectation. Overall the bias is similar across all models and for mixtures of different species and organs, and can be well described by equation \@ref(eq:eq-dd-gamma-bias-2), with smaller deviations caused by prediction errors for $\gamma_\text{MIRS}$ and variation between litter types. The bias estimated with model 2 matches the expected relation best.  

(ref:dd-plot-27-caption) Values of $\gamma$ minus $\gamma_\text{MIRS}$ (predicted by all three models) for mixtures of two litter types (components) versus the mass fraction of the first component. Columns indicate the second component of the mixture (indicated by $\gamma$ and species) and rows indicate the first component of the mixture (indicated by $\gamma$). Point colors indicate the species of the first component. Point shapes indicate the model with which $\gamma_\text{MIRS}$ was predicted. The black lines are calculated with equation \@ref(eq:eq-dd-gamma-bias-2) using average coefficients from model 1.    

```{r dd-plot-27, fig.cap='(ref:dd-plot-27-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_27))
```


<!--## Application to peat cores

Application of the models to peat cores from the drained and partly rewetted Venner Moor supports the sensitivity analysis results (Fig. \@ref(fig:dd-plot-10-2) (a)): Models `r dd_model_names(1, dd_model_info)` and `r dd_model_names(3, dd_model_info)` predict a large degree of decomposition for the decomposed (based on C/N, HI$_{1630/1090}$, peat structure) pre-drainage peat, and a small degree of decomposition for the peat formed since rewetting. Similarly, for the drained and forested core, both models correctly predict a small degree of decomposition for the top 5 cm, followed by more decomposed peat from the drainage acrotelm. Also, a smaller degree of decomposition in the mid part of the profile, which is still water saturated, and an increase in the degree of decomposition in the bottom part agree with C/N, HI$_{1630/1090}$, and peat structure. This analysis suggests that, even though many of the peat samples are not within the prediction domain of the models and we did not correct $\gamma_\text{MIRS}$ for multi-component systems, the suggested decomposition indicator currently can at least identify qualitative differences in the degree of decomposition.  
Application of the models to peat cores with volcanic ash layers from Patagonia also supports results from the sensitivity analysis (Fig. \@ref(fig:dd-plot-10-2) (b)): Whenever mineral contents are large (as indicated by smaller C$_\text{MIRS}$ compared to adjacent layers), prediction errors are large. In addition, for the ash layer at 170 cm depth in core PBR 2, model `r dd_model_names(1, dd_model_info)` and `r dd_model_names(3, dd_model_info)` predict a larger median degree of decomposition than model `r dd_model_names(2, dd_model_info)` which is consistent with results from the sensitivity analysis and indicates that model `r dd_model_names(2, dd_model_info)` may less reliably extrapolate for samples with large mineral contents. An interesting observation is that each model can predict different degree of decomposition values for samples with similar mineral content in core PBR 2 (Fig. \@ref(fig:dd-plot-10-2) (b)): For the ash layer at 170 cm depth, all models predict a larger median degree of decomposition whereas they predict a smaller median degree of decomposition for the bottom mineral layer. This may indicate that the suggested decomposition indicator or an improved version of it can predict the degree of decomposition for samples with moderate mineral contents, unlike humification indices (HI$_{1630/1090}$ also shown in the plots) that are also computed from MIRS, but suggest a smaller degree of decomposition in both cases.  
Since model `r dd_model_names(1, dd_model_info)` underestimates the degree of decomposition for training samples with intermediate degree of decomposition and N content smaller than 0.02 g g$^{-1}$, smaller estimates for $\gamma_\text{MIRS}$ by model `r dd_model_names(1, dd_model_info)` for cores from the Venner Moor and Patagonian cores are probably underestimates.



(ref:dd-plot-10-2-caption) Degree of decomposition ($\gamma_\text{MIRS}$) predicted by all three models, C content (measured or predicted from MIRS), HI$_{1630/1090}$ for two cores from the drained peatland Venner Moor (a) and cores from three Patagonian peatlands with volcanic ash layers (b). In (b), approximate position of ash layers from @Broder.2012 are shown as grey bands.

```{r dd-plot-10-2, fig.cap='(ref:dd-plot-10-2-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_10_2))
```
-->

<!--## Underestimation of the degree of decomposition in multi-component systems {#dd-methods-bias-simulations-results}-->

```{r}
dd_simulation_2_threshold_m_1_1_f_1 <- 
  dd_simulation_2 |>
  dplyr::filter(! is.na(bias)) |>
  dplyr::mutate(
        m_1_1_f = m_1_1 / m_1
  ) |>
  dplyr::group_split(gamma_1) |>
  purrr::map(function(.x) {
    .x |>
      dplyr::filter(abs(gamma_1 - gamma_2) == max(abs(gamma_1 - gamma_2)) & bias > -0.1) |>
      dplyr::filter(m_1_1_f > 0.05 & m_1_1_f < 0.95)
  }) |>
  dplyr::bind_rows() |>
  dplyr::pull(m_1_1_f) |> 
  min() |> 
  magrittr::multiply_by(100) |>
  round(0)

dd_simulation_2_threshold_bias_gamma_1 <- 
  dd_simulation_2 |>
    dplyr::filter(! is.na(bias) & abs(gamma_1 - gamma_2) == max(abs(gamma_1 - gamma_2))) |>
    dplyr::mutate(
        bias = abs(gamma_1 - gamma_mirs)
    ) |>
    dplyr::filter(bias <= 0.1) |>
    dplyr::pull(mf_1_1) |> 
  min() |>
  round(2)
```





<!--## Relation of other decomposition indicators to the degree of decomposition

For the litter samples used in litterbag experiments (*Sphagnum capillifolium*, *Typha latifolia*, and *Phragmites australis*), Pearson correlations between other decomposition indicators and the measured degree of decomposition have the expected sign (positive for  HI$_{1630/1090}$, $\delta^{15}$N, $\Delta$G$_\text{f}^0$; negative for C/N, H/C, O/C, and NOSC; mixed for $\delta^{13}$C depending on the relative importance of isotope fractionation versus enrichment of $^{13}$C-depleted compounds), but they confound differences in litter chemistry between species and plant organs with decomposition (Fig. \@ref(fig:dd-plot-12)). For undecomposed litter from all species, HI$_{1630/1090}$ varies in a larger range than for all peat samples with a C content larger than 0.3 g g$^{-1}$, whereas all models accurately reproduce that the samples are not decomposed (Fig. \@ref(fig:dd-plot-11)).

(ref:dd-plot-12-caption) Relation of decomposition indicators to the measured degree of decomposition for samples from litterbag experiments. For *P. australis* rhizomes, we additional differentiate between initial N contents (low, medium, and high N content) because of the contrasting response of high-N rhizomes. Parameter $\rho$ is the Pearson correlation computed with average values (not considering measurement or prediction errors). "Reliability" refers to how reliable the decomposition indicator values are. MIRS predicted variables are less reliable when they are not in the prediction domain, making some of our analyses more speculative than others. For *Phragmites australis*, $\delta^{13}$C and $\delta^{15}$N values are predicted from MIRS. All values for O/C, H/C, $\Delta$G$_\text{f}^0$, and NOSC are predicted from MIRS.--><!-- The last panel shows the degree of decomposition fitted by model `r dd_model_names(3, dd_model_info)` as comparison [---todo: better than showing fitted values is to show predictions for hold out folds during cross-validation. I should add that!].--> <!--Error bars are standard deviations, lines and shaded areas are ordinary least squares linear regression model averages and 95% confidence intervals. 

```{r dd-plot-12, fig.cap='(ref:dd-plot-12-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_12))
```

(ref:dd-plot-11-caption) HI$_{1630/1090}$, measured (or known) degree of decomposition, and degree of decomposition predicted by all models for undecomposed litter samples in stratified CV. The orange bar gives the range of HI$_{1630/1090}$ (times 100) for all peat samples analyzed in this study with a C content larger than 0.3 g g$^{-1}$. This allows to compare differences in HI$_{1630/1090}$ caused by differences in the chemistry of undecomposed litter to differences in HI$_{1630/1090}$ caused by differences due to decomposition. Taxa are arranged by their median N$_\text{MIRS}$ content (largest median N$_\text{MIRS}$ at the top) and points for HI$_{1630/1090}$ are scaled by their N$_\text{MIRS}$.

```{r dd-plot-11, fig.cap='(ref:dd-plot-11-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_11))
```
-->

<!--For peat samples, scatterplots and Pearson correlations agree with the analysis for litter samples, bulk density shows the expected positive correlation to the degree of decomposition and dependence on mineral contents (e.g. bottom layers of the core from the forested subsite of the Venner Moor), and the relations vary between cores and core sections, as expected from previous studies (Fig. \@ref(fig:dd-plot-6-1-2)). Also as expected, samples with larger mineral contents (C$_\text{MIRS}$ < 0.3 g$_\text{C}$ g$_\text{sample}^{-1}$) are outliers for HI$_{1630/1090}$ and decomposition indicators predicted from models not trained with mineral-rich samples (O/C, H/C, $\Delta$G$_\text{f}^0$). Samples with small estimated degree of decomposition in cores with macrofossil observations indicate, similarly to the litter samples, that the decomposition indicators confound differences in litter chemistry with decomposition (Fig. \@ref(fig:dd-plot-6-1-2)).

(ref:dd-plot-6-1-2-caption) Relation of decomposition indicators to the  degree of decomposition of peat samples predicted by model `r dd_model_names(2, dd_model_info)`. Columns are cores and rows are decomposition indicators. $\rho$ is the Pearson correlation computed with average values (not considering measurement or prediction errors) across all cores. "Reliability" refers to how reliable the decomposition indicator values are. MIRS predicted variables are less reliable when they are not in the prediction domain, making some of our analyses more speculative than others. All values for O/C, H/C, $\Delta$G$_\text{f}^0$, and NOSC are predicted from MIRS. For all other decomposition indicators, variable numbers of values were predicted from MIRS. Vertical error bars are standard deviations and horizontal error bars (for the predicted degree of decomposition are 50% prediction intervals).

```{r dd-plot-6-1-2, fig.cap='(ref:dd-plot-6-1-2-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_6_1)[[2]])
```

\clearpage-->

<!--Using all peat samples, we analyzed how the suggested decomposition indicator is related to various other decomposition indicators measured for the same samples or predicted from MIRS. Pearson correlations have the expected sign (positive for bulk density, HI$_{1630/1090}$, $\delta^{15}$N, $\Delta$G$_\text{f}^0$; negative for C/N, H/C, O/C, $\delta^{13}$C, and NOSC) [---todo: ref figure]. Scatter plots show that the relations vary between cores and core sections, as expected from previous studies. Also as expected, samples with larger mineral contents (C$_\text{MIRS}$ < 0.3 g$_\text{C}$ g$_\text{sample}^{-1}$) are outliers for HI$_{1630/1090}$ and decomposition indicators predicted from models not trained with mineral-rich samples (O/C, H/C, $\Delta$G$_\text{f}^0$, $\delta^{13}$C, $\delta^{15}$N).

Conducting the same analysis with the *Sphagnum capillifolium*, *Typha latifolia*, and *Phragmites australis* litter samples from @Arsenault.2024a and @Reuter.2020 while using the measured degree of decomposition produces patterns similar to the peat samples for the decomposed litter, and suggests that all of the analyzed decomposition indicators confound differences in litter chemistry between species and plant organs with decomposition. Similar patterns are visible for peat samples with macrofossil observations.-->  


<!--Assuming that the suggested decomposition indicator (model ...) accurately describes the degree of decomposition, the scatterplots, together with macrofossil data, suggest that many decomposition indicators confound decomposition with differences in litter chemistry between species, in line with suggestions of previous studies. 


For the two cores with macrofossil data, the decomposition indicator values of _undecomposed_ peat vary gradually between hummock *Sphagnum* peat at the one extreme and sedge or lawn and hollow *Sphagnum* peat at the other extreme (models ... and ...). Model ... predicts an on average larger degree of decomposition for many of these samples and therefore suggests a closer relation to the analyzed decomposition indicators except for C/N. We speculate that model ...  overestimates the degree of decomposition here, but currently our only clues for this hypothesis are that, as mentioned above, the prediction domain of model ... is more representative for undecomposed peat samples, and that we used a diverse set of undecomposed litter samples to compute the models.-->

<!--The relation of H/C and O/C to the suggested decomposition indicator are consistent with current understanding of peat decomposition: O/C and H/C decrease as decomposition starts because of the preferential decomposition of hemicellulose that has more O and H relative to C than most other compounds.-->

<!--## Reconstruction of past average *Sphagnum* NPP

```{r, eval=FALSE}
dd_reconstruction_initial_mass_rate_1_summary <- 
  dd_reconstruction_initial_mass_rate_1 |>
  dplyr::select(taxon_rank_value, core_label, id_model, degree_of_decomposition, initial_mass_rate, degree_of_decomposition_in_pd) |>
  dplyr::group_by(taxon_rank_value, core_label, id_model) |>
  dplyr::summarise(
    dplyr::across(
      dplyr::all_of(c("degree_of_decomposition", "initial_mass_rate")),
      posterior::rvar_mean
    ),
    degree_of_decomposition_in_pd = sum(degree_of_decomposition_in_pd)/length(degree_of_decomposition_in_pd),
    .groups = "drop"
  )
```


Table \@ref(tab:dd-tab-dd-reconstruction-initial-mass-rate-1) summarizes estimates of past *Sphagnum* NPP for each core and species or species group, using the degree of decomposition estimated by each of the three models, and Fig. \@ref(si-fig:dd-plot-20) shows estimates for individual layers versus time. Since all analyzed layers have a relative small degree of decomposition, according to all models, the estimates do not differ much between models when their errors are considered. However, since model `r dd_model_names(3, dd_model_info)` estimates a somewhat larger average degree of decomposition with larger prediction errors for more decomposed peat samples, and we divide by one minus the degree of decomposition, NPP estimates have larger errors for model `r dd_model_names(3, dd_model_info)` with long right tails and therefore upper 95% interval values up to ` dd_reconstruction_initial_mass_rate_1_summary |> dplyr::filter(id_model == 3) |> dplyr::mutate(initial_mass_rate = initial_mass_rate |> magrittr::divide_by(1000) |> posterior::quantile2(probs = 0.975)) |> dplyr::pull(initial_mass_rate) |> max() |> round(1)` kg m$^{-2}$ yr$^{-1}$.

(ref:tab-dd-tab-dd-reconstruction-initial-mass-rate-1-caption) Estimated average NPP (kg m$^{-2}$ yr$^{-1}$) for different *Sphagnum* species or species groups in layers of three peat cores and using the degree of decomposition predicted by the three models. Values are given as averages and boundaries of 95% prediction intervals and are averages of of the average NPP for individual layers for the same model, species or species group, and core. Superscripts are levels for the fraction of layers within the prediction domain for the prediction models: $^1$: $[0, 50]$%, $^2$: $(50, 80]$%, $^3$: $(80, 100]$%.

```{r dd-tab-dd-reconstruction-initial-mass-rate-1, tab.cap='(ref:tab-dd-tab-dd-reconstruction-initial-mass-rate-1-caption)', message=FALSE, warning=FALSE, eval=FALSE}
dd_reconstruction_initial_mass_rate_1_summary |>
  dplyr::mutate(
    degree_of_decomposition = degree_of_decomposition * 100,
    initial_mass_rate = initial_mass_rate / 1000,
    degree_of_decomposition_in_pd =
      dplyr::case_when(
        degree_of_decomposition_in_pd <= 0.5 ~ "1", 
        degree_of_decomposition_in_pd <= 0.8 ~ "2", 
        TRUE ~ "3"
      ),
    dplyr::across(
      dplyr::all_of(c("degree_of_decomposition", "initial_mass_rate")),
      function(.x) {
        digits =
          switch(
            dplyr::cur_column(),
            "degree_of_decomposition" = 0,
            "initial_mass_rate" = 1
          )
        paste0(.x |> mean() |> round(digits), " (", .x |> posterior::quantile2(probs = 0.025) |> round(digits), ", ", .x |> posterior::quantile2(probs = 0.975) |> round(digits), ")$^{", degree_of_decomposition_in_pd, "}$")
        }
    ),
    core_label =
      dplyr::case_when(
        core_label == "eb1005_MH1" ~ "MH1",
        core_label == "eb1006_MK1" ~ "MK1",
        core_label == "eb1018_OD2" ~ "OD2"
      ) |>
      factor(levels = c("MH1", "MK1", "OD2")),
    taxon_rank_value =
      paste0("$", taxon_rank_value, "$") |>
      stringr::str_replace_all(" ", "~")
  ) |>
  dplyr::arrange(taxon_rank_value, core_label, id_model) |>
  dplyr::select(-degree_of_decomposition, -degree_of_decomposition_in_pd) |>
  tidyr::pivot_wider(
    names_from = "id_model",
    values_from = "initial_mass_rate"
  ) |>
  kableExtra::kable(
    format = "latex",
    booktabs = TRUE, 
    longtable = FALSE,
    escape = FALSE, 
    row.names = FALSE,
    col.names = c("Species group", "Core", 1:3)
  ) |>
  kableExtra::kable_styling(latex_options = c("scale_down")) |>
  kableExtra::add_header_above(c(" " = 2, "Model" = 3)) |>
  kableExtra::collapse_rows(1, latex_hline = "none", valign = "top") |>
  kableExtra::row_spec(0, align = "c")
```
-->

## Degree of decomposition and reconstructed NPP for peat cores

Figure \@ref(fig:dd-plot-28) shows $\gamma$ and NPP of individual litter types estimated with $\gamma_\text{MIRS}$ and the average $\gamma$ and total aboveground NPP for the three peat cores versus time. For MH1, $\gamma$ was smaller than 50% throughout the last 500 years; larger values values were only estimated for two short phases around 1700 and in the last decades, where *S. cuspidatum* and *E. vaginatum* dominated and TA-WTD indicate drier conditions. The reconstructed median NPP at MH1 is larger than at the two other sites, the median was stable during the last 500 years and larger than 0.1 kg m$^{-2}$ yr$^{-1}$. Recent and past drying events coincide with a decrease in NPP. For MK1, the oldest core, the median $\gamma$ is ca. 75% for the fen phase and peat formed during 500 to 1800 CE, but similarly small as for MH1 for the intermediate layers dominated by hummock *Sphagnum* species (model 3 predicts a larger $\gamma$ here than the other two models). The NPP was similarly large as in MH1 only in the fen phase and much smaller (median smaller than 0.15 kg m$^{-2}$ yr$^{-1}$) afterwards. No fen phase is covered by OD2 and $\gamma$ was small and more stable than in the other cores, with only a small peak roughly during the period of the Little Ice Age. The NPP was similar to MK1. In contrast to all other cores the NPP strongly increased during the last decades in OD2. In many cases, the apparent mass accumulation rate (AMAR) is nearly identical to the reconstructed NPP, but AMAR underestimates the median NPP and decomposition mass losses for more decomposed layers. Many decomposition peaks appear to be related to peaks in reconstructed WTD (drying events) (for example MH1: ca. 1700 CE and in recent decades, MK1: ca. 1700, 1500, 1000 BCE, OD2: ca. 2700, 1700, 700, 400 BCE), but there are also WTD peaks that do not coincide with an increase in $\gamma$ (also not when possible time lags between WTD and $\gamma$ are considered, for example: MK1: ca. 3000 BCE, OD2: ca. 700 CE).       
Common patterns between all models and cores are that dominance of *Eriohporum* sp. in peat often coincides with a large $\gamma$, that a small $\gamma$ is estimated for the topmost layers (as expected), and that a rather small $\gamma$ ($<50$%) is estimated for many *S. fuscum*, *S.rubellum*, and *S. magellanicum* layers even though they are several thousands of years old (and similarly, the overall $\gamma$ for these layers). The dominance in *Eriohporum* is often preceded by an increase in $\gamma$ of *Sphagnum* peat (e.g., MK1: ca. 1300 BCE, OD2: ca. 700 and 0 BCE) and $\gamma$ often only decreases after *Eriohporum* volume fractions are smaller again. This pattern is compatible with secondary decomposition and colonization by *Eriohporum* caused by drier conditions.  
While the models indicate overall similar trends in $\gamma$ and NPP, values of $\gamma$ and NPP vary comparatively much between the three models: model 1 estimates the smallest $\gamma$ and NPP, model 3 the largest values, and model 2 estimates intermediate values. Model 3 does, for example, not estimate a decrease in $\gamma$ for MK1 during ca. 5000 to 4500 BCE. However, the estimated maximum values of $\gamma$ per core are similar. Overall, model 3 suggests a smaller average difference in $\gamma$ between phases with more decomposition losses and phases with less decomposition losses than the other two models. Large errors in $\gamma$ are mainly due to extrapolation for more decomposed samples. Large errors in NPP are mainly due to dating errors (unlike many studies, we considered these errors here).  
Our estimates for *S. fuscum* NPP cover the range of NPP from different studies compiled by @Bona.2018 ($`r dd_data_bona2018_moss_npp_summary_1 |> dplyr::mutate(dplyr::across(dplyr::starts_with("npp"), function(.x) .x/1000), npp = paste0(round(npp_mean, 2), " \\pm ", round(npp_sd, 2))) |> dplyr::pull(npp)`$ kg$^{-1}$ yr$^{-1}$, average $\pm$ standard deviation) [`r dd_data_bona2018_moss_npp_summary_1 |> dplyr::mutate(reference_id = purrr::map(reference_id, function(.x) paste0("@", .x))) |> dplyr::pull(reference_id) |> unlist() |> unique() |> paste(collapse = "; ")`], and are contained in the wider NPP range measured in monospecific *S. fuscum* patches in @Bengtsson.2021 ($`r dd_data_bengtsson2021_moss_npp_summary_1 |> dplyr::mutate(dplyr::across(dplyr::starts_with("npp"), function(.x) .x/1000), npp = paste0(round(npp_mean, 2), " \\pm ", round(npp_sd, 2))) |> dplyr::pull(npp)`$ kg$^{-1}$ yr$^{-1}$) [@Bengtsson.2020] or other studies [@Wieder.2016; @Wieder.2019]. While it is not surprising that our estimates are within these broad ranges, it is reassuring to see that median estimates produced with the suggested approach do not obviously contradict expected ranges for NPP. 
<!--
5. Possible confounding factors for large $\gamma$ of S. fuscum in model 3: (1) model 3 overestimates $\gamma$. (2) The assumption that macrofossil abundances equal mass fractions is wrong and macrofossil abundances for S. fuscum are larger than its mass fractions in the samples; this would imply that in the mixing model for samples with large $\gamma$, S. fuscum needs to explain much more of the large $\gamma$ than warranted based on its mass fraction.
-->

(ref:dd-plot-28-caption) Reconstructions of $\gamma$ and the NPP from $\gamma_\text{MIRS}$ for the three peat cores analyzed here. Columns show results for peat cores and rows for the three prediction models for $\gamma_\text{MIRS}$. The x axis is the age of the upper boundary of the layers, measured since the coring date. Panel (a) shows $\gamma$ estimated for taxa with a volume fraction of at least 40% in at least 15 samples (colored points are medians and error bars 90% confidence intervals) and the median $\gamma$ of each layer (grey shaded areas and grey line; the grey shaded areas are 50% and 90% confidence intervals). Panel (b) shows the NPP (or, in presence of roots, the initial mass rate) for taxa with a volume fraction of at least 40% in at least 15 samples (colored points are medians and error bars 90% confidence intervals), the median NPP of each layer (grey shaded areas and grey line), and the median apparent mass accumulation rate (AMAR) (brown line). The y axis is clipped at 1 kg m$^{-2}$ yr$^{-1}$. Blue lines are average reconstructed WTD (normalized). The color of the small rugs along the x-axis ("Is in prediction domain?") indicate whether the MIRS for the peat samples are within the prediction domain of the models for $\gamma_\text{MIRS}$.   

```{r dd-plot-28, fig.cap='(ref:dd-plot-28-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_28))
```


## Relation between saturated hydraulic conductivity and $\gamma$

The estimated relation between MIRS-estimated saturated hydraulic conductivity ($K_\text{sat}$) and $\gamma$ for the three peat cores and with $\gamma$ estimated from $\gamma_\text{MIRS}$ predicted by the three models is shown in Fig. \@ref(fig:dd-plot-29), together with $K_\text{sat}$-$\gamma$ relations suggested in earlier studies. At least for the cores analyzed here, default parameterizations differ from estimated averages (albeit one has to consider that both $K_\text{sat}$ and $\gamma$ estimates have relative large errors): In contrast to the relation implemented in the Holocene Peatland Model (equation \@ref(eq:eq-dd-ksat-gamma-frolking2010-1)), the data suggest that $K_\text{sat}$ decreases continuously with $\gamma$ and our data do not cover ranges of $\gamma$ large enough to test the suggested step-like decrease for strongly decomposed peat. The relation suggested in @Morris.2011b and used in the DigiBog model overestimates $K_\text{sat}$ and estimates a too strong decrease of $K_\text{sat}$ with $\gamma$. The relation suggested in @Morris.2015a fits our estimates well, but the estimated depth dependence is too strong (for example, if a depth of 50 cm is assumed for all samples, $K_\text{sat}$ would be underestimated, as shown in Fig. \@ref(fig:dd-plot-29); in contrast our estimates suggest a less strong control of depth ($\approx$ total stress) on $K_\text{sat}$ for our cores).

(ref:dd-plot-29-caption) Relation between MIRS-predicted saturated hydraulic conductivity ($K_\text{sat}$) and $\gamma$ for the three peat cores analyzed here (columns) and with $\gamma_\text{MIRS}$ predicted by one of the three models (rows). Points are average estimates and error bars are 90% prediction intervals. Point colors indicate the depth of the peat layers. Point shapes indicate whether MIRS from which $K_\text{sat}$ were predicted are within the training prediction domain of the model ("yes") or not ("no"). Predictions outside the prediction domain may have larger errors than estimated. The lines are $K_\text{sat}$-$\gamma$ relations suggested in the literature (see the text for details). 

```{r dd-plot-29, fig.cap='(ref:dd-plot-29-caption)', out.width='100%'}
knitr::include_graphics(tar_read(dd_plot_29))
```

# Discussion

To address the problem of estimating $\gamma$ for peat samples, we developed and evaluated a novel decomposition indicator, $\gamma_\text{MIRS}$, and developed the mathematical tools necessary to estimate $\gamma$ with $\gamma_\text{MIRS}$ even for complex mixtures of several litter types with different $\gamma$. This enabled us to analyze in detail $\gamma$ and the past NPP for three peat cores based on bulk MIRS. Overall, this suggests that the approach developed here is a promising first step towards accurate estimation of peat $\gamma$ and reconstruction of past NPP. As is the case for all novel and complex prediction models, intensive testing is required to address remaining limitations and reduce prediction errors and some of this work remains to be done by future studies. The diverse palette of tests used here both shows that $\gamma$ can be estimated with $\gamma_\text{MIRS}$ with reasonably small errors for the various litter types considered here, and helps to identify limitations and possible future improvements.  
In the next subsections, we synthesize the mixing model approach to compensate the bias in $\gamma_\text{MIRS}$ for mixtures of litter types (subsection \@ref(dd-discussion-mixtures-1)), we synthesize the evaluation of the prediction models, and we identify limitations which can be addressed in future studies (subsection \@ref(dd-discussion-model-evaluation-1)). We also discuss how the approach developed here can be even further improved to estimate $\gamma$ and reconstruct the NPP, how $\gamma_\text{MIRS}$ can be linked to process models and what other research questions may be addressed and illustrate this with our analysis of the three peat cores (subsection \@ref(dd-discussion-implications-1)).    

<!--In order to improve peatland models and the interpretation of peat cores, we suggested properties of an ideal decomposition indicator and based on this developed models that predict the degree of decomposition ($\gamma$) in litterbag experiments from MIRS ($\gamma_\text{MIRS}$). In addition, we analyzed how other commonly used decomposition indicators (C/N, H/C, O/C, $\delta^{13}$C, $\delta^{15}$N, HI$_{1630/1090}$, $\Delta$G$_\text{f}^0$, and NOSC) relate to $\gamma$. This allows us to better illustrate limitations of existing decomposition indicators (subsection \@ref(dd-discussion-other-di-1)), in particular how they confound differences in litter chemistry with decomposition, than previous studies and to evaluate our prediction models and advantages and disadvantages of $\gamma_\text{MIRS}$ (subsections \@ref(dd-discussion-gamma-mirs-advantages-1) and \@ref(dd-discussion-gamma-mirs-model-evaluation-1)). We argue that $\gamma_\text{MIRS}$ has advantages as qualitative and also as absolute decomposition indicator. At the same time, improvements of the prediction models is desirable to better achieve these advantages and we outline opportunities to do so (subsection \@ref(dd-discussion-improvements)).  
We also identified a new source of error of many decomposition indicators, including $\gamma_\text{MIRS}$: the underestimation of $\gamma$ in multi-component systems, and analyzed in detail the magnitude of this bias depending on the degree of decomposition and relative proportions of the components. When the aim is to estimate the degree of decomposition on an absolute scale, this bias cannot be neglected, and it will also distort interpretation of qualitative decomposition indicators because the bias is not linear. In subsection \@ref(dd-discussion-gamma-mirs-bias), we present strategies how to avoid or correct this bias in different settings. Using one of these strategies, we estimated the degree of decomposition of *Sphagnum* species in layers up to several thousand years old from three peat cores and also reconstruct their NPP (subsection \@ref(dd-discussion-naive-reconstruction)). In addition, we discuss how a multi-component analysis of $\gamma_\text{MIRS}$ can be used to estimate litter-specific decomposition rates and the slow-down of decomposition due to increasing litter recalcitrance (subsection \@ref(dd-discussion-analysis-of-long-term-decomposition)), both factors that introduce large errors into peat accumulation simulated with process models [@Frolking.2001; @Quillet.2013; @Teickner.2025f] and that are assumed to not be measurable with litterbag experiments or peat accumulation data [@Clymo.1998; @Frolking.2001; @Teickner.2025]. Thus, $\gamma_\text{MIRS}$ can be useful to improve process models, and $\gamma_\text{MIRS}$ may be an alternative or supplement to process models for the reconstruction of past NPP from peat cores.-->  


<!--When estimating the degree of decomposition of a peat sample, the bias can be avoided by estimating the degree of decomposition for each component of a peat sample with $\gamma_\text{MIRS}$; when testing peatland models, the bias can simply be included in the model without losing information; with bulk peat measurements, the degree of decomposition can be estimated with small errors under some conditions that we identify based on our simulation analyses (subsection ... [---todo]). Based on this, we estimate the degree of decomposition of *Sphagnum* species in layers up to several thousand years old from three peat cores and also reconstruct their NPP (subsection ... [---todo]).-->  


<!--We will discuss the results of our analysis of this issue and give recommendations adding to those already given in the Introduction for using $\gamma_\text{MIRS}$ and how better approximations may be achieved with future research. Using this knowledge, we estimated the degree of decomposition of *Sphagnum* species and their NPP during several periods in the Holocene.

We have developed models that predict the degree of decomposition from MIRS ($\gamma_\text{MIRS}$) that were preprocessed in different ways and conducted quantitative and qualitative checks for predictive accuracy and robustness under extrapolation, considering differences in nitrogen contents, mineral contents, and mixtures of undecomposed and very decomposed samples. Here, we discuss to what extent these models fulfill properties of an ideal decomposition indicator, how their predictions compare among each other, how they can be improved, and what the advantages and disadvantages of $\gamma_\text{MIRS}$ are compared to other decomposition indicators.

To explore how well other commonly used decomposition indicators (C/N, H/C, O/C, $\delta^{13}$C, $\delta^{15}$N, HI$_{1630/1090}$, $\Delta$G$_\text{f}^0$, and NOSC) are related to the degree of decomposition, we used scatterplots and Pearson correlation of these decomposition indicators versus measured degree of decomposition for the litter data and explored how the relations vary between litter types. For HI$_{1630/1090}$, we also estimated its variability for undecomposed litter compared to the peat samples analyzed in this study. We discuss, what these results imply for the suitability of these peat properties as decomposition indicators. 

In the Introduction, we described the issue that decomposition indicators whose signal depends on the amount of organic matter remaining will underestimate the degree of decomposition of a peat sample even when the same decomposition indicator perfectly fits the degree of decomposition for individual components of the sample. $\gamma_\text{MIRS}$ is not immune to this issue, but there are strategies that either avoid this problem or allow to use $\gamma_\text{MIRS}$ to estimate the degree of decomposition of the sample or a dominant component. We will discuss the results of our analysis of this issue and give recommendations adding to those already given in the Introduction for using $\gamma_\text{MIRS}$ and how better approximations may be achieved with future research. Using this knowledge, we estimated the degree of decomposition of *Sphagnum* species and their NPP during several periods in the Holocene.-->


<!--To develop and evaluate a decomposition indicator that directly estimates the degree of decomposition of individual litter types, is testable against litterbag experiments, does not confound differences in the chemistry of undecomposed litter with decomposition, and is easy to measure, we computed models that predict the degree of decomposition from MIRS that were preprocessed in different ways and conducted quantitative and qualitative checks for predictive accuracy and robustness under extrapolation, considering differences in nitrogen contents, mineral contents, and mixtures of undecomposed and very decomposed samples. These analyses suggest that the models make qualitatively sensible predictions for peat, do not confound differences in the chemistry of undecomposed litter with decomposition, and are robust against mineral additions in the sense that their prediction errors increase. Our analyses also suggest differences between the models and identify limitations that we elaborate in the following sections. We expect that many of these limitations can be addressed with additional training data.  
 
To explore how well other commonly used decomposition indicators (bulk density, C/N, H/C, O/C, $\delta^{13}$C, $\delta^{15}$N, HI$_{1630/1090}$, $\Delta$G$_\text{f}^0$, and NOSC) are related to the degree of decomposition, we used scatterplots and Pearson correlation of these decomposition indicators versus measured or predicted degree of decomposition for the litter data and peat data and explored how they vary between litter types. For HI$_{1630/1090}$, we also estimated its variability for undecomposed litter compared to the peat samples analyzed in this study. We discuss, how reliable our findings are and what these results imply for the suitability of these peat properties as decomposition indicators.  
We warn against a naive reconstruction of past NPP from peat masses and the suggested decomposition indicator and discuss approaches that avoid these problems. An error analysis suggests, however, that the naive approach can be used, under special conditions, to reconstruct *Sphagnum* NPP from whole layer masses and the degree of decomposition, with absolute bias smaller than 20%, which allowed us to estimate the NPP of *Sphagnum* species or species groups for three peat cores for some periods during the Holocene. Finally, we provide a list of future investigations to improve the suggested decomposition indicator. -->


## Estimating $\gamma$ from $\gamma_\text{MIRS}$ for mixtures of litter types {#dd-discussion-mixtures-1}

As outlined in the Introduction, $\gamma_\text{MIRS}$ is biased when litter types with different degree of decomposition are mixed. Our results suggests that equation \@ref(eq:eq-dd-gamma-bias-2) reasonably well fits $\gamma_\text{MIRS}$ for mixtures of litter types; larger deviations only occur for litter mixtures that are probably rare under natural conditions (mixtures of one to two thirds of an undecomposed litter with a two to one thirds of a strongly decomposed litter). Moreover, at least for models 2, these deviations are small compared to prediction errors. Thus, since  equation \@ref(eq:eq-dd-gamma-bias-2) reasonably well fits $\gamma_\text{MIRS}$ for mixtures of litter types and since $\gamma_\text{MIRS}$ fits $\gamma$ of individual litter types well, the bias can be avoided or corrected.  
There are three possible approaches to avoid or correct this bias: The first approach is to use bulk $\gamma_\text{MIRS}$ measurements and information on the relative mass fractions of each litter type in the sample (for example from macrofossil analysis) to estimate $\gamma$ from $\gamma_\text{MIRS}$ with the mixing model (equation \@ref(eq:eq-dd-gamma-mixing-model-1)) developed here; this corrects the bias. The second approach is to measure MIRS for the separate litter types and weigh the litter types (component-specific analysis). In these cases, $\gamma_\text{MIRS}$ can be estimated for each litter type and, with reasonable prediction errors, equals $\gamma$. This can be used to estimate the initial mass of the litter types, $m_i(t_0)$, from $\gamma_\text{MIRS}$ with: $m_i(t_0) = m(t_0) / (1 - \gamma_\text{MIRS})$. Then, equation \@ref(eq:eq-dd-gamma-bias-2) can be used to compute $\gamma$ of the entire sample; this approach avoids the bias. The third option is a combination of the first and section approach, where a component-specific analysis is conducted only for some litter types.  
The second approach can be used to estimate $\gamma$ more accurately, both for the sample and individual litter types, than is possible with the first approach This is particularly relevant for litter types that have only a small mass fraction because the mixing model cannot estimate $\gamma$ very accurately for these litter types (for this reason we only analyzed litter types with a macrofossil abundance of at least 40% in the peat samples). Even though the sample $\gamma$ is less sensitive to the values of litter types with small mass fraction than $\gamma$ of these litter types, also estimates for the sample $\gamma$ will have smaller errors. The third option shares this advantage and may be an option when it is difficult to separate all litter types (for example when litter types are strongly decomposed).  
This points to a possible practical difficulty: it is not yet clear how to exactly define litter types and mass fractions such as to avoid biases. The mixing model and the component-specific analysis require litter types to be separated, but future studies will need to test which plant organs need to be separated to avoid biases. For example, it is quite clear that tree branches will decompose at a different rate than leaves and therefore both need to be treated as separate litter types to avoid biases. But does the same apply, for example, also to *Sphagnum* branches versus stems? Here, we assumed that macrofossil volume fractions are directly proportional to mass fractions, as is done in process models, but this assumption is not well tested. Peat layers may also consist of peat formed over a long time range that is strongly compressed and may therefore consist of not much decomposed parts and strongly decomposed parts. In such cases, it may not be sufficient to define litter types only via taxonomic information because litter from the same species may have a large range of $\gamma$ values in such layers. On the other hand, component-specific analysis can be used to determine whether $\gamma$  of two litter types is sufficiently similar throughout the decomposition process such that the bias in $\gamma_\text{MIRS}$ is only small and both may therefore be treated as same litter type.  
Overall, it is exciting to see that the possibility to estimate $\gamma$ raises these issues and offers opportunities to test them and therefore to better understand peat chemistry and degree of decomposition. 


## Which properties of an ideal decomposition indicator does $\gamma_\text{MIRS}$ fulfill and where are there still limitations? {#dd-discussion-model-evaluation-1}

The main question we want to address here is to what extent the prediction models for $\gamma_\text{MIRS}$ fulfill the properties of an ideal decomposition indicator suggested in the Introduction and compared to previous studies: 

1. **Linear relation to $\gamma$:** Our model evaluation suggests that all three models fit $\gamma$ of individual litter types with reasonable errors and no overall bias (Fig. \@ref(fig:dd-plot-15)). For litter similar to the training data, prediction errors are comparable to or slightly larger than errors between litterbag replicates [e.g., @Teickner.2025] and therefore currently acceptable for many applications. Predictions for novel litter types should be treated with caution since there might be biases (Fig. \@ref(fig:dd-plot-15)), however the dataset here already covers quite many different litter types and therefore we expect that this risk is small, except for woody plant organs which are currently underrepresented in the training data. Overall, $\gamma_\text{MIRS}$ is linearly related to $\gamma$ of individual litter types.

2. **Tested against litterbag data:** While there are gaps in our training data, we are not aware of other studies that test prediction of $\gamma$ with decomposition indicators with an as diverse set of litter samples as done here.  

3. **No confounding by differences in litter chemistry:** The models successfully fit $\gamma$ for many diverse litter types. Whilst decomposed litter samples are underrepresented in the training data, undecomposed samples from many different taxa and all relevant plant functional types (non-*Sphagnum* mosses, *Sphagnum* mosses, sedges, herbs, shrubs, and trees) have a good fit in the model and, in most cases, also in the stratified CV. Also mixtures of different litter types are well predicted when the bias in $\gamma_\text{MIRS}$ is compensated. Overall, this indicates that our models only have a small bias for different litter types.

4. **No confounding by minerals:** MIRS are strongly impacted by silicates. It will therefore be unlikely, at least with the simple linear prediction models developed here, to accurately predict $\gamma$ from MIRS. For this reason, it is encouraging that $\gamma_\text{MIRS}$ has large prediction errors when mineral contents increase; conversely to $\gamma_\text{MIRS}$, humification indices indicate a smaller $\gamma$ for samples with mineral admixtures and therefore is biased by silicates.

5. **No confounding when mixing litter types:** While $\gamma_\text{MIRS}$ is a biased estimator for $\gamma$ for litter mixtures, this bias can be corrected, as discussed in the previous subsection. We also emphasize that such a bias is not specific to $\gamma_\text{MIRS}$, but is a property of any decomposition indicator that measures a property where weights when computing averages for samples depend on $\gamma$, which is the case for nearly all decomposition indicators (all we are aware of except for $\gamma_\text{ARM}$).

Predicting $\gamma$ is a complex task and as any complex method, $\gamma_\text{MIRS}$ has limitations, many of which can be addressed in future studies:

1. **Variability between models:** Currently, it is difficult to recommend a best model because the models have trade-offs that need to be evaluated with additional data. Our recommendation is to consider the variability of model predictions while taking into account known biases (see the next point) until further tests have been performed. While confidence intervals of $\gamma$ estimated by the models overlap even for samples outside of the prediction domain, the large variability indicates that at least some of these predictions are biased. Overall, we assume that predictions of model 2 are most accurate because model 1 underestimates $\gamma$ for litter with small N content and model 3 probably overestimates $\gamma$ for more decomposed samples. 

2. **Limitations of individual models:** Model `r dd_model_names(1, dd_model_info)` fits the data worst and underestimates $\gamma$ most for *T. latifolia* and *S. capillifolium* samples with N content < 0.02 g g$^{-1}$ and therefore probably underestimates $\gamma$ for many peat samples and compared to the other two models. However, it has the best average predictive performance in the stratified CV, which we think best approximates possible applications of the models to peat samples, and it is also the model that confounds the least differences in initial litter chemistry with decomposition in the CV (Fig. \@ref(fig:dd-plot-15)).  
Model `r dd_model_names(3, dd_model_info)` overfits the data most, based on the stratified CV, and confounds the most differences in initial litter chemistry with decomposition in the CV (Fig. \@ref(fig:dd-plot-15)), but when trained on all data, it fits $\gamma$ best and does not confound difference in the chemistry of decomposed litter with decomposition, and it does not underestimate the degree of decomposition for *T. latifolia* and *S. capillifolium* samples with N content < 0.02 g g$^{-1}$. For more decomposed samples, model `r dd_model_names(3, dd_model_info)` may overestimate $\gamma$, as indicated by high $\gamma$ estimates in the analysis of the peat cores (Fig. \@ref(fig:dd-plot-28)) and underestimation of the theoretically expected bias for mixtures of undecomposed and very decomposed components (Fig. \@ref(fig:dd-plot-27)).  
In all aspects mentioned in the previous two points, model `r dd_model_names(2, dd_model_info)` is intermediate between the other two models and in addition has the best average predictive performance in the stratified CV (Fig. \@ref(fig:dd-plot-15)) and best compensates the bias for mixtures of litters (Fig. \@ref(fig:dd-plot-27)). However, even though 95% prediction intervals also cover the complete possible range for samples with large mineral contents, the median predicted degree of decomposition and bulk of the posterior distribution decrease the more silicates are within a sample, which means that $\gamma$ is underestimated for samples with large silicate contents.

3. **Prediction domain coverage:** It would be useful to include material from more litter types (in particular woody plant organs) and samples with larger $\gamma$ than is currently the case in order to reduce extrapolation errors of the models for decomposed peat samples.

4. **Confounding by carbonates:** We could not test whether $\gamma_\text{MIRS}$ is biased by admixtures of carbonates which can be relevant in fen peat. It is likely that carbonate-rich samples confound $\gamma_\text{MIRS}$ because peaks caused by carbonates interfere with peaks caused by aromatics and proteins [@Tatzber.2007].

5. **Modeling approach:** We used a modeling approach that is robust against overfitting with relatively small sample sizes, but there are more flexible modeling approaches [e.g., @ViscarraRossel.2024] that may have smaller prediction errors or may better fit more diverse litter types or even make more accurate predictions for samples with minerals or carbonates. These modeling approaches need larger and more homogeneous datasets to avoid overfitting. 

6. **Limitations of the mixing model:** We developed a simple mixing model that corrects the bias in $\gamma_\text{MIRS}$ and also estimates $\gamma$ of individual litter types. Beside the open questions mentioned in the previous subsection, we suggest that this model can be extended to incorporate additional information. For example, we did not incorporate prior information on differences in decomposition rates between litter types and this may lead to more accurate estimates for $\gamma$. A natural extension of the mixing model is to link it to a process model that constrains $\gamma$ of individual litter types.

Overall, this evaluation suggests that $\gamma_\text{MIRS}$ fulfills many properties of an ideal decomposition indicator, even if not applicable to all peat samples. Despite many opportunities for improvements, we are not aware of a decomposition indicator that has the potential to predict $\gamma$ as well as does $\gamma_\text{MIRS}$ and that can be as easily measured, even with small sample amounts, as is the case for $\gamma_\text{MIRS}$.


## Implications for peatland research and peatland restoration {#dd-discussion-implications-1}

Here, we discuss some problems in peatland research and restoration that may be addressed with $\gamma_\text{MIRS}$ and illustrate some of these points with our analysis of the three mountain bog cores.    

1. **Estimation of long-term decomposition rates and improvement of process models**: Our analysis of the peat cores suggests a quite large variation in $\gamma$ for hummock *Sphagnum* species, with several peaks across the entire depth, even if the peat is several thousands of years old. These differences were probably caused by differences in aerobic decomposition losses since catotelm decomposition rates are assumed to vary less with depth and over time and since at least some layers with larger hummock *Sphagnum* $\gamma$ coincide with or precede WTD peaks (Fig. \@ref(fig:dd-plot-28)) [@Morris.2015]. Preservation of such variations in peat may allow to better estimate long-term decomposition rates of plant functional types and reconstruct past environmental conditions with process models.   
It is interesting that a small $\gamma$ (smaller than $50$%, median values smaller than 25%) is suggested by all three models (even model 3 that may overestimate $\gamma$, Fig. \@ref(fig:dd-plot-28)) for hummock *Sphagnum* litter in some of the peat samples several hundreds or even thousands of years old and from all three peat cores. Some studies suggest much larger figures for $\gamma$ [@Clymo.1984; @Shao.2022a], whereas other studies suggest values for $\gamma$ similar to those found here [@Young.2017; @Ramirez.2023]. A litterbag synthesis suggests that under aerobic conditions, *S. fuscum* may lose 50% of its initial mass after 20 to 50 years [@Teickner.2025f]. Assuming a catotelm decomposition rate of 0.0004 yr$^{-1}$, this would lead to an additional mass loss of ca. 15% over 1000 years and thus a mass loss ca. 20% larger than the upper 90% confidence interval of $\gamma$ predicted for these samples. This may suggest either that litterbag experiments overestimate decomposition rates (at least under these conditions), that the time until incorporation into the catotelm was even shorter than 20 to 50 years, or that other factors than moisture limited decomposition during these periods (e.g., thermodynamic constraints in the catotelm [e.g., @Beer.2007; @Blodau.2011], colder temperatures, etc.). When linked to process models, $\gamma_\text{MIRS}$ may help to address these questions and therefore to better understand long-term decomposition losses.  
Other open problems that may be addressed with $\gamma_\text{MIRS}$ are: (1) A more precise quantification of differences in decomposition rates between litter types. If (aboveground) litter from the same layer decomposed under the same environmental conditions, current decomposition models [e.g., @Frolking.2001] assume that differences in $\gamma$ rates are only controlled by differences in the maximum possible decomposition rates of these litter types. This allows to estimate differences or ratios of litter type-specific decomposition rates from $\gamma_\text{MIRS}$ of individual litter types (supporting information \@ref(si-dd-rate-differences-1)). (2) It is commonly assumed that decomposition rates slow down the more of the initial mass has already been decomposed and this is an important control of long-term peat accumulation rates [e.g. @Frolking.2001; @Clymo.1998]. However, this slow-down is difficult to quantify [@Frolking.2001; @Clymo.1998] because the effect only manifests over long time periods not covered by litterbag experiments and because of difficulties to accurately estimate decomposition parameters from peat cores; measurements of $\gamma_\text{MIRS}$ from many peat cores and samples across large time gradients may help to address this problem (supporting information \@ref(si-dd-rate-differences-1)). (3) Similarly to the slow-down due to decreasing litter quality, it is assumed that anaerobic decomposition rates decrease with depth due to thermodynamic limitations [e.g., @Beer.2007; @Blodau.2011], which can have large effects on long-term peat accumulation [@Quillet.2013]. Measurements of $\gamma_\text{MIRS}$ may complement targeted litterbag experiments to better estimate this decrease in anaerobic decomposition rates [@Teickner.2025f; @Frolking.2010]. 

2. **Analysis of the WTD-decomposition feedback**: Saturated hydraulic conductivity ($K_\text{sat}$) is an important control of peatland WTD and depends on the pore size distribution which in turn depends on the stiffness of the organic matter [e.g., @Mahdiyasa.2022]. Decomposition decreases the stiffness and therefore facilitates pore collapse, implying that $K_\text{sat}$ and $\gamma$ are negatively related [e.g., @Morris.2011b]. Several studies hypothesized relations between saturated hydraulic conductivity ($K_\text{sat}$) and $\gamma$ [@Frolking.2010 with corrected formulas as included in the supplement to @Treat.2022; @Morris.2011b; @Morris.2015a; @Mahdiyasa.2022; @Mahdiyasa.2023], but these hypotheses could only be evaluated against not well tested proxies of $\gamma$, such as ratios of C/N values [@Morris.2015a], if at all. To test these hypotheses, we predicted $K_\text{sat}$ for the peat samples analyzed here from MIRS and plotted these predictions versus $\gamma$ together with the suggested $K_\text{sat}$-$\gamma$ relations (Fig. \@ref(fig:dd-plot-29)). The results indicate that, at least for the cores analyzed here, default parameterizations differ from estimated averages (albeit one has to consider that both $K_\text{sat}$ and $\gamma$ estimates have relatively large errors). Estimation of $\gamma$ from peat MIRS would also allow to test the relation between $\gamma$ and Young's modulus suggested in @Mahdiyasa.2022 and @Mahdiyasa.2023, and thus also the suggested relations to porosity, bulk density, and $K_\text{sat}$ suggested therein, if measurements for Young's modulus of peat would be available. Our analysis indicates that $\gamma_\text{MIRS}$ may be useful to estimate processes relevant to better understand the ecohydrological feedback in peatlands and to improve existing peatland models.

3. **Definition of reference states for restoration and long-term monitoring of peatland states**: Peatland restoration usually aims to restore the peat accumulation function by increasing NPP and reducing decomposition losses. Here, $\gamma_\text{MIRS}$ can be useful to define baselines for the NPP of individual species and the overall NPP, and what fraction of the initial mass should be transferred from the acrotelm to the catotelm. For example, based on our results, one may target a NPP of ca. 0.2 kg m$^{-1}$ yr$^{-1}$ at MH1 as historical reference state for restoration since this corresponds to the past average NPP over several centuries where the small $\gamma$ indicates that much of this sequestered C is stored in the catotelm. Our analysis also suggests that a drying trend during the last decades halved the NPP (Fig. \@ref(fig:dd-plot-28)). Here, $\gamma_\text{MIRS}$ allows to distinguish between a reduction in NPP versus an increase in decomposition losses. While a process model analysis is necessary to estimate net mass or carbon balances [@Young.2021], a decrease in NPP and increase in $\gamma$ of surface peat very likely indicate a decrease in the net mass or carbon balance since the acrotelm contributes most to decomposition losses. Considering that no long-term measurements of gas fluxes are required to obtain this information, but only standard analyses of peat cores, this may be a cost-efficient approach to estimate how much the current state of a peatland differs from more pristine states in the past and to evaluate restoration over longer time periods. Our analysis also illustrates the possibility to define site-specific reference states for NPP of target communities that were realized at the site over long time periods because $\gamma_\text{MIRS}$ allows to infer that NPP has probably always been smaller at MK1 and OD2 than at MH1 (Fig. \@ref(fig:dd-plot-28)). 
The rather large errors in the reconstruction here can be reduced not only by improving prediction of $\gamma_\text{MIRS}$, but also by more precise dating and by replacing MIRS-predicted bulk densities with bulk density measurements.  
In summary, this suggests that $\gamma_\text{MIRS}$ may be useful to better understand peat accumulation, improve process models, define baselines for peatland restoration, and monitor long-term restoration trajectories.      


\conclusions[Conclusions]

To address the question whether it is possible to estimate $\gamma$ of peat samples, we developed models that predict $\gamma$ measured in litterbag experiments from MIRS ($\gamma_\text{MIRS}$) and a mixing model that compensates biases in $\gamma_\text{MIRS}$ for mixtures of litter types and therefore allows to estimate $\gamma$ from $\gamma_\text{MIRS}$ with bulk MIRS measurements.  
Our model evaluation suggests that the models accurately fit $\gamma$ of several litter types, that the models do not confound differences in litter chemistry or additions of minerals with decomposition, and that our mixing model can correct underestimation of $\gamma$ by $\gamma_\text{MIRS}$ for litter mixtures with known mass fraction of each component, which allows to estimate $\gamma$ of dominant components, to estimate the average $\gamma$ of the sample, and to reconstruct the aboveground NPP (in absence of roots) of dominant litter types and the entire vegetation community at the sampling location.  
As is the case for all novel and complex prediction models, intensive testing is required to address remaining limitations and reduce prediction errors and some of this work remains to be done by future studies. Main limitations are that there is a relatively large variability in median predictions of $\gamma$ between prediction models using differently preprocessed MIRS due to biases for litter with small N contents (model 1), possible overfitting (model 2 and 3), extrapolation for more decomposed peat samples and therefore larger prediction errors (all models), large prediction errors for mineral-rich samples, missing woody plant organs in the training data, and possible limitations we could not test here (for example the influence of carbonates). Despite these limitations, $\gamma_\text{MIRS}$ fulfills many of the properties of an ideal decomposition indicator and is a promising first step towards accurate estimation of peat $\gamma$ and reconstruction of past NPP, in particular compared to other the ash residue method, which has large errors, and other decomposition indicators, which cannot estimate $\gamma$ quantitatively.  
Our analysis of three mountain bog peat cores illustrates that, if the limitations of the models are considered, $\gamma_\text{MIRS}$ may be useful to estimate $\gamma$ and reconstruct past aboveground NPP for individual litter types and peat samples. We discussed that these estimates of $\gamma$ may be useful to address a number of questions relevant to understand peatland processes and to plan and monitor peatland restoration (in particular when studies measure MIRS for individual litter types): Estimating long-term rates of peat decomposition, using peat cores as long-term litterbag experiments to understand environmental controls (thermodynamics limitations, decrease of litter quality), testing hypotheses about long-term decomposition processes (e.g., the time taken to incorporate peat into the catotelm), estimating the feedback between decomposition and peat hydraulic properties (e.g., the $K_\text{sat}-\gamma$ relation), and definition of past reference states for NPP and decomposition losses. Preliminary results indicate that *Sphagnum* peat can have a small $\gamma$ even if it is several hundreds to thousands years old, and that at, least for the cores analyzed here, suggested $K_\text{sat}-\gamma$ relations differ from the averages estimated here and vary between peat cores.  
Thus, $\gamma_\text{MIRS}$ may be useful to improve understanding of peat accumulation dynamics, process models, and may be used for a cost-efficient definition of peatland reference states and monitoring restoration trajectories.


<!--The suggested decomposition indicator --- the degree of decomposition predicted from mid-infrared spectra, $\gamma_\text{MIRS}$ --- has several advantages over existing decomposition indicators: It is a direct estimate for the degree of decomposition of litter types, only requires transmission mid-infrared measurements, and is independently testable against litterbag experiments. We illustrated that other decomposition indicators (C/N, humification indices, $\delta^{13}$C, $\delta^{15}$N, O/C, H/C, NOSC, and $\Delta$G$_\text{f}^0$) confound differences in the chemistry of undecomposed litter between species with decomposition mass losses and we hypothesize that many of them can vary for undecomposed litter as much as they vary due to decomposition. In contrast, $\gamma_\text{MIRS}$ is much less confounded by differences in the chemistry of undecomposed litter because the prediction model was trained on a diverse set of litter samples. Moreover, unlike humification indices, it does not confound large mineral contents with a small degree of decomposition or differences in N contents or differences in aromatic contents of undecomposed litter as decomposition mass losses. A main limitation of our implementation of the decomposition indicator is that we do not use a systematic sample collection, but use litterbag samples with *Sphagnum capillifolium*, *Typha latifolia*, and *Phragmites australis* available from previous studies. Our training data therefore lack a more complete representation of decomposed litter from other species (in particular sedges, shrubs, and hollow and lawn mosses), samples with high mineral contents, and more decomposed samples, and future studies can address these limitations. The model to estimate the suggested decomposition indicator is available from the irpeatmodels package [---todo:ref] and can be used with the irpeat package [@Teickner.2022d]. We expect that the suggested decomposition indicator, or an improved version of it, is a suitable replacement for currently used decomposition indicators and that it will be useful to reconstruct past NPP and distinguish between environmental conditions that lead to similar net peat accumulation but with different water table depths and matter fluxes in peatland models.-->


---
abstract: |
  The degree of decomposition of peat ($\gamma$) is useful to understand peatland degradation and peat accumulation, to reconstruct past net primary production (NPP), and to improve peatland models. None of the available decomposition indicators allows to estimate $\gamma$ with sufficient accuracy. We suggest prediction of $\gamma$ measured in litterbag experiments from mid-infrared spectra (MIRS) as a novel decomposition indicator, $\gamma_\text{MIRS}$, and compute prediction models for $\gamma_\text{MIRS}$ with available litterbag experiments and litter data from diverse species from the Peatland Mid-Infrared Database. For individual litter samples, the prediction models fit the data well, have reasonable prediction errors (average RMSE between `r dd_stan_1_kfold_2_rmse$rmse |> mean() |> range() |> round(2) |> knitr::combine_words(and = " and ")` g g$^{-1}$), and neither confound differences in litter chemistry nor differences in silicate contents with decomposition losses. We show that an underestimation of $\gamma$ by $\gamma_\text{MIRS}$ matches theoretical expectations; it can therefore be compensated, using plant macrofossil analysis data as a first approximation to mass fractions of peat components and a simple mixing model, or it can be avoided with component-specific measurements instead of bulk measurements. This allows to estimate $\gamma$ of peat samples and of dominant litter types and therefore also to reconstruct past NPP. To illustrate the approach, we analyze three cores from European mountain bogs and discuss how it can be used to improve process models and support restoration of peatlands. In particular, we test previously suggested relations between the saturated hydraulic conductivity and $\gamma$, illustrate how $\gamma$ measured on individual litter types may allow to use peat cores as natural litterbag experiments, and define reference states for $\gamma$ and NPP for the three analyzed peat cores. Improvements to reduce prediction errors of the approach require more diverse litterbag data, especially woody species and more decomposed litter. Further improvements can be achieved with measurements of MIRS on individual macrofossil types instead of bulk measurements, and an improved estimation of mass fractions of macrofossil types in peat samples instead of assuming that macrofossil abundances equal macrofossil masses. 
---