A comprehensive porewater isotope model for simulating benthic nitrogen cycling: description, application to lake sediments, and uncertainty analysis

The combination of various nitrogen (N) transformation pathways (mineralization, nitrification, denitrification, DNRA, anammox) modulates the fixed-N availability in aquatic systems, with important environmental consequences. Several models have been developed to investigate specific processes and estimate their rates, especially in benthic habitats, known hotspots for N-transformation reactions. Constraints on the N cycle are often based on the isotopic composition of N species, which integrates signals from various reactions. However, a comprehensive benthic N-isotope model, encompassing all canonical pathways in a stepwise manner, and including nitrous oxide, was still lacking. Here, we introduce a new diagenetic N-isotope model to analyse benthic N processes and their N-isotopic signatures, validated using field data from the porewaters of the oligotrophic Lake Lucerne (Switzerland). As parameters in such a complex model cannot all uniquely be identified from sparse data alone, we employed Bayesian inference to integrate prior parameter knowledge with data-derived information. For parameters where marginal posterior distributions considerably deviated from prior expectations, we performed sensitivity analyses to assess the robustness of these findings. Alongside developing the model, we established a methodology for its effective application in scientific analysis. For Lake Lucerne, the model accurately replicated observed porewater N-isotope and concentration patterns. We identified aerobic mineralization, denitrification, and nitrification as dominant processes, whereas anammox and DNRA played a less important role in surface sediments. Among the estimated N isotope effects, the value for nitrate reduction during denitrification was unexpectedly low (2.8 ± 1.1 ‰). We identified the spatial overlap of multiple reactions to be influential for this result.

1 Introduction

Nitrogen (N) is an essential element for all living organisms (Xu et al., 2022) and often limits primary production in aquatic systems (Kessler et al., 2014). In order to meet the global demand for fixed N (nitrate, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, and ammonium, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$), industrial fixation of atmospheric dinitrogen (N₂) through the Haber-Bosch process now exceeds biological N₂ fixation, with unforeseeable consequences regarding the ability of the environment to remove the excess fixed N, leaving the global N cycle imbalanced (Kessler et al., 2014). High fixed-N in aquatic systems has detrimental environmental consequences (Denk et al., 2017; Yuan et al., 2023), including eutrophication, ecosystem deterioration, and greenhouse gas emissions (e.g., nitrous oxide, N₂O). Thus, understanding the fate of fixed N in aquatic ecosystems and quantifying N fluxes are crucial for global budget estimates (Pätsch and Kühn, 2008).

In aquatic systems, benthic habitats are important hotspots in the transformation of large amounts of fixed N (Dale et al., 2019; Pätsch and Kühn, 2008; Xu et al., 2022), owing to sharp oxyclines and the co-occurrence of aerobic and anaerobic processes. The active N cycle in these sediments is driven by the flux of organic matter (OM) from the photic zone along with elevated concentrations of other electron donors (Ibánhez and Rocha, 2017; Wankel et al., 2015). Aerobic reactions, such as nitrification (stepwise ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ oxidation to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ via nitrite, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, with N₂O as by-product), are usually restricted to the top few millimetres in OM-rich sediments (e.g., in small lakes) or extend several centimetres deep in OM-poor sediments (e.g., in large oligotrophic lakes and the ocean) (Pätsch and Kühn, 2008; Wankel et al., 2015). The fate of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, produced via nitrification either locally in the sediments or in the water column, determines a system's capacity to function as an efficient N sink (Wankel et al., 2015). Denitrification, the stepwise reduction of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ to N₂ (via ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂O), has been identified as a key pathway for anaerobic N removal. Additionally, anammox, the anaerobic oxidation of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ to N₂ using ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, can contribute to N loss (Ibánhez and Rocha, 2017; Kampschreur et al., 2012; Wankel et al., 2015), especially in oligotrophic lake sediments (Crowe et al., 2017). In anammox, partial oxidation of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ generates ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ as a by-product to provide reducing equivalents for the fixation of inorganic carbon (C) (Brunner et al., 2013; Strous et al., 1999). Counteracting N removal by anammox and denitrification, the dissimilatory ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction to ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ (DNRA) contributes to N retention (Denk et al., 2017; Ibánhez and Rocha, 2017; Rooze and Meile, 2016). The balance between these N-transforming reactions is strongly influenced by environmental conditions, particularly the ratio of organic C to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and oxygen (O₂) availability. For instance, DNRA may be predominant under high $\mathrm{C}:{{\mathrm{NO}}_{\mathrm{3}}}^{-}$ ratios (Ibánhez and Rocha, 2017; Kraft et al., 2011; Wang et al., 2020). Oxygen is a central regulator in this context: it controls the coupling of nitrification with denitrification, anammox and DNRA, and modulates N₂O production and consumption, with peak N₂O yields typically occurring at the oxic-anoxic interface (Ni et al., 2011). The spatial overlap of aerobic and anaerobic N cycling processes at this transition zone in sediments often results in very low concentrations of metabolic intermediates (e.g., N₂O) in porewater, complicating their measurements in natural benthic environments. This is particularly true for the analysis of natural-abundance DIN isotopologues, which provide critical insights into N-cycling reactions and pathways. However, measuring these isotopologues, especially low-concentration intermediates in porewater, is technically challenging, if not impossible at present. To overcome these limitations, isotope modelling has become an essential tool for quantifying rapid N turnover at the oxic-anoxic interface, and for evaluating environmental controls on N dynamics and isotope signatures across diverse settings (Denk et al., 2017; Wankel et al., 2015).

Natural abundance stable isotope measurements provide insights into the N cycle, and the fluxes within its pathways, as microbial processes impart unique isotopic imprints on the involved N pools (Lehmann et al., 2003; Rooze and Meile, 2016; Wankel et al., 2015). In most microbial processes, the isotopically lighter molecules are preferentially consumed, yielding ¹⁵N-depleted products and ¹⁵N-enriched substrates (normal N-isotopic fractionation) (Kessler et al., 2014), with few exceptions, such as ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation, which occurs with an inverse N isotope fractionation (Casciotti, 2009; Martin et al., 2019). The isotopic composition of a given N pool is expressed in δ-notation, δ¹⁵N (‰ vs. std) = $\left[\right(R_{\text{sample}}/R_{\text{std}})-\mathrm{1}]\times \mathrm{1000}$, where R is the isotope ratio ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$, and the internationally recognized standard is atmospheric N₂ (Denk et al., 2017; Martin et al., 2019). The extent of the isotopic fractionation for a reaction is quantified using the isotope effect, ε, defined as $\mathit{\epsilon }\left(\mathrm{‰}\right)=[\mathrm{1}-{(}^{\mathrm{H}}k{/}^{\mathrm{L}}k)]\times \mathrm{1000}$, where ^Hk and ^Lk are the specific reaction rates for the isotopically heavy and light molecules, respectively (Sigman and Fripiat, 2019). For instance, δ¹⁵N-${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ analysis can help differentiate reductive and oxidative pathways of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ consumption, as they are characterised by a normal and an inverse kinetic isotope effect, respectively (Dale et al., 2019; Martin et al., 2019; Rooze and Meile, 2016). Despite considerable efforts to estimate isotope effects for most N-transformation processes (Denk et al., 2017), isotope effects estimated in batch cultures often differ from in situ measurements (Martin et al., 2019). To date, only limited efforts have been made to develop comprehensive benthic isotope models that integrate multiple N-transformation processes in a stepwise manner, and assess the expression of their isotope effects in the porewater of aquatic sediments, validated with observational data (Denk et al., 2017; Rooze and Meile, 2016).

Existing N-isotope models address specific aspects of the N cycle (Denk et al., 2017), such as denitrification (Kessler et al., 2014; Lehmann et al., 2003; Wankel et al., 2015), ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation and reduction (Buchwald et al., 2018) or N₂O dynamics (Ni et al., 2011; Wunderlin et al., 2012). As denitrification is the primary pathway for fixed-N loss in many aquatic systems, models integrating dual ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ isotopes (Lehmann et al., 2003; Wankel et al., 2015) have been used for example, to constrain its partitioning between water-column and benthic denitrification (Lehmann et al., 2005), as well as the contribution of regenerated ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ supporting denitrification (Lehmann et al., 2004). Rooze and Meile (2016) combined isotope data with a reaction-transport model to investigate the influence of hydrodynamics on fixed-N removal, highlighting enhanced coupling of nitrification-N₂ production by benthic infauna. Buchwald et al. (2018) used dual ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ isotope analyses and a reaction-diffusion model to demonstrate the tight coupling of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation near oxic-anoxic interfaces, emphasizing the central role of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ in N recycling. In contrast, most N₂O modelling efforts (primarily concentration-based models) to date have focused on engineered systems such as wastewater treatment, where they have been used to assess N₂O production pathways under variable conditions, and to minimize its emissions (Ni et al., 2011; Wunderlin et al., 2012). Challenges in measuring N₂O isotopologues in natural settings, especially in sediment porewaters, have limited the broader application of N₂O isotopic approaches and led to the exclusion of N₂O from benthic N-isotope modelling efforts so far. Nonetheless, given the key role of N₂O in the N cycle, and its sensitivity to redox conditions, there is a growing need for modelling frameworks that integrate multi-species N-isotope dynamics, even in the absence of direct measurements of N-cycle intermediate like ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂O to more accurately capture the interconnected nature of N transformations in natural systems.

With this study, we introduce a comprehensive 1-D diffusion-reaction model, encompassing all canonical N-transformation processes and most DIN isotopologues, to assess the role of distinct environmental factors (e.g., OM reactivity, bioturbation) in shaping porewater N dynamics and the N isotopic signatures the different N transformations (and combinations thereof) generate. Furthermore, by considering the stepwise nature of the N-cycling pathways, the model quantifies and isotopically characterizes key intermediates (i.e., N₂O, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$), which serve as substrates for subsequent reactions (Martin et al., 2019). Moreover, the model acts as a valuable research tool for analysing process couplings (e.g., DNRA-anammox interactions) (Dale et al., 2019; Hines et al., 2012), which are crucial for accurately estimating N removal and recycling, and can influence the apparent isotope effects of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$. Incorporating N₂O isotopologues as state variables enables the model to resolve the relative importance of N₂O producing mechanisms across small-scale benthic oxic-anoxic interfaces, and to quantify their contribution to sedimentary N₂O emissions.

The application of a comprehensive diagenetic N isotope model to measured porewater profiles of selected inorganic N compounds often results in parameter identifiability issues. Specifically, similar fits to the observed data might be achieved with comparable accuracy using different parameter sets, each yielding distinct transformation rates. To reduce the risk of drawing erroneous conclusions from such identifiability problems, we employed the following modelling strategies:

• Use of prior knowledge. Prior knowledge informed both the development of the model structure and the selection of parameter values. The model parameterization was adapted as deemed necessary to effectively integrate this prior knowledge. This approach aims to produce a plausible representation of the mechanisms governing the data.

• Consideration of uncertainty. Uncertainty in model parameters was explicitly accounted for using epistemic probability distributions. Bayesian inference (Bernardo and Smith, 1994; Gelman et al., 2013; Robert, 2007) was employed to combine prior knowledge with information obtained from observational data. The resulting posterior distribution of the parameters and calculated results provide a comprehensive uncertainty description, which is, however, still conditioned on prior information about the model structure and parameters.

• Sensitivity analysis. To test the robustness of key results against modelling assumptions, we assessed their sensitivity to the choice of prior probability distribution of the model parameters and to the inclusion of specific active processes within the model.

Since the numerical implementation of Bayesian inference requires the computationally intensive Markov Chain Monte Carlo (MCMC) sampling technique (Andrieu et al., 2003), an efficient model implementation is required. To meet this need, we implemented the model in Julia (Bezanson et al., 2017) (https://julialang.org, last access: 11 July 2024), a high-performance programming language. This choice also enables the use of automatic differentiation, which supports advanced MCMC techniques like Hamiltonian Monte Carlo (HMC) (Betancourt, 2017; Neal, 2011). The model was tested using field measurements from oligotrophic Lake Lucerne. It is important to emphasize that this isotope model is designed as a research tool, rather than a predictive instrument. Its primary purpose is to test hypotheses and assumptions related to the biogeochemical controls on N isotope signatures in natural environments, and to assess the identifiability of process rates and N isotope effects from observational data.

Figure 1 Simplified scheme of the N-transformation reactions considered for the diagenetic isotope model described in this paper. Continuous lines identify aerobic processes, while dashed lines indicate anaerobic processes. The state variables explicitly modelled as substrates for the considered reactions are highlighted with outlined boxes; O₂ is modelled as a state variable and as a regulator of aerobic and anaerobic processes; organic matter (OM) is not a state variable per se within the framework of this model, but acts as a source of N for the remaining processes. The isotopic fractionation of each process is shown using + and − signs to represent the ¹⁵N-enriching and ¹⁵N-depleting effects of the respective reactions.

Table 1 Chemical equations and reaction rate formulations for ¹⁴N and ¹⁴N¹⁴N compounds across all modelled processes. The rates for ¹⁵N, ¹⁵N¹⁴N, and ¹⁵N¹⁵N are formulated analogously by replacing the concentration of the isotopologue of interest as needed. The turnover rates for ¹⁵N-containing species are scaled by a factor of ($\mathrm{1}-\mathit{\epsilon }/\mathrm{1000}$), as outlined in the text. The complete set of equations including all isotopic compositions, and the process stoichiometry is provided in Appendix A: Model processes and stoichiometry. Anaerobic mineralization encompasses OM degradation coupled to iron and manganese reduction, as well as through methanogenesis.

2 Model description

2.1 Model formulation

A one-dimensional diffusion-reaction model was developed to simulate the concentrations of inorganic N compounds (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, N₂, N₂O), distinguishing between ¹⁴N and ¹⁵N isotopes (${}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${}^{\mathrm{14}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ${}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ¹⁴N₂, ¹⁴N¹⁵N, ¹⁵N₂, ¹⁴N₂O, ¹⁴N¹⁵NO, ¹⁵N₂O), as well as for O₂ and sulfate (${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$) concentrations. Their production and consumption rates are described by incorporating key processes of the canonical N cycle: aerobic mineralization, denitrification, nitrification, anammox, DNRA, mineralization by ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction, and anaerobic mineralization (other than ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ driven) (Fig. 1). All reactions (Table 1) are described using the general formula:

$$\begin{array}{}\text{(1)}& \text{rate}=k_{\text{max}}\cdot \text{limitation}\cdot \text{inhibition}\end{array}$$

where k_max represents the maximum conversion rate under ideal conditions (in µM d⁻¹). The terms for limitation by substrate X and inhibition by substance Y for the process i are defined following Michaelis–Menten kinetics (Martin et al., 2019):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }\text{limitation}={\displaystyle \frac{\left[X\right]}{K_{X,i}+\left[X\right]}}\text{(3)}& {\displaystyle }& {\displaystyle }\text{inhibition}={\displaystyle \frac{K_{Y,i}}{K_{Y,i}+\left[Y\right]}}\end{array}$$

where [X] and [Y] are the concentrations (in µM) of substances X and inhibitor Y, respectively, while K_X,i and K_Y,i are their respective half-saturation and inhibition constants (in µM) for process i, respectively. While the model supports exponential equations for limitation and inhibition terms, Michaelis–Menten kinetics were chosen for this study, as they are more commonly employed in N models (Rooze and Meile, 2016). The specific reaction rate equations are implemented taking into account the concentrations of ¹⁴N, ¹⁵N, ¹⁴N¹⁴N, ¹⁴N¹⁵N, and ¹⁵N¹⁵N species separately for the limitation term. For ¹⁵N-containing species, specific reaction rates are reduced by ($\mathrm{1}-\mathit{\epsilon }/\mathrm{1000}$) relative to ¹⁴N-containing species, reflecting the isotope effect associated with a given reaction (detailed descriptions of the model processes are provided in Appendix A: Model processes and stoichiometry).

Molecular diffusion is modelled taking into account the reduced solute movement due to tortuosity (Burdige, 2007). Additionally, bioturbation is included as a transport term enhancing diffusion, with its influence exponentially decreasing with depth. Boundary conditions are set based on observed concentrations of N compounds, O₂, ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ at the upper boundary, and by zero fluxes at the lower boundary, except for ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$. The ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ flux (and its ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{F}}_{{\mathrm{NH}}_{\mathrm{4}}}}$) was jointly estimated with the model parameters, as the field data display a clear ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentration gradient at 5 cm. Total N, ¹⁴N and ¹⁵N concentrations, along with their fluxes, are used for model parameterization (see Appendix B: Reaction-diffusion model for details).

The model is formulated as a dynamic model, but simulated to steady-state for comparison with observational data. Concentrations of ¹⁴N- and ¹⁵N-containing compounds are converted to total concentrations and δ¹⁵N.

2.2 Description of modelled transformation processes

This section outlines the modelled processes for ¹⁴N and ¹⁴N¹⁴N compounds (Table 1). A comprehensive overview of the transformation processes for all isotopologues, and stoichiometric relations is provided in Appendix A: Model processes and stoichiometry.

Mineralization of OM, the sole external N source, is differentiated in the model according to the specific electron acceptor involved: aerobic mineralization (O₂), denitrification and DNRA (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$), ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction, and anaerobic mineralization. The latter encompasses all remaining redox species (i.e., other than O₂, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$) below the nitracline (e.g., manganese oxides, iron oxides, carbon dioxide).

Denitrification is modelled as a three-step process: (1) ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$; (2) ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ to N₂O; and (3) N₂O to N₂. The first step, typically regarded as the rate-limiting step (Kampschreur et al., 2012), is the primary control on the overall expression of the N isotope effect (Kessler et al., 2014; Rooze and Meile, 2016). To prevent unrealistic rates, subsequent steps are constrained by setting $k_{\text{Den2}}=f_{\text{Den2}}\times k_{\text{Den1}}$ and $k_{\text{Den3}}=f_{\text{Den3}}\times k_{\text{Den1}}$, and specifying priors for f_Den2 and f_Den3. The re-parameterization of the second and third steps using the f_Den2Den1 and f_Den3Den1 factors corresponds to exactly the same model without any approximation or simplification. It serves solely to facilitate the specification of priors, as more knowledge is typically available about ratios of maximum rates (i.e., $f_{\text{Den2Den1}}=k_{\text{Den2}}/k_{\text{Den1}}$) than about the absolute maximum rates themselves. The ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ N isotope effect during benthic denitrification is known to be suppressed in the overlying water due to diffusion limitation (Dale et al., 2022; Kessler et al., 2014; Lehmann et al., 2003), though its expression at the porewater level remains less well constrained (Wankel et al., 2015). Transiently accumulating intermediates, such as N₂O, that can escape to the overlying water and alter benthic N fluxes (Rooze and Meile, 2016), are also considered. Lastly, to ensure mass balance, the model accounts for clumped (doubly substituted; e.g., ¹⁵N¹⁵NO and ¹⁵N¹⁵N) isotopocules, but does not distinguish between isotopomers (i.e., ¹⁴N¹⁵NO and ¹⁵N¹⁴NO) due to lack of N₂O isotope data needed for model validation. For the purpose of comparison with previous N models, a simplified one-step denitrification pathway (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ to N₂ with no release of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ or N₂O into the environment) approach is also implemented in the model code.

Nitrification is modelled as a two-step process: (1a) ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$; (1b) ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ to N₂O; (2) ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$. As for denitrification, the second step of nitrification is constrained to prevent unrealistic rates: $k_{\text{Nit2}}=f_{\text{Nit2}}\times k_{\text{Nit1}}$, with specifying a prior for f_Nit2. N₂O production yield during the first step is O₂-dependent, and is modelled accordingly:

$$\begin{array}{}\text{(4)}& f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\text{\_Nit1}}={\displaystyle \frac{ba}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+a}}\end{array}$$

where b and a are empirical parameters derived from Ji et al. (2018). N₂O production also occurs via nitrification-denitrification, implicitly modelled by allowing reaction coupling via the intermediate ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$. The expression of isotope effects depends on substrate availability and reaction completion. For instance, incomplete nitrification has been shown to result in isotopically heavy ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ efflux from the sediments (Dale et al., 2022; Lehmann et al., 2004; Rooze and Meile, 2016). However, similar phenomena for N₂O and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ remain poorly understood.

The limited understanding of porewater N isotope dynamics, especially for processes other than denitrification, hinges on the scarcity of isotope data for crucial N species like ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ in natural settings (Martin et al., 2019; Wankel et al., 2015). In the present model, we investigated the importance of these solutes, and how N-turnover processes like DNRA and anammox shape the distribution of their N isotopes. DNRA is modelled as a two-step process: (1) ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$; and (2) ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ to ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$. This approach separates the impact of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction on ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, and allows comparison of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ isotopic signatures induced by denitrification, DNRA, and anammox. Anammox is modelled to include both the comproportionation of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ to N₂ (main reaction, “m”), and the ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ production via ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation (side reaction, “s”) (0.3 $\mathrm{mol}{{\mathrm{NO}}_{\mathrm{3}}}^{-}$ produced per 1 $\mathrm{mol}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and 1.3 $\mathrm{mol}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$) (Tables 1 and A1) (Martin et al., 2019), which imparts a strong inverse isotope fractionation (Brunner et al., 2013; Magyar et al., 2021).

The relative importance of reductive ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ pathways is constrained by altering maximum conversion rates, k, as: $k_{\text{DNRA1}}=f_{\text{DNRA1,Den1}}\times k_{\text{Den1}}$; $k_{\text{DNRA2}}=f_{\text{DNRA2,Den2}}\times k_{\text{Den2}}$; $k_{\text{Anam}}=f_{\text{Anam,Den2}}\times k_{\text{Den2}}$, where prior information on f factors was obtained from experimental rate measurements (see below). Altogether these reactions provide a comprehensive overview of N isotope dynamics in porewater and enable the assessment of influential environmental conditions in shaping them.

2.3 Model assumptions

The model builds on the following considerations and assumptions:

• i.

  The inputs of sinking OM and associated advective transport relative to the sediment surface are not explicitly modelled, as the dissolved O₂ and N-compound profiles tend to reach quasi-steady state on short timescales (days to weeks). This simplification may not be valid for continental shelf sediments, where advection dominates solute movement due to high sediment permeability (Rooze and Meile, 2016). Therefore, in our model, porewater profiles are shaped primarily by molecular diffusion and bioturbation (the latter approximated as enhanced diffusion), along with reaction processes.

• ii.

  Hinging on assumption (i), the rates of OM-degrading processes are assumed to be limited by the availability of oxidants and not of OM, as in Kessler et al. (2014), an assumption that holds for sediments with sufficient readily degradable OM, but may break down at great depths. As OM is neither a state variable nor a limiting substrate, its production and consumption rates are not tracked and are considered uninfluential within the current model.

• iii.

  Microorganisms involved in N-transformation pathways are not explicitly modelled, meaning that maximum conversion rates, k, represent a combination of bacterial maximum specific growth rates and abundance. These parameters likely vary significantly across systems, due to differences in OM loading. Variabilities in cell-specific rates, and consequently in isotope effects, over depth and substrate availability were not considered.

• iv.

  N assimilation is not included, which is plausible if the turnover rates of the modelled processes are considerably higher than the N assimilation rates.

• v.

  Maximum specific conversion rates for all reactions are constant with depth, implying uniform bacterial abundance and activity across the sediment layer affected by any given process.

• vi.

  Limitation and inhibition kinetics are modelled using Michaelis–Menten functions, as they are commonly employed in N-cycle models (Rooze and Meile, 2016); exponential equations are provided within the code as an alternative approach, depending on user preference.

• vii.

  OM composition is approximated by the Redfield ratio ($\mathrm{C}:\mathrm{N}:\mathrm{P}$ = $\mathrm{106}:\mathrm{16}:\mathrm{1}$), used to estimate the fraction of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ released during OM mineralization, γ.

• viii.

  Anaerobic mineralization includes all processes involving redox species below the nitracline (e.g., manganese, iron, and carbon dioxide) with the exception of ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction, with no distinction in reaction rate for different oxidants. Reduction of ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ is modelled separately, as it can occur at faster rates than oxidation by iron(III), Fe³⁺, and manganese, Mn⁴⁺, in some lacustrine systems (Steinsberger et al., 2020), and is the dominant anaerobic mineralization process in marine settings.

• ix.

  Re-oxidation of reduced species other than ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (e.g., Fe²⁺, Mn²⁺, H₂S, CH₄) is neglected in the O₂ budget for the modelled interval; this is appropriate where their upward fluxes are minor, but may underestimate O₂ demand in settings with substantial reduced-species fluxes. Future users are encouraged to adapt the model to their research questions and dataset, including adding processes and state variables, provided that they can be constrained.

• x.

  OM mineralization occurs with no N isotopic fractionation; that is, the released ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ has the same N isotopic composition of OM, which is a model parameter considered for estimation.

• xi.

  Diffusivities of isotopologues are considered identical, as their differences have been reported to be minimal (Lehmann et al., 2007; Wankel et al., 2015).

• xii.

  Bioturbation enhances diffusion equally for all modelled species. As no solid was included as a state variable of the model, the impact of bioturbation on solid phase mixing was neglected.

• xiii.

  The yield of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ during anammox is fixed at 0.3 mol ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ per 1 $\mathrm{mol}{\mathrm{NH}}_{\mathrm{4}}^{+}$, although reported values range from 0.26 to 0.32 (Brunner et al., 2013).

• xiv.

  The ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ equilibrium during anammox has been previously reported to occur under environmental stress conditions with a strong isotopic fractionation (up to −60.5 ‰) (Brunner et al., 2013). Since it leads to the production of ¹⁵N-enriched ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, similarly to the kinetic isotopic fractionation during ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, variable values of ε_Anam,side (−15 ‰ to -45 ‰) can encompass both kinetic and equilibrium fractionation.

• xv.

  ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ adsorption and desorption rates are assumed to be comparable, and to occur with negligible isotopic fractionation, resulting in no net effect on the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool concentration or isotopic composition.

The model incorporates deliberate simplifications to reduce complexity, while remaining adaptable to new data or insights; however, it is acknowledged that these assumptions may significantly influence model outcomes and should be carefully considered when interpreting results.

2.4 Prior knowledge about model parameters

Model parameter values were derived from an extensive literature review, and formulated as prior distributions, as detailed and referenced in Appendix C: Prior values for inference. Positive parameters were parameterized as Lognormal priors, while priors of positive or negative parameters were parameterized as Normal distributions. Mean values were derived from the provided references, standard deviations were assigned either as absolute values or as percentages of the mean, depending on the class of variables. For parameters that are lake-specific (see model assumption iii.) and expected to be well identifiable from data, such as the maximum conversion rates of various processes (i.e., aerobic mineralization, the first step of nitrification, the first step of denitrification, mineralization by ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction, anaerobic mineralization) and the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ flux from deeper sediment layers, only limited prior knowledge is available, making the use of uniform priors preferable. As their interpretability can be questionable, uniform priors were applied only to parameters expected to be well-identifiable, ensuring that prior variations within the marginal posterior range would remain small, even with alternative broad priors. This approach avoids specifying typical expected values, while maintaining robust identifiability. The maximum conversion rates for anammox, DNRA, as well as the second step of nitrification and the second and third steps of denitrification (Anam, DNRA1, DNRA2, Nit2, Den2 and Den3) were more challenging to identify from data, as the sensitivity of model results to these parameters becomes very low when the concentration of the converted substance becomes small. Additionally, prior specification for these rates was difficult, due to the expected variability among different lakes, similar to other maximum conversion rate parameters. Therefore, their priors were formulated as ratios relative to the better-constrained maximum conversion rate of the first nitrification (i.e., k_Nit1) or denitrification step (i.e., k_Den1). This approach allowed for the characterization of the relative importance of each process without requiring absolute rate values. The joint prior for all parameters was assumed to be an independent combination of their respective marginal prior distributions.

2.5 Model-based analysis process

To partially reduce structural uncertainty of the model and to account for parameter non-identifiability, Bayesian inference was applied, considering all uncertain parameters listed in Appendix C: Prior values for inference. Some parameters were excluded from this analysis, including molecular diffusion coefficients, compound concentrations at the sediment surface, zero fluxes from deeper sediment layers (except for the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ flux, which was inferred jointly with other parameters) and bioturbation. These values are considerably less uncertain than the other model parameters, except for bioturbation, which was addressed separately through a scenario analysis, following Bayesian inference under the Base scenario.

The posterior distribution (probability density) of the model parameters, f_post, is expressed as

$$\begin{array}{}\text{(5)}& f_{\text{post}}\left(\mathit{\theta }\right)={\displaystyle \frac{f_{\mathrm{L}}\left(C\mathrm{|}\mathit{\theta }\right)f_{\text{pri}}\left(\mathit{\theta }\right)}{\int f_{\mathrm{L}}\left(C\mathrm{|}{\mathit{\theta }}^{\prime }\right)f_{\text{pri}}\left({\mathit{\theta }}^{\prime }\right)\mathrm{d}{\mathit{\theta }}^{\prime }}}\end{array}$$

where f_pri is the prior distribution (probability density) of the model parameters, f_L(C|θ) is the likelihood function of the model, C represents the observed compound concentrations, or δ¹⁵N values, and θ denotes the model parameters. The likelihood function f_L(C|θ) is defined as a multivariate, uncorrelated Normal distribution with constant variances (standard deviation, σ_δ) for δ¹⁵N values, and variances increasing linearly with concentration, leading to a standard deviation ${\mathit{\sigma }}_C=\sqrt{{\mathit{\sigma }}_{C,a}C+{\mathit{\sigma }}_{C,b}^{\mathrm{2}}}$ for O₂, ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$, and N compound concentrations. This formulation incorporates the combined uncertainties in model structure, sampling, and concentration measurements. To account for the unknown magnitude of these uncertainties, the coefficients of these relationships, σ_C,a, σ_C,b, and σ_δ, were inferred alongside the model parameters.

The marginal posteriors of individual parameters were compared with their priors to evaluate whether observational data provided information about these parameters, and whether this information was in conflict with the priors. In addition, two-dimensional marginals were examined to identify potential identifiability issues. Finally, uncertainty in the model results was calculated by propagating parameter uncertainty to the model results under consideration of their uncertainty for given parameter values as formulated in the likelihood function:

$$\begin{array}{}\text{(6)}& f_{\text{post}}\left(C\right)=\int f_{\mathrm{L}}\left(C\mathrm{|}\mathit{\theta }\right)f_{\text{post}}\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\end{array}$$

For the parameters with marginal posteriors in conflict with prior information, we conducted additional scenario analyses, fixing parameters, and narrowing or widening prior distributions. These analyses evaluated the model's compatibility with observational data if parameters better aligned with prior information and assessed changes in posterior distribution with weaker priors. These scenario analyses complemented the assessment of bioturbation uncertainty mentioned above.

2.6 Discretization and numerical algorithms

The partial differential equations outlined in Appendix B: Reaction-diffusion model were solved using the Method of Lines. For spatial discretization, a grid was employed with cell thickness increasing progressively from the sediment surface toward deeper layers. This adaptive grid design reduced the total number of cells required, while still maintaining high resolution near the sediment-water interface, where steep concentration gradients typically occur (Appendix D: Model discretization). The resulting system of ordinary differential equations (ODE) was solved by a standard ODE solver. Parameter inference was conducted using two advanced Bayesian inference algorithms: Metropolis (Andrieu et al., 2003; Vihola, 2012) and Hamiltonian Monte Carlo (Betancourt, 2017; Neal, 2011) algorithms.

2.7 Model implementation

The model was implemented in Julia (Bezanson et al., 2017) (https://julialang.org, last access: 11 July 2024) to achieve high-performance and facilitate automatic differentiation. The DifferentialEquations.jl package (Rackauckas and Nie, 2017) was used to solve the system of ODEs; performance testing of several ODE solvers identified the FBDF solver (adaptive order and adaptive time-step backward-differencing solver) as the most suitable for handling the stiffness of the ODE system. The ForwardDiff.jl package (Revels et al., 2016) was used for automatic differentiation; Bayesian inference was conducted using the adaptive Metropolis sampler from the AdaptiveMCMC package (Vihola, 2020), and the Hamiltonian Monte Carlo algorithm implemented in the AdvancedHMC.jl package (Xu et al., 2020). Further implementation details are provided in Appendix E: Model implementation. Simulations were performed at sciCORE (https://scicore.unibas.ch, last access: 26 February 2025), the scientific computing centre at the University of Basel.

3 Sample collection and analyses

3.1 DIN concentrations and isotopes

Sediment cores were retrieved at the deepest location of the Kreuztrichter basin in Lake Lucerne, a large oligotrophic lake in Switzerland (Baumann et al., 2024), in April 2021 using a gravity corer with PVC liners. The sediment cores were stored at 4 °C and processed using two porewater-sampling methods: whole-core squeezing (WCS; Bender et al., 1987) for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ samples, and Rhizon samplers (Rhizosphere research products, Wageningen, NL) for ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ samples. The WCS technique provides a high depth resolution near the sediment-water interface (0–5 cm, resolution: ∼ 0.7–1 mm), where ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ is present in porewaters, while the Rhizon sampling method allows collecting samples at greater sediment depths (> 5 cm, resolution: ≥ 0.5 cm). ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentrations were measured using ion chromatography (940 Professional IC Vario, Metrohm). δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ were determined using the denitrifier method (Casciotti et al., 2002; Sigman et al., 2001), and the hypobromite-azide method (Zhang et al., 2007), respectively. In both methods, sample N from ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ or ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ is converted into N₂O, which is then purified and analysed by isotope ratio mass spectrometry (Delta V Plus, Thermo Fisher Scientific). The typical analytical precision is ∼ 0.25 ‰ (McIlvin and Casciotti, 2010).

3.2 Process rate measurements

For model parameterization, reaction rates for denitrification, DNRA, and anammox were determined using established protocols for ¹⁵N-tracer incubations (Holtappels et al., 2011). After recovery and sectioning of the core into 1 cm intervals, 1 g of sediment was placed into 12 mL gas-tight glass vials (Exetainers^®, Labo, UK). These Exetainers were then filled with anoxic, sterilized bottom water, amended with the following tracers: (Exp1) ${}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}$, (Exp2) ${}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}+{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$. Exetainers were incubated at 6 °C in the dark, and terminated at designated time points (0, 6, 12, 24, and 36 h) by adding ZnCl₂. Gas headspace samples were analysed for the production of ¹⁴N¹⁵N and ¹⁵N¹⁵N using gas-chromatography isotope ratio mass spectrometry (GC-IRMS; Isoprime, Manchester, UK). Linear regression of ¹⁴N¹⁵N and ¹⁵N¹⁵N production over time was used to calculate N₂ production rates, with standard errors derived from deviations in the regression slopes across the five-time points. For the determination of ${}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$ production from ${}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}$ additions, ${}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$ was chemically converted to N₂ gas using the alkaline-hypobromite method (Jensen et al., 2011). The resulting ¹⁴N¹⁵N was quantified by GC-IRMS. Linear regression of ¹⁴N¹⁵N production over time was used to calculate potential rates of ²⁹N₂ (i.e., ${}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}$) production. Rates of denitrification, DNRA, and anammox were calculated according to Holtappels et al. (2011) and Risgaard-Petersen et al. (2003). Only data from the upper 1 cm were used to parameterize the model, as the investigated sediments displayed a shallow nitracline and the highest anammox contribution at 0–0.5 cm depth.

4 Results and Discussion

The developed diagenetic N isotope model addresses existing knowledge gaps in understanding porewater N dynamics, and aims to clarify the roles of distinct N-transformation processes in shaping the distribution of N isotopes to be potentially used to constrain benthic N (isotope) fluxes across different environments. Here, we present (1) the results of Bayesian inference applied to a large number (∼ 60) of model parameters (see prior definition in Appendix C: Prior values for inference), with a focus on assessing their uncertainty, (2) a detailed scenario analysis, focusing on parameters that exhibit significant shifts in their marginal posterior distributions relative to their prior, as well as on the effect of variable contributions from different ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction pathways, and the impact of enhanced bioturbation on model outcomes, (3) a sensitivity analysis, evaluating the importance of individual model processes in shaping benthic N isotope dynamics, (1) the importance of process coupling in benthic N cycling, with a particular focus on the role of intermediate ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ in influencing δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ dynamics. All results are based on porewater concentration, isotope, and rate measurement data from a sampling campaign conducted in Lake Lucerne in April 2021. Additionally, we performed (2) a sensitivity analysis examining model output responses to modifications of selected parameters using artificially simulated settings (e.g., variable contributions of denitrification/anammox/DNRA); this analysis demonstrates the model's capability for addressing diverse research questions.

4.1 Bayesian inference

The model implementation was highly efficient, achieving simulation times of about 12 s on a 13th Gen Intel^® Core™ i9–13 900 K processor with 3.00 GHz and 64 GB of memory (of which only a small fraction was needed) for a 100 d simulation starting from constant concentration profiles. This efficiency enabled the execution of Markov chains of 20 000 iterations within a few days on the scientific computing centre at the University of Basel (https://scicore.unibas.ch, last access: 26 February 2025). By combining these chains, samples of 100 000 iterations were generated. The Hamiltonian Monte Carlo algorithm outperformed the adaptive Metropolis algorithm during burn-in to the core of the posterior distribution. However, for final posterior sampling with about 60 parameters, adaptive Metropolis sampling proved more efficient in terms of effective sample size per unit of simulation time. Despite these efforts in getting computational efficiency, and the use of advanced MCMC algorithms, reaching convergence of the Markov chains remained challenging. We got five consistent Markov chains without discernible trends for each scenario; however, some widening of the chains and the resulting effective sample size on the order of 500 indicate that we were not able to get a good coverage of the tails of the posterior distribution. This outcome demonstrates that incorporating so many uncertain model parameters pushes the limits of Bayesian inference in terms of numerical tractability. However, the resulting uncertainty estimates are certainly more realistic than those obtained by fixing many poorly constrained parameters to unique values to reduce the dimension of the parameter space.

Figure 2 Vertical porewater profiles of concentrations (a, b) and isotopic composition (δ¹⁵N) (c) of the state variables for the Base scenario. Continuous lines represent model simulations, while symbols represent observational data from Lake Lucerne. For ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentrations, filled diamonds represent low-resolution data from Rhizon sampling, while open diamonds represent the high-resolution WCS data, adjusted to align with absolute concentrations measured in the low-resolution dataset. Dashed lines enclose 95 % credibility intervals resulting from parametric uncertainty, while thin solid lines represent total uncertainty.

The simulation results of solute concentration and δ¹⁵N profiles in the most plausible Base scenario (Fig. 2) integrate prior knowledge (Appendix C: Prior values for inference) with observational data through Bayesian inference. The profiles closely reproduce the available, albeit limited, data, and conform to expected depth-related trends: oxidants (i.e., O₂, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$) are readily consumed via aerobic mineralization and nitrification (O₂), denitrification (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$), and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction. While mineralization is assumed to involve negligible N isotopic fractionation, the first step of nitrification causes significant enrichment in ¹⁵N of the residual ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool, yielding δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ values up to 11.2 ‰ at 0.15 cm, due to strong N isotope fractionation, estimated at ε_Nit1 = 12.0 ‰ (to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$) and 36.4 ‰ (to N₂O). Unfortunately, extremely low ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentrations measured in the top 2 cm hindered the determination and verification of the modelled δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ in this zone with field data. Both ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂O accumulate in the upper 0.5 cm, reaching up to 0.4 and 2 µM, respectively. Below 0.3 cm, denitrification leads to the progressive ¹⁵N enrichment of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂O, while N₂-producing mechanisms (i.e., denitrification and anammox) cause only minimal changes to the modelled δ¹⁵N-N₂ profile, due to the dominance of a large pre-existing N₂ pool. For concentrations, the 95 % credibility intervals of parametric uncertainty are rather narrow, whereas the much broader total uncertainty is dominated by the lumped uncertainty term in the likelihood function, which primarily reflects the model's structural uncertainty. The error, beyond the parameter error, is parameterized using the two sigma values (σ_C,a and σ_C,b; see Sect. 2.5), and exceeds what would arise from measurement and sampling alone. This suggests that the larger error is attributable to the model's structural limitations. Conversely, δ¹⁵N profiles exhibit small total uncertainty, as model results for δ¹⁵N closely match observational data, with minimal random and systematic deviations (parameterized using the sigma value σ_δ, see Sect. 2.5).

Figure 3 Vertical profiles of transformation rates for distinct N-cycling processes affecting the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, and N₂O pools. Dashed lines enclose 95 % credibility intervals resulting from parametric uncertainty. Positive reaction rate values indicate production, negative values indicate consumption of a given DIN species.

Figure 4 Prior (dashed line) and posterior marginal distributions (continuous line) for illustrative parameters, which could be identified and showed (a) good (f_Anam,side) and (b) poor agreement (ε_Den1) with prior knowledge, and (c) for parameters that could not be identified ($a_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}$); 2D correlation plot for ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}},\text{MinSulfRed}}$ versus $F_{{\mathrm{NH}}_{\mathrm{4}}}$ (d).

The model provides insights into the underlying process rates (Fig. 3) that shape the simulated profiles (Fig. 2). Vertical profiles of transformation rates for ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂O clearly illustrate the sequential dominance of different N-transformation processes with increasing sediment depth and decreasing O₂ availability. Aerobic processes, namely aerobic mineralization and nitrification, primarily control ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ transformation rates, peaking at 450 and 350 µM d⁻¹, respectively (Fig. 3a). Nitrification sustains denitrification by producing both ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (up to 350 µM d⁻¹) and ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ (up to 275 µM d⁻¹) in the upper 0.4 cm (Fig. 3b and c). A strong spatial overlap of nitrification and denitrification emerges in the depth distribution of processes affecting the ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ pool, suggesting a potential interplay between these pathways (Fig. 3c).

A key strength of this model is the incorporation of N₂O as a state variable. Our model results reveal that, although N₂O production via nitrification is minimal (not visible in Fig. 3d), the strong isotopic fractionation associated with this reaction (${\mathit{\epsilon }}_{\text{Nit1},{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ = 36.4 ‰) generates N₂O with δ¹⁵N values of −1.2 ‰ to −2.2 ‰ in the top 0.2 cm (Fig. 2c). At a depth of approximately 0.35 cm, up to 2.1 µM of N₂O accumulate, coinciding with the highest rates of N₂O production through denitrification. Conversely, N₂O consumption by the last denitrification step peaks at 0.5 cm, leading to a progressive increase in δ¹⁵N-N₂O with depth. This zonation likely reflects the O₂ sensitivity of the distinct N₂O-producing and -consuming processes. Specifically, N₂O reductases are known to be strongly inhibited by O₂, and therefore exhibit greater activity below the oxycline (Wenk et al., 2016). Although the model does not explicitly include the enzymes responsible for N-transformation pathways, the chosen and estimated kinetic parameters reflect substrate affinity and inhibition strength. Consequently, inhibition constants like $K_{{\mathrm{O}}_{\mathrm{2}},\text{Den2}}$ and $K_{{\mathrm{O}}_{\mathrm{2}},\text{Den3}}$ provide indirect insights into the O₂ dependency of these enzyme-mediated reactions, effectively shaping the modelled redox zonation.

The model adequately captures the concentration and isotopic composition of the state variables, in agreement with field measurement and the expected patterns of underlying N-transformation processes and reaction coupling (Figs. 2 and 3). One key strength of the step-wise model is its ability to quantify reaction coupling, which is challenging to infer directly from state variable pools (i.e., reactive intermediates), if they are rapidly turned over.

To address the variable ranges for the model parameters found in the literature, and to reduce structural uncertainty imposed by fixed parameter values, we estimated a large set of parameters using Bayesian inference. The obtained joint posterior distribution of model parameters enabled us to assess the knowledge acquired from data. Marginal posterior distributions of individual parameters, and two-dimensional marginal distributions of parameter pairs, were particularly useful in this context (Fig. 4 shows examples for the four categories defined below; Fig. S1 in the Supplement provides an overview of all marginal prior and posterior parameter distributions). By comparing marginal posterior distributions with their corresponding priors, parameters were classified as well identifiable or poorly identifiable. While this classification involves some subjectivity in determining how much narrower a posterior distribution should be compared to its prior distribution to classify such parameter as well identifiable, some clear patterns emerged:

• 1.

  Well identifiable parameters: The marginal posterior distribution is clearly narrower than the prior, indicating that data provide meaningful information about the parameter's value. Two cases were observed:

  • a.

    The marginal posterior distribution is within the prior range, suggesting that the information from the data is in agreement with prior knowledge (Fig. 4a). Examples include: f factors for anammox (f_Anam,Den2 = 0.2) and both DNRA steps (f_DNRA1,Den1 = 0.005, f_DNRA2,Den2 = 0.005), estimated using ¹⁵N-tracer incubation experiments for the investigated system, and parameters such as $K_{{\mathrm{NO}}_{\mathrm{3}},\text{Den1}}$ and $K_{{\mathrm{O}}_{\mathrm{2}},\text{MinOx}}$, constrained from clearly defined oxidant declines. Maximum conversion rates for aerobic mineralization, denitrification, ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction, and anaerobic mineralization, as well as the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ flux from deeper sediment layers, also belong to this category, although we approximated very wide priors by uniform priors (see Sect. 2.4), making it less visible in the plot.

  • b.

    The marginal posterior distribution significantly deviated from the prior range (Fig. 4b), suggesting that the information from the data is in conflict with prior knowledge. The most striking example is ε_Den1, estimated at 2.8 ± 1.1 ‰ for the Lake Lucerne dataset, far lower than the typical 15 ‰–25 ‰ reported in the literature for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction (Lehmann et al., 2003; Rooze and Meile, 2016), suggesting a reduced N-isotopic fractionation (or at least, of its expression) at the porewater level. This finding contrasts with model-derived values for the cellular isotope effect of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction observed in the porewater of marine sediments (ε_Den > 10 ‰) (Lehmann et al., 2007). While a detailed investigation of the biological mechanisms behind such reduced expression across benthic environments is beyond the scope of this study and will be addressed separately by the authors, the potential role of reaction couplings in modulating benthic N isotope dynamics is discussed in Sect. 4.4.

• 2.

  Poorly identifiable parameters: The marginal posterior distribution resembles the prior distribution, suggesting poor identifiability. This can occur for two possible reasons:

  • a.

    The parameter exerts negligible influence on the model output that corresponds to observational data (Fig. 4c). For example, parameters like the N₂O yield during nitrification, $a_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}$ and $b_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}$, could not be constrained without specific data on N₂O production. The current model encompasses several processes and state variables, which, at times, were hard to corroborate with the limited dataset in hand (a situation that may apply regularly to environmental studies, particularly in benthic environments). Therefore, their values were taken from previous studies (Ji et al., 2018). For other parameters, such as ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}},\text{DNRA1}}$ and ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}},\text{DNRA2}}$, little knowledge was acquired from the data in hand, due to the relatively low maximum rates of DNRA compared to other processes. In such cases, the posterior distribution may remain close to the prior, not because the prior range was incorrect, but because the available data could not further constrain it.

  • b.

    Although data are available and the model output is sensitive to the parameter, other parameters influence the output similarly. This leads to parameter correlation in the posterior distribution and reduces identifiability, as observed for ${\mathit{\gamma }}_{{\mathrm{NH}}_{\mathrm{4}},\text{MinSulfRed}}$ and $F_{{\mathrm{NH}}_{\mathrm{4}}}$ (Fig. 4d), which exhibit correlation, making their estimates interdependent (Guillaume et al., 2019). Here, the estimate of the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ flux from the lower boundary of the model depends on the estimate of the amount of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ released via OM mineralization coupled to ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction.

The comparison of marginal priors and posteriors of the parameters (Fig. S1) demonstrates that excellent agreement between model outputs and observational data (Fig. 2) can be achieved for 54 of the 58 estimated parameters compatible with their priors. Exceptions include: the higher-than-expected rate for the second denitrification step relative to the first (expressed by the factor f_Den2,Den1), the large half-saturation constant for ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction ($K_{{\mathrm{SO}}_{\mathrm{4}},\text{MinSulfRed}}$), and smaller-than-expected N isotope effects for the first steps of denitrification and nitrification (ε_Den1 and ${\mathit{\epsilon }}_{\text{Nit1},{{\mathrm{NO}}_{\mathrm{2}}}^{-}}$, respectively). The largest deviation is observed for ε_Den1, which is further examined in the next subsection.

Notably, the seven parameters, for which a uniform prior was chosen to approximate a very wide prior (k_MinOx, k_Den1, k_MinSulfRed, k_MinAnae, k_Nit1, $F_{{\mathrm{NH}}_{\mathrm{4}}}$, δ¹⁵N, $F_{{\mathrm{NH}}_{\mathrm{4}}}$), were identifiable, indicating that highly system-specific prior knowledge is not crucial for these estimates. Most of the other model parameters showed limited narrowing of the marginal posterior relative to the prior, reflecting the rather limited information gain that can be obtained from data. The three model error parameters (σ_C,a, σ_C,b, σ_δ) were well identifiable and will be used in the following sections to compare the fit quality across different modelling scenarios.

4.2 Scenario analysis

Building on the findings discussed in the previous section, we explored the apparent prior-data conflict regarding ε_Den1 in greater detail. Additionally, we assessed whether the estimated process rates overlooked potential reaction coupling, which might go undetected through ¹⁵N-tracer incubation experiments, by exploring the variability in contributions of anammox and DNRA (i.e., f_Anam, f_DNRA1 and f_DNRA2). Lastly, given the uncertainty regarding solute-diffusion enhancement by bioturbation, we investigated a scenario with increased bioturbation. These considerations led to four key scenarios:

• A.

  Narrow priors for ε. This scenario investigated the effects of restricting ε variability to a narrower range (prior standard deviation of 1 ‰ instead of 5 ‰). The aim was to test whether the marked reduction in the marginal posterior of ε_Den1 persisted under stricter prior assumptions, and whether this decreased flexibility significantly impacted the quality of the model fit.

• B.

  Fixed ε. Here, the model output was assessed under the assumption that the literature data regarding N isotope effects are correct (i.e., ε values not estimated). This scenario complemented Scenario (A) by testing whether a good fit to the data could still be achieved by fixing the ε_Den1 value (and all other isotope effects) at its prior mean.

• C.

  Wider priors for f. In this scenario, greater variability in DNRA and anammox contributions (prior standard deviation of 100 % instead of 25 %) was allowed to test the impact of relaxed prior assumptions on the relative contributions of these processes in the model output.

• D.

  Enhanced bioturbation. This scenario simulated a faster solute-diffusive transport due to higher infaunal activity by doubling the bioturbation coefficient (D_bio = 2 cm² d⁻¹ instead of 1 cm² d⁻¹), to investigate the sensitivity of the results to this uncertain parameter, which was not included in the Bayesian analysis. In the model, the bioturbation strength at the sediment surface is defined by the parameter D_bio, and it decreases exponentially with depth, with the typical bioturbation depth parameter, depth_bio. As the diffusion enhancement by bioturbation is highly uncertain, this scenario aims to assess solely the sensitivity of the model output to changing bioturbation magnitude.

The results demonstrate a strong dependence of the estimated parameters on the chosen prior assumptions (Fig. 5). Across all scenarios, marginal posterior distributions for the selected parameters are generally narrower than the prior distributions, though results vary substantially. In Scenario (A) (Narrow priors for ε), restricting the prior range significantly constrained ε_Den1, limiting its deviation from the prior (Fig. 5m; note that the prior for Scenario (A) is 5 times narrower than the one shown, which represents the prior for all other scenarios). These results closely resemble those from Scenario B (Fixed ε), where no deviation was possible (Figs. 5 and S2 in the Supplement). Both scenarios exhibit lower denitrification rates than the Base scenario (Fig. 5b), but comparable fit quality for total (¹⁴N+¹⁵N) concentration, quantified by σ_C,a (i.e., the dominant term of standard deviation of the model error for concentrations, see Sect. 2.5) (Fig. 5x). On the other hand, scenarios (A) and (B) display poorer fit quality for δ¹⁵N profiles, indicated by a large value of σ_δ (Fig. 5z), suggesting that the model structure cannot adequately reproduce the δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ profiles without adapting the ε_Den1 value. While biological isotope effects of 15 ‰–30 ‰ are typical for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction (Lehmann et al., 2007), lower values under almost-complete ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ consumption have been reported (Thunell et al., 2004; Wenk et al., 2014). This finding is further confirmed by comparable marginal posteriors for ε_Den1 across all scenarios considered in this study, besides scenarios (A) and (B). To test the robustness of our model, we ran a base scenario simulation for marine sediments in the Bering Sea (station MC16) (Lehmann et al., 2007) (data not shown). Moreover, a manuscript currently in preparation presents an extensive comparison of model application across different sites and demonstrates a much wider range of ¹⁵ε_Den1 values, exceeding 20 ‰.

Figure 5 Marginal probability densities across the five considered scenarios for selected estimated parameters, showing both prior (dashed line) and posterior distributions (continuous lines): Base scenario (SD_f = 25 %, SD_ε = 5 ‰, D_bio = 1 cm² d⁻¹), Narrower ε (SD_ε = 1 ‰), Fixed ε (i.e., ε taken from bibliography), Wider  f (SD_f = 100 %) and Enhanced bioturbation (D_bio = 2.0 cm² d⁻¹). Of the ∼ 60 estimated parameters, those shown here were selected for their relevance to the discussion. See main text for further details.

In Scenario (C) (Wider  f), allowing greater variability in anammox and DNRA contributions results in the lowest f_Anam,Den2 values, although such deviation is not substantial compared to the Base scenario output (Fig. 5i). The estimated f_DNRA1,Den1 and f_DNRA2,Den2 values in Scenario (C) mostly align with those of the Base scenario, corroborating the marginal role of DNRA in Lake Lucerne. Such findings confirm the accuracy of the rate measurements performed with ¹⁵N tracer incubations.

Scenario (D) (Enhanced bioturbation) stands out with the highest conversion rates (i.e., k_MinOx, k_MinSulfRed, and k_Nit1) (Fig. 5a, e, and g) to ensure sufficient oxidant consumption at higher supply/flux rates (reproducing the observed gradient despite higher diffusivity). Despite these changes, bioturbation had negligible effects on porewater N isotope dynamics, with estimated isotope effects and fit quality for δ¹⁵N profiles (σ_δ) comparable to those of the Base scenario.

The obtained concentration depth profiles for the four scenarios are generally comparable, as newly estimated parameters ensured good fitting of the data (Fig. S2). However, in Scenarios A and B, stricter constraints on prior knowledge for parameter estimation result in little to no suppression of all isotope effects (i.e., relatively strong N isotopic fractionation), leading to great variability in the δ¹⁵N profiles. Poor fits to the δ¹⁵N data are observed under these conditions, as evidenced by the greater ¹⁵N enrichment of the ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ pool compared to the measured-data profiles (Fig. S2). Similarly, the δ¹⁵N-N₂O profiles exhibit sharp declines to approximately −15 ‰ in the upper 0.5 cm under scenarios (A) and (B), driven by the strong expression of ${\mathit{\epsilon }}_{\text{Nit1},{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}$ (40.1 ‰ and 40.0 ‰, respectively). In contrast, Scenarios (C) and (D) closely resemble the Base scenario, with only minor δ¹⁵N-N₂O variations.

Figure 6 Vertical concentration (a–d) and isotopic composition (e–h) profiles for state variables. Model output obtained with all processes included (a, e) are compared with model simulations where individual processes are switched off: nitrification (b, f), anammox (c, g), and DNRA (d, h), without running inference again. Continuous lines represent the model output, while symbols represent measured data from Lake Lucerne. For ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, open diamonds represent the high-resolution dataset, adjusted to align with absolute concentrations measured in the low-resolution dataset (filled diamonds).

Figure 7 Posterior marginal probability distributions of modelled sediment-water interface fluxes (in $\mathrm{nmol}{\mathrm{cm}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) for all state variables, generated from inference runs, across the four scenarios considered for model validation against experimental data from Lake Lucerne.

4.3 Importance of modelled processes and their impact on porewater N isotope signatures

The importance of modelled processes and their impact on N isotope signatures were investigated by selectively deactivating individual processes and comparing the model outputs to the Base scenario. Aerobic mineralization, denitrification, and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$ reduction were considered essential to preserve redox zonation (e.g., sequential decline of O₂, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, and ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$) and N dynamics. The following processes were individually turned off: (a) nitrification (“NitOff”); (b) anammox (“AnamOff”); and (c) DNRA (“DNRAOff”). Initially, each process was simply inactivated to assess its impact on model outputs (Fig. 6). Subsequently, inference was conducted after deactivating each process, to investigate their importance for model performance, parameter and flux estimation, and for the identifiability of rate parameters by evaluating the quality of the fit to the data, especially on the δ¹⁵N profiles (Fig. 7, Figs. S3 and S4 in the Supplement).

Switching off nitrification significantly alters the model output compared to the Base scenario (Fig. 6a, b, e, and f), indicating its central role in the benthic N dynamics. Key effects include ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ accumulation throughout the investigated depths, with a flattening of the δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ profile (i.e., less curvature towards higher δ¹⁵N values) in the upper 0.5 cm, as the only other source of ¹⁵N-enriched ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ besides nitrification would be anammox, which is inhibited under oxic conditions. Furthermore, nitrification-denitrification coupling via ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ weakens in this scenario, resulting in lower overall N₂ production (as indicated by the lower maximum N₂ concentration of 734 µM compared to 745 µM in the Base scenario). These results suggest that partially reducing, or fully eliminating, nitrification lowers the system's capacity to act as an efficient N sink. In other words, the findings confirm that nitrification is a critical process that, when closely coupled to denitrification, helps to enhance the ecosystem's potential to remove fixed N. All other N-isotopic state variables also show a flatter δ¹⁵N profile, with only a progressive enrichment in ¹⁵N below 0.5 cm, primarily driven by denitrification (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, and N₂O). The impact of disabling nitrification is clearly reflected in the δ¹⁵N-N₂O profile across the upper 0.3 cm, where the typical nitrification-induced dip is absent, and δ¹⁵N-N₂O values remain relatively constant (∼ 7 ‰–8 ‰). In contrast, the effects of turning off anammox or DNRA are more subtle, owing to their generally lower reaction rates in Lake Lucerne (Fig. 6c, d, g, and h). Notably, in the absence of anammox, N₂O exhibits lower δ¹⁵N values in the upper 0.3 cm compared to the Base scenario, likely due to higher N₂O yields via nitrification, as reduced competition for ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ with anammox provides more substrate for nitrification.

Upon running inference for each case, concentration and N isotope profiles for the NitOff, AnamOff, and DNRAOff scenarios are generally similar to those of the Base scenario (Fig. S3), with notable exceptions in the NitOff case. In the absence of nitrification, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ accumulates and the δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ profile remains largely flat, since anammox, the only other ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$-consuming process, is minimal under oxic conditions. No δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ measurements are available for the top 1 cm, so the model output could not be verified with field data. The N₂O pool systematics also diverge between the NitOff and Base scenarios. Specifically, in the NitOff case, no nitrification-derived N₂O accumulates in the upper 0.4 cm, and consequently, the δ¹⁵N-N₂O profiles lacks the typical nitrification-associated decline in this layer. Instead, N₂O becomes progressively enriched in ¹⁵N below 0.4 cm. While most estimated parameters and fluxes are consistent across the four scenarios, the NitOff scenario stands out again, exhibiting strong effects on the anammox rates and associated isotope effects (e.g., f_Anam,Den2, ${\mathit{\epsilon }}_{\text{Anam},{\mathrm{NH}}_{\mathrm{4}}}$) (Fig. S4), as well as on benthic fluxes of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and N₂O (Fig. 7). Nonetheless, the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentration profile is well-captured, as indicated by a low σ_C,a, reflecting a good match between model and concentration data even in the absence of nitrification. This finding implies that the model cannot resolve the relative contributions of nitrification versus anammox to ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ consumption based on the concentration and isotope data, highlighting the importance of prior knowledge regarding f_Anam,Den2.

The comparison of process rates across these four scenarios provides insights, unveiling the extent of process coupling and competition (Fig. S5 in the Supplement) (Hines et al., 2012). For instance, anammox and nitrification compete for both ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ as substrates, causing the rate of one process to be enhanced, when the other is switched off. The ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ oxidation and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ production rates via nitrification (Nit1) are higher (∼ 0.2 cm depth) in the AnamOff scenario than in the Base scenario. Even more obviously, enhanced rates of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ oxidation, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ consumption, and ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ production via anammox are observed in the NitOff scenario than in the Base scenario. Process coupling, specifically nitrification-denitrification, is further confirmed by lower rates for ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction via denitrification (Den2) in the absence of nitrification. In general, the influence of DNRA on production and consumption rates of the considered state variable appears minimal, owing to the limited environmental relevance of DNRA in Lake Lucerne. Overall, the similarly good fits obtained across these three scenarios and the Base scenario reflect the poor identifiability of the switched off processes. This suggests that the data can be well-fitted even without these three processes, emphasizing the importance of prior knowledge about their environmental relevance.

4.4 The role of process coupling via ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$

Previous models of benthic N isotope dynamics have focused on individual reactions or overlooked the role of intermediate species, such as ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (Kessler et al., 2014; Lehmann et al., 2007). Our study confirms that ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ plays a critical role in coupling multiple N-transformation processes and shaping benthic N isotope dynamics, including that of δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$. While such process coupling has been examined in the water column (Frey et al., 2014), it remains, to our knowledge, largely unexplored in sedimentary environments.

Figure 8 Depth profiles of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ concentrations and N isotopic composition (a, c), and rates of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$-producing and -consuming processes (b, d), as simulated by the Base scenario (a, b), and the one-step denitrification approach (c, d). In the one-step approach, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ is reduced directly to N₂, omitting ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ as an intermediate; thus, no ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ is produced or consumed through denitrification. Dashed lines enclose 95 % credibility intervals resulting from parametric uncertainty.

To assess the significance of this coupling, we implemented a one-step denitrification approach that bypasses ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ as an intermediate, replacing the three-step pathway used throughout this paper (Fig. 8). In this simplified model, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ concentrations and isotopic signatures are shaped solely by nitrification (and to a marginal extent, DNRA and anammox), as denitrification no longer contributes to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ production. This modification leads to significantly reduced ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ accumulation, restricted to the upper 0.3 cm, and lower anammox activity, due to a lack of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ substrate below the oxycline. The absence of denitrification-derived ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ has profound effects on the N isotope dynamics. First, a consistent ∼ 15 ‰ offset between δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and δ¹⁵N-${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ is evident across all modelled depths (Fig. 8c). This offset is ascribed to the isotope effect of the second nitrification step (ε_Nit2 = −13.7 ‰), and the lack of ¹⁵N enrichment in the ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ pool from denitrification. Second, the estimated isotope effect for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction (ε_Den) increases to 5.5 ± 0.9 ‰, nearly double than in the Base scenario, indicating that elevated δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ values in the field data may, to some extent, reflect ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ isotope dynamics, rather than solely the effect of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction (Fig. 1).

These findings emphasise the importance of both ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$-producing and -consuming processes in modulating δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, and consequently, estimates of ε_Den1. Although nitrification is typically aerobic and denitrification anaerobic, evidence exists that indicates spatial overlap of these two processes at the bottom of oxyclines in natural aquatic environments (Frey et al., 2014; Granger and Wankel, 2016). In this transition zone, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ produced by either pathway can be oxidised to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ or reduced to N₂O, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ or N₂ (Fig. 3), significantly affecting its δ¹⁵N signature (depending on the N-branching). For instance, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction to N₂O enriches the residual ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ pool in ¹⁵N. If this ¹⁵N-enriched ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ is subsequently oxidized to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ (a reaction that exhibits an inverse kinetic isotope effect), the resulting ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ will be markedly enriched in ¹⁵N (Fig. 1). Such interactions have been shown to influence apparent isotope effects for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ in the water column (Frey et al., 2014), and likely exert similar effects in sediments, where sharp redox gradients create overlapping zones of nitrification and denitrification. This coupling may explain the discrepancy in estimated ε_Den1 values between the Base scenario (2.8 ± 1.1 ‰) and the one-step denitrification model approach (5.5 ± 0.9 ‰).

Anammox further complicates these dynamics, as it depends on ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ excreted into the environment. Without denitrification, which releases ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (Sun et al., 2024), anammox is substrate limited (Fig. 8). Thus, while previous benthic studies estimated denitrification isotope effects using one-step denitrification approaches (Lehmann et al., 2007), our findings call for the adoption of a stepwise modelling approach (Sun et al., 2024) that better captures the interdependence of N-transformation pathways, and their integrated effects on ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ isotope dynamics. A more detailed examination of these interactions is essential for refining our understanding and quantification of isotope effects associated with ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction in sedimentary systems.

Figure 9 Depth profiles of process rates, solute concentrations and δ¹⁵N values for the two idealized case scenarios investigated: (i) ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction via DNRA and denitrification (a–d), (ii) N₂ production via anammox and denitrification (e–h). Shadings represent different model scenarios within each case, as defined in the legend. For case (i), colour shading lightens with increasing contribution of DNRA (relative to denitrification) to total ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction. DNRA accounts for 0 % (f_DNRA=0), 33 % (f_DNRA=0.5), 50 % (f_DNRA=1) and 66 % (f_DNRA=2) of total ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction (a). The resulting effects on the production rates of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and N₂ (b), as well as on their concentrations and N isotopic composition (c, d), are shown. For case (ii), colour shading lightens with increasing contribution of anammox (relative to denitrification) to total ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ consumption and associated N₂ production. Anammox contributes 0 % (f_Anam=0), 33 % (f_Anam=0.5), 50 % (f_Anam=1) and 66 % (f_Anam=2) of total ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ consumption (e, f). The resulting impacts on N₂ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ concentrations and δ¹⁵N values are shown in panels (g) and (h).

4.5 Model applicability in distinct scenarios

Beyond applying and testing the developed diagenetic N isotope model at our site of interest (Lake Lucerne), we believe its strength hinges on its versatility to address distinct research questions and objectives. We explored two scenarios as examples of how the model can be adapted to provide insights into the N cycle in benthic environments and the N isotopic fingerprints that the combined N-cycling processes leave behind (Fig. 9). Understanding these fingerprints and how they might be modulated in natural environments (e.g., through the variable balance between individual processes constrained by environmental conditions) is important for correctly interpreting the distribution of ${}^{\mathrm{15}}\mathrm{N}{/}^{\mathrm{14}}\mathrm{N}$ ratios in N species as biogeochemical tracer, helping to pinpoint and disentangle individual N-turnover processes where they co-occur.

For comparison purposes, we used the estimated parameters from the Base scenario and modified the relative importance of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ or ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction via (i) denitrification vs. DNRA, and (ii) denitrification vs. anammox. This was done by progressively increasing the factors that define the contributions of DNRA (f_DNRA1,Den1 and f_DNRA2,Den2) and anammox (f_Anam,Den2) from 0 (i.e., no DNRA/anammox) to 2 (corresponding to DNRA and anammox accounting for $\mathrm{2}/\mathrm{3}$ of the total ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction, respectively). Simultaneously, the rates of the first two steps of denitrification (k_Den1 and f_Den2,Den1) were adjusted to maintain consistent overall ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction rates across scenarios. These model results were not validated against observational data and should therefore be considered as illustrative examples of the model's sensitivity to selected parameters, rather than as predictions with direct environmental relevance.

• i.

  N removal versus N retention. The model results confirm the spatial co-occurrence of DNRA and denitrification, with peak ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ (data not shown) and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (Fig. 9a) reduction activities localized between 0.4–0.6 cm depth. In contrast, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and N₂ production exhibit subtle differences in depth distribution: ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ production via DNRA extends across a broader sediment layer than N₂ production via denitrification (Fig. 9b). This pattern likely reflects the inhibitory effect of O₂ on N₂O reduction, the final denitrification step, pushing N₂ production to deeper, anoxic layers below the oxycline.

  Reduction of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ exhibits distinct isotope effects depending on the pathway: denitrification (ε_Den1≈ 2.8 ± 1.1 ‰) and DNRA (ε_DNRA1≈ 20.0 ± 2.9 ‰), according to our model estimates (Fig. 5m and v). This large difference reflects the difficulty of constraining DNRA isotope effects through Bayesian inference, due to its low environmental relevance in the top 1 cm of Lake Lucerne sediments. Although not proven so far, this isotope offset implies that ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reducers impart distinct isotopic fractionation depending on the pathway, which is rather implausible. However, if true, increasing DNRA activity would lead to a stronger ¹⁵N enrichment in the residual ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ pool (Fig. S6d in the Supplement), with downstream impacts on the product pools (N₂ and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$) (Fig. 9c and d).

  Denitrification-derived N₂ mixes with a large ambient N₂ pool (717 µM; δ¹⁵N ∼ 0 ‰), resulting in slightly elevated δ¹⁵N-N₂ values in the top 1 cm. While this increase is subtle (Δδ¹⁵N < 0.1 ‰), it becomes more pronounced as a larger fraction of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ (and subsequently ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$) is reduced to N₂ (denitrification) rather than to ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ (DNRA) (Fig. 9c) due to the distinct isotope effects associated with ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction via denitrification and DNRA. Under full expression of the denitrification isotope effect (i.e., ε_Den1≈ 20 ‰), δ¹⁵N-N₂ much lower than 0 ‰ would be expected; in contrast, ε_Den1≈ 2.8 ‰ likely suppresses such isotopic dynamics, resulting in only subtle δ¹⁵N-N₂ changes. As more ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ is reduced via DNRA (ε_DNRA1≈ 20.0 ‰) than via denitrification (ε_Den1≈ 2.8 ‰), a stronger ¹⁵N depletion is expected in the ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ pool; if this ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ is then reduced to N₂ will lead to lower δ¹⁵N-N₂ than in a purely-denitrifying case. Such interaction can explain the shift toward lower δ¹⁵N-N₂ values as ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ is increasingly reduced via DNRA with a strong isotope effect recorded in our model. Thus, the slightly elevated δ¹⁵N-N₂ values observed in our model confirms that denitrification dominates over DNRA, and operates with a reduced isotope effect (2.8 ‰), likely due to diffusive limitation.

  In contrast, enhanced DNRA activity leads to ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ accumulation and a progressive decrease in δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ in the upper 0.5 cm, consistent with strong isotopic fractionation during DNRA (Fig. 9d). This ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool appears to promote nitrification, as indicated by higher ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation rates (Fig. S6a and b), resulting in the production of ¹⁵N-depleted ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ (Fig. S6c). Notably, if this isotopically light ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ is subsequently reduced via denitrification, it can lead to the formation of N₂ with unusually low δ¹⁵N values, even if denitrification itself operates with a modest isotope effect. This secondary effect underscores how DNRA not only alters substrate availability but also indirectly influences the isotopic composition of denitrification end products. The strong spatial overlap of DNRA, denitrification and nitrification highlights the central role of DNRA in fuelling internal N recycling (Wang et al., 2020) with implications that extend to the δ¹⁵N of both intermediate and terminal N pools.

  Thus, if ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction via DNRA and denitrification occurs with distinct isotope effects, our model has the potential to disentangle their respective contributions based on δ¹⁵N profiles of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, and to a lesser extent of N₂ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$. Importantly, our results underscore a potentially critical, yet underappreciated, coupling between DNRA and nitrification in benthic environments. If verified, this interaction, largely invisible in concentration profiles alone, can significantly influence isotopic signatures and must be considered when interpreting sediment N dynamics through an isotope lens.

• ii.

  N removal via denitrification versus anammox. The results for this case scenario reveal, somewhat unexpectedly, some similarities between denitrification and anammox with respect to ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction to N₂ and associated N isotope signatures. The isotope effects associated with denitrification are low (2.8 ‰ for ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ reduction and 7.9 ‰ for ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction), whereas anammox imparts stronger isotopic fractionation (14.4 ‰ for ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction to N₂ and −30.0 ‰ for its oxidation to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$). These values reflect parameter estimations specific to Lake Lucerne's surface sediments (upper 1 cm), where anammox activity is low. Both ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction and N₂ production peak around 0.5 cm depth, with minor differences in the thickness of the active layer due to variations in substrate affinity between modelled processes (Fig. 9e and f). The total rate of ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction to N₂, via either anammox or denitrification, remains consistent across all case scenarios. Nonetheless, slight differences can be observed in some N pools as anammox becomes the dominant fixed-N loss path. Increased anammox activity leads to elevated N₂ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ concentrations (Fig. 9g and h), likely due to the use of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ as a substrate, which mitigates substrate limitation under low ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ availability (i.e., 1.3 $\mathrm{mol}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$ needed to produce 1 mol  N₂ via anammox versus 2 $\mathrm{mol}{{\mathrm{NO}}_{\mathrm{2}}}^{-}$ via denitrification). When anammox prevails, δ¹⁵N-${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ values increase due to the stronger isotope effect associated with ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ reduction via anammox relative to denitrification. This enrichment is partially counterbalanced by the inverse kinetic isotope effect during ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ oxidation to ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ (Brunner et al., 2013), leading to ¹⁵N-enriched ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ below 0.8 cm; notably, this isotopic shift occurs without significant changes in total ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ concentrations (Fig. S6g and h). Lastly, substantial differences emerge in the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool: higher anammox activity correlates with lower ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ concentrations and elevated δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ values throughout most of the sampled depths (Fig. S6e and f). This isotopic enrichment likely overlaps with the effect of nitrification on the ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ pool in the upper 0.3 cm.

  While some differentiation between denitrification and anammox is evident in the isotope signatures of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, the expected contrasts in the ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ and N₂ pools are surprisingly muted. This near-indistinguishability in isotopic outcomes suggests a degree of functional and isotopic redundancy between the two pathways under the modelled conditions. These results highlight the need for further investigation, particularly through refined isotope-based methods (e.g., inclusion of NO_x O-isotopes or clumped nitrate isotopes) and more mechanistic modelling, to distinguish the respective contributions of denitrification and anammox to N removal in sedimentary systems.

5 Conclusions

We developed a comprehensive diagenetic N isotope model that integrates multiple N transformations in benthic environments. The model's complexity requires the use of prior knowledge in addition to the observed data, in order to achieve the most plausible descriptions of the ongoing processes. To address uncertainty in prior knowledge, and to reduce structural errors associated with fixed parameter values, we applied Bayesian inference for a large parameter set (∼60) for data analysis. The computational demands of this approach were met by implementing the model in Julia, with compatibility for automatic differentiation to allow for advanced Markov chain Monte Carlo algorithms needed for Bayesian inference. Despite these optimization efforts to enhance efficiency, inference runs still took 2–3 weeks of computation time (in addition to preceding simulations to reduce burn-in) to achieve sufficiently good convergence of the Markov chains of the posterior parameter distribution. Alongside concentrations and δ¹⁵N values for different N species, the model provides depth profiles of process rates and all fluxes, including their uncertainties. These outputs enable a detailed assessment of the processes shaping N cycling (i.e., concentration profiles) and isotope patterns in sediments.

Application of the developed model to a test dataset from Lake Lucerne successfully reproduced measured profiles of O₂, ${{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}$, ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$, δ¹⁵N-${{\mathrm{NH}}_{\mathrm{4}}}^{+}$, and δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$. The model also produced realistic vertical distributions of conversion rates, revealing clear depth-dependent zonation. Most marginal posterior distributions of estimated parameters were in good agreement with their priors. Yet, strong deviations were observed for the N isotope effect associated with the first step of denitrification, ε_Den1, which was estimated at ∼ 2.8 ± 1.1 ‰, significantly lower than the expected ∼ 20 ‰. These findings were confirmed by additional simulations performed using narrower priors and a fixed ε_Den1 value of 20 ‰, both of which resulted in a substantial deterioration in the model's ability to reproduce δ¹⁵N-${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ profiles. This, in turn, can be taken as indication for a suppressed denitrification ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ isotope effect at the porewater level in Lake Lucerne, potentially due to process coupling via ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$. The model's ability to quantify such interactions, which can be difficult to discern in situ or from field data alone, is a key strength of this stepwise model framework. A manuscript assessing such dynamics across distinct sites is currently being prepared to further corroborate these findings.

Further sensitivity tests highlighted that the model could still achieve good fits to the observational data even when certain individual processes were excluded, demonstrating the critical role of prior knowledge regarding estimated parameters and their associated uncertainties.

Overall, this study presents one of the first comprehensive diagenetic N isotope models that explicitly incorporate multiple N transformation pathways in a stepwise manner and are validated against field measurements. Rather than serving as a purely predictive tool, this model is intended to stimulate scientific discussion on the quantification of N transformations and isotope dynamics in sediments based on observed data. Future developments could focus on improving identifiability through additional, targeted observations, expanding model validation across distinct benthic environments, and incorporating additional isotope tracers, such as δ¹⁸O of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, to further strengthen the model structure and improve its reliability.

Appendix A: Model processes and stoichiometry

The following equations were developed based on the information reported in Table A1.

Table A1 Overview of all modelled N-transformation pathways, including their stoichiometry and governing equations. R denotes the ${}^{\mathrm{15}}\mathrm{N}/{(}^{\mathrm{14}}\mathrm{N}{+}^{\mathrm{15}}\mathrm{N}$) ratio derived from OM. The γ parameter defines the fraction of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ released during OM mineralization for each reaction. Anammox encompasses both the comproportionation of ${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$ to N₂, defined as the main (“m”) reaction, and the production of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ from ${{\mathrm{NO}}_{\mathrm{2}}}^{-}$, defined as the side (“s”) reaction.

$$\begin{array}{c}{\displaystyle }{r}_{\text{MinOx}}=k_{\text{MinOx}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}},\text{MinOx}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r}_{\text{MinAnae}}=k_{\text{MinAnae}}{\displaystyle \frac{K_{{\mathrm{NO}}_{\mathrm{3}},\text{MinAnae}}}{K_{{\mathrm{NO}}_{\mathrm{3}},\text{MinAnae}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{MinAnae}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{MinAnae}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r}_{\text{MinSulfRed}}=k_{\text{MinSulfRed}}{\displaystyle \frac{\left[{{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}\right]}{K_{{\mathrm{SO}}_{\mathrm{4}},\text{MinSulfRed}}+\left[{{\mathrm{SO}}_{\mathrm{4}}}^{\mathrm{2}-}\right]}}{\displaystyle \frac{K_{{\mathrm{NO}}_{\mathrm{3}},\text{MinSulfRed}}}{K_{{\mathrm{NO}}_{\mathrm{3}},\text{MinSulfRed}}+\left[{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{MinSulfRed}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{MinSulfRed}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Anam}}=k_{\text{Anam}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NH}}_{\mathrm{4}},\text{Anam}}+\left[{}^{\mathrm{14}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NO}}_{\mathrm{2}},\text{Anam}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{Anam}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Anam}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Nit1a}}=k_{\text{Nit1}}\left(\mathrm{1}-f_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}\right){\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NH}}_{\mathrm{4}},\text{Nit1}}+\left[{}^{\mathrm{14}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Nit1}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Nit1b}}=k_{\text{Nit1}}{f}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}{\displaystyle \frac{\mathrm{1}}{{\left(K_{{\mathrm{NH}}_{\mathrm{4}},\text{Nit1}}+\left[{}^{\mathrm{14}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NH}}_{\mathrm{4}}}^{+}\right]\right)}^{\mathrm{2}}}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Nit1}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Nit2}}=k_{\text{Nit2}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NO}}_{\mathrm{2}},\text{Nit2}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Nit2}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Den1}}=k_{\text{Den1}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NO}}_{\mathrm{3}},\text{Den1}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den1}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den1}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Den2}}=k_{\text{Den2}}{\displaystyle \frac{\mathrm{1}}{{\left(K_{{\mathrm{NO}}_{\mathrm{2}},\text{Den2}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]\right)}^{\mathrm{2}}}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den2}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den2}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{Den3}}=k_{\text{Den3}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Den3}}+\left[{}^{\mathrm{1414}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]+\left[{}^{\mathrm{1415}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]+\left[{}^{\mathrm{1515}}{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den3}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{Den3}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{DNRA1}}=k_{\text{DNRA1}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NO}}_{\mathrm{3}},\text{DNRA1}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{3}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{DNRA1}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{DNRA1}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{r^{\prime }}_{\text{DNRA2}}=k_{\text{DNRA2}}{\displaystyle \frac{\mathrm{1}}{K_{{\mathrm{NO}}_{\mathrm{2}},\text{DNRA2}}+\left[{}^{\mathrm{14}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]+\left[{}^{\mathrm{15}}{{\mathrm{NO}}_{\mathrm{2}}}^{-}\right]}}{\displaystyle \frac{K_{{\mathrm{O}}_{\mathrm{2}},\text{DNRA2}}}{K_{{\mathrm{O}}_{\mathrm{2}},\text{DNRA2}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{f}_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}=b_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}{\displaystyle \frac{a_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}}{a_{{\mathrm{N}}_{\mathrm{2}}\mathrm{O},\text{Nit1}}+\left[{\mathrm{O}}_{\mathrm{2}}\right]}}\\ {\displaystyle }{k}_{\text{Den2}}=f_{\text{Den2,Den1}}{k}_{\text{Den1}}{k}_{\text{Den3}}=f_{\text{Den3,Den1}}{k}_{\text{Den1}}{k}_{\text{Nit2}}=f_{\text{Nit2,Nit1}}{k}_{\text{Nit1}}\\ {\displaystyle }{k}_{\text{Anam}}=f_{\text{Anam,Den2}}{k}_{\text{Den2}}{k}_{\text{DNRA1}}=f_{\text{DNRA1,Den1}}{k}_{\text{Den1}}{k}_{\text{DNRA2}}=f_{\text{DNRA2,Den2}}{k}_{\text{Den2}}\end{array}$$

Appendix B: Reaction-diffusion model

Nomenclature

t	time [d]
z	depth coordinate within sediment (0 at the sediment surface, d at the lower boundary of the modelled sediment layer) [cm]
d	depth of the modelled sediment layer [cm]
C(z,t)	substance concentration (mass per volume of water) as a function of depth and time
p(z)	porosity of the sediment (water volume divided by sediment volume) as a function of sediment depth
D(z)	diffusivity of the substance in the water as a function of depth (usually constant and equal to the molecular diffusion coefficient; however, bioturbation could be modelled as an increase in diffusivity close to the sediment surface)
r(C)	transformation rate of the substance (mass per volume of water per unit of time)
C0	substance concentration at the sediment surface
Fd	substance flux from deep sediment into the modelled sediment layer at the lower boundary of the modelled sediment layer (mass per unit of total sediment surface and per unit of time)

Partial Differential Equation for Sediment Layer

Mass balance within the sediment layer:

$$p{\displaystyle \frac{\partial C}{\partial t}}-{\displaystyle \frac{\partial }{\partial z}}\left(Dp{\displaystyle \frac{\partial C}{\partial z}}\right)=pr$$

Differential equation for concentration:

$${\displaystyle \frac{\partial C}{\partial t}}={\displaystyle \frac{\mathrm{1}}{p}}{\displaystyle \frac{\partial }{\partial z}}\left(Dp{\displaystyle \frac{\partial C}{\partial z}}\right)+r$$

Diffusion (molecular diffusion corrected for tortuosity, and bioturbation):

$$D={\displaystyle \frac{D_{\text{mol}}}{a_{\text{tort}}{p}^{\mathrm{1}-m_{\text{tort}}}}}+D_{\text{bio}}{e}^{-\frac{z}{d_{\text{bio}}}}$$

Boundary conditions:

$$C(\mathrm{0},t)=C_{\mathrm{0}},D(d,t)p(d,t){\displaystyle \frac{\partial C}{\partial z}}(d,t)=F_{\mathrm{d}}$$

For N compounds with a single N atom, the boundary conditions are calculated from total concentrations, C_tot, and δ¹⁵N as follows:

$$\begin{array}{rl}& r=\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}{\mathrm{1000}}}+\mathrm{1}\right)R_{\text{std}}\\ & C_{{}^{\mathrm{14}}\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+r}}{C}_{\text{tot}}\\ & C_{{}^{\mathrm{15}}\mathrm{N}}={\displaystyle \frac{r}{\mathrm{1}+r}}{C}_{\text{tot}}\end{array}$$

For N compounds with two N atoms, the boundary conditions are calculated from total concentrations, C_tot, and δ¹⁵N as follows (Drury et al., 1987):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }r=\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}\mathrm{N}}{\mathrm{1000}}}+\mathrm{1}\right)R_{\text{std}}\\ {\displaystyle }& {\displaystyle }{C}_{{}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\mathrm{2}r+r^{\mathrm{2}}}}{C}_{\text{tot}}\\ {\displaystyle }& {\displaystyle }{C}_{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}}={\displaystyle \frac{\mathrm{2}r}{\mathrm{1}+\mathrm{2}r+r^{\mathrm{2}}}}{C}_{\text{tot}}\\ {\displaystyle }& {\displaystyle }{C}_{{}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{15}}\mathrm{N}}={\displaystyle \frac{r^{\mathrm{2}}}{\mathrm{1}+\mathrm{2}r+r^{\mathrm{2}}}}{C}_{\text{tot}}\end{array}$$

Appendix C: Prior values for inference

Table C1 Model parameters estimated using Bayesian inference, alongside their prior values and associated uncertainties. The posterior values (estimated mean with their standard deviation) for the base scenario (Sect. 4.1) are also reported. Parameters are grouped into three categories: (A) reaction rates parameters (i.e., defining process kinetics), (B) isotope parameters (i.e., isotope effects for the modelled processes and the N isotopic composition of OM), and (C) parameters used in the one-step denitrification approach (${{\mathrm{NO}}_{\mathrm{3}}}^{-}\to {\mathrm{N}}_{\mathrm{2}}$ instead of ${{\mathrm{NO}}_{\mathrm{3}}}^{-}\to {\mathrm{NO}}_{\mathrm{2}}^{-}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}\to {\mathrm{N}}_{\mathrm{2}}$). Where a wide range of values was reported in the literature, the most relevant value for benthic environments was selected, and the corresponding reference is reported.

Appendix D: Model discretization

We discretize the partial differential equations outlined in Appendix B using the Method of Lines. This approach involves explicit discretization in space, followed by the application of an ODE solver to the resulting system of ODEs.

Spatial discretization

Numerical discretization of sediment layer (n cells, cell expansion factor f):

Visualization:

Cell boundaries ($i=\mathrm{1},\mathrm{\dots },n+\mathrm{1}$):

$$z_i^b=\left\{\begin{array}{ll}\frac{i-\mathrm{1}}{n}d& \text{for }f<\mathrm{1.1}(i=\mathrm{1},\mathrm{\dots },n+\mathrm{1})\\ \frac{f^{\frac{i-\mathrm{1}}{n}}-\mathrm{1}}{f-\mathrm{1}}d& \text{for }f\ge \mathrm{1.1}(i=\mathrm{1},\mathrm{\dots },n+\mathrm{1})\end{array}\right.$$

Cell midpoints $(i=\mathrm{1},\mathrm{\dots },n)$:

$$z_i^m={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(z_i^b+z_{i+\mathrm{1}}^b)$$

Explanation for the cell expansion factor:

The cell size is approximately (the larger n the closer) proportional to

$${\displaystyle \frac{\partial z_i^b}{\partial i}}={\displaystyle \frac{\partial }{\partial i}}\left({\displaystyle \frac{f^{\frac{i-\mathrm{1}}{n}}-\mathrm{1}}{f-\mathrm{1}}}d\right)={\displaystyle \frac{\log \left(f\right)}{f-\mathrm{1}}}{\displaystyle \frac{\mathrm{1}}{n}}{f}^{\frac{i-\mathrm{1}}{n}}d$$

Comparing these cell sizes at the lower and upper boundaries leads to

$${\displaystyle \frac{\frac{\partial z_i^b}{\partial i}{|}_{i=n+\mathrm{1}}}{\frac{\partial z_i^b}{\partial i}{|}_{i=\mathrm{1}}}}=f$$

This expression clarifies the meaning of the cell expansion factor (approximately equal to the ratio of cell size of lowest to uppermost cell).

Discretized Ordinary Differential Equations

Mass balance within sediment layer cells ($i=\mathrm{2},\mathrm{\dots },n-\mathrm{1}$):

$$\begin{array}{rl}& p\left(z_i^{\text{m}}\right){\displaystyle \frac{\partial C}{\partial t}}\left(z_i^{\text{m}}\right)(z_{i+\mathrm{1}}^{\text{b}}-z_i^{\text{b}})\\ & =-p\left(z_i^{\text{b}}\right)D\left(z_i^{\text{b}}\right){\displaystyle \frac{C\left(z_i^{\text{m}}\right)-C\left(z_{i-\mathrm{1}}^{\text{m}}\right)}{z_i^{\text{m}}-z_{i-\mathrm{1}}^{\text{m}}}}+p\left(z_{i+\mathrm{1}}^{\text{b}}\right)D\left(z_{i+\mathrm{1}}^{\text{b}}\right)\\ & \times {\displaystyle \frac{C\left(z_{i+\mathrm{1}}^{\text{m}}\right)-C\left(z_i^{\text{m}}\right)}{z_{i+\mathrm{1}}^{\text{m}}-z_i^{\text{m}}}}+p\left(z_i^{\text{m}}\right)r\left(z_i^{\text{m}}\right)(z_{i+\mathrm{1}}^{\text{b}}-z_i^{\text{b}})\end{array}$$

Differential equation for concentrations at cell midpoints of inner cells ($i=\mathrm{2},\mathrm{\dots },n-\mathrm{1}$):

$$\begin{array}{rl}{\displaystyle \frac{\partial C}{\partial t}}\left(z_i^{\text{m}}\right)=& {\displaystyle \frac{\begin{array}{c}-p\left(z_i^{\text{b}}\right)D\left(z_i^{\text{b}}\right)\frac{C\left(z_i^{\text{m}}\right)-C\left(z_{i-\mathrm{1}}^{\text{m}}\right)}{z_i^{\text{m}}-z_{i-\mathrm{1}}^{\text{m}}}\\ +p\left(z_{i+\mathrm{1}}^{\text{b}}\right)D\left(z_{i+\mathrm{1}}^{\text{b}}\right)\frac{C\left(z_{i+\mathrm{1}}^{\text{m}}\right)-C\left(z_i^{\text{m}}\right)}{z_{i+\mathrm{1}}^{\text{m}}-z_i^{\text{m}}}\end{array}}{p\left(z_i^{\text{m}}\right)(z_{i+\mathrm{1}}^{\text{b}}-z_i^{\text{b}})}}\\ & +r\left(z_i^{\text{m}}\right)\end{array}$$

Boundary conditions:

$$\begin{array}{rl}C\left(z_{\mathrm{1}}^{\text{b}}\right)& =C_{\mathrm{0}},\\ & D(z_{n+\mathrm{1}}^{\text{b}},t)p(z_{n+\mathrm{1}}^{\text{b}},t){\displaystyle \frac{C\left(z_{n+\mathrm{1}}^{\text{b}}\right)-C\left(z_n^{\text{m}}\right)}{z_{n+\mathrm{1}}^{\text{b}}-z_n^{\text{m}}}}\\ & =F_{\mathrm{d}}\\ & \to C\left(z_{n+\mathrm{1}}^{\text{b}}\right)=C\left(z_n^{\text{m}}\right)+F_{\mathrm{d}}{\displaystyle \frac{z_{n+\mathrm{1}}^{\text{b}}-z_n^{\text{m}}}{D(z_{n+\mathrm{1}}^{\text{b}},t)p(z_{n+\mathrm{1}}^{\text{b}},t)}}\end{array}$$

Differential equations for concentrations at cell midpoints of top and bottom cell $(i=\mathrm{1},i=n)$:

$$\begin{array}{c}{\displaystyle \frac{\partial C}{\partial t}}\left(z_{\mathrm{1}}^{\text{m}}\right)={\displaystyle \frac{\begin{array}{rl}{\displaystyle }& {\displaystyle }-p\left(z_{\mathrm{1}}^{\text{b}}\right)D\left(z_{\mathrm{1}}^{\text{b}}\right){\displaystyle \frac{C\left(z_{\mathrm{1}}^{\text{m}}\right)-C\left(z_{\mathrm{1}}^{\text{b}}\right)}{z_{\mathrm{1}}^{\text{m}}-z_{\mathrm{1}}^{\text{b}}}}\\ {\displaystyle }& {\displaystyle }+p\left(z_{\mathrm{2}}^{\text{b}}\right)D\left(z_{\mathrm{2}}^{\text{b}}\right){\displaystyle \frac{C\left(z_{\mathrm{2}}^{\text{m}}\right)-C\left(z_{\mathrm{1}}^{\text{m}}\right)}{z_{\mathrm{2}}^{\text{m}}-z_{\mathrm{1}}^{\text{m}}}}\end{array}}{p\left(z_{\mathrm{1}}^{\text{m}}\right)(z_{\mathrm{2}}^{\text{b}}-z_{\mathrm{1}}^{\text{b}})}}+r\left(z_{\mathrm{1}}^{\text{m}}\right)\\ {\displaystyle }{\displaystyle \frac{\partial C}{\partial t}}\left(z_n^{\text{m}}\right)={\displaystyle }{\displaystyle \frac{\begin{array}{rl}{\displaystyle }& {\displaystyle }-p\left(z_n^{\text{b}}\right)D\left(z_n^{\text{b}}\right){\displaystyle \frac{C\left(z_n^{\text{m}}\right)-C\left(z_{n-\mathrm{1}}^{\text{m}}\right)}{z_n^{\text{m}}-z_{n-\mathrm{1}}^{\text{m}}}}\\ {\displaystyle }& {\displaystyle }+p\left(z_{n+\mathrm{1}}^{\text{b}}\right)D\left(z_{n+\mathrm{1}}^{\text{b}}\right){\displaystyle \frac{C\left(z_{n+\mathrm{1}}^{\text{b}}\right)-C\left(z_n^{\text{m}}\right)}{z_{n+\mathrm{1}}^{\text{b}}-z_n^{\text{m}}}}\end{array}}{p\left(z_n^{\text{m}}\right)(z_{n+\mathrm{1}}^{\text{b}}-z_n^{\text{b}})}}+r\left(z_n^{\text{m}}\right)\end{array}$$

Appendix E: Model implementation

The model was implemented in Julia (Bezanson et al., 2017) (https://julialang.org, last access: 11 July 2024). The implementation is available with open access at https://gitlab.com/p.reichert/Nsediment (last access: 19 January 2025). The version used for this study corresponds to commit 7afecdf1af871e8f8030360d658ec1cf54d20716.

The partial differential equations described in Appendix B were spatially discretized according to the approach outlined in Appendix D. The resulting ordinary differential equations were then numerically solved by the Method of Lines using the package DifferentialEquations.jl (Rackauckas and Nie, 2017). Discretizing the modelled sediment layer into 50 cells, and considering 14 state variables, resulted in a system of 700 ordinary differential equations. The performance of several ODE solvers was compared, resulting in the use of the adaptive order and adaptive time step backward-differencing solver FBDF to account for the stiffness of the ODE system.

Maintaining compatibility with automatic differentiation while allowing flexible parameter selection for inference was a key implementation challenge. This was addressed by using separate arrays for parameter values and names, and by prepending the parameters to be estimated, ensuring a contiguous array of the parameters. To avoid inefficiencies related to the search of parameter names, the association of parameter names to array indices was resolved within the differential equation solver function. This solver, which includes the function to calculate the right-hand side of the differential equation as an internal function, ensures that the index resolution has to be done only once and remains available for all calls of the integrator by the solver. This approach enabled compatibility of our implementation with the automatic differentiation package ForwardDiff.jl (Revels et al., 2016).

Bayesian inference was implemented with both an adaptive Metropolis sampler from the AdaptiveMCMC package (Vihola, 2020) and the Hamiltonian Monte Carlo algorithm from the AdvancedHMC.jl package (Xu et al., 2020).

All model outputs were written to text files and post-processed using R version 4.4.2 (https://www.r-project.org, last access: 31 October 2024).