the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 29 Aug 2025
A comprehensive porewater isotope model for simulating benthic nitrogen cycling: Description, application to lake sediments, and uncertainty analysis
Abstract. 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.
- Preprint
(3279 KB) - Metadata XML
-
Supplement
(2731 KB) - BibTeX
- EndNote
Status: open (until 10 Oct 2025)
- RC1: 'Comment on egusphere-2025-4089', Anonymous Referee #1, 29 Sep 2025 reply
-
RC2: 'Comment on egusphere-2025-4089', Prof. Jay Brandes, 29 Sep 2025
reply
Overall, I find this to be a worthy contribution to the field and suggest publication with minor changes. It does a good job of examining where measurement limitations are important in a complicated, under-determined system. I am not a modeling expert so I will not address the validity of their approach, other than to say that they do an impressive job of examining the multitude of process in sediment N cycling.
Specific comments- Line 61. Replace ‘all’ with ‘present’. This century has seen tremendous advances in the ability to measure N cycle species, and it would be folly to state that these advances will not continue into the future.
Fig 2 and other color plots. -These should be redone with an eye to increase legibility and distinction between parameters. It is quite difficult to distinguish between yellow, light orange and other similar colors, why not use a wider range?, and please consider those who are colorblind! The best practice is to assume that the reader might only have a greyscale printoff of the figure and make sure that your images are legible in grayscale.
The concept of diffusion in sediments influencing the effective isotopic fractionation expressed in sediment (but not the intrinsic isotopic fractionation of denitrifiers) has been discussed widely in the literature, it is not at all a surprise that they find this as a requirement in their model. They may wish to better acknowledge this in their discussion/conclusions.
Citation: https://doi.org/10.5194/egusphere-2025-4089-RC2 -
RC3: 'Comment on egusphere-2025-4089', Anonymous Referee #3, 06 Oct 2025
reply
The authors present a diagenetic N-isotope model, for use in aquatic sediments. The model is fitted to data and used to estimate the magnitudes of various sedimentary processes. The manuscript is well-written, the model is novel, and the sensitivity analysis and model fitting is state of the art.
In contrast to other diagenetic models, this model considers only dissolved nitrogen species, imposing the mineralisation rates not by modelling organic matter, but rather by imposing the maximum rates of the separate processes. While the mineralisation processes comprise oxic mineralisation, denitrification and anoxic processes, the sulphate reduction is modelled separately. It is not clear why the authors have distinguished sulphate reduction from the anaerobic mineralisation (but I guess this is because sulphate was measured). However, one could think of a simpler model where sulphate reduction would be part of the anaerobic mineralisation.
Ignoring organic matter in the model assumes that the mineralisation is only dependent on the availability of oxidants and not on organic matter. The anaerobic mineralisation is the closure term here and it is not limited by any substrate: it has only inhibition components (p.9). Hence, below the layers where oxygen and nitrate are present, anoxic mineralisation will continue at the same rate for all depths, and integrated anaerobic mineralisation will be infinitely large (theoretically). In the model, this is overcome by imposing an ammonium flux at the lower boundary, which effectively represents a *finite* ammonium production by anaerobic mineralisation. This means that the depth of the model is also an important model parameter, and so one cannot simply extend the model domain, and obtain the same results, as one could do in other diagenetic models. This should be mentioned in the model assumptions section.
More seriously, is that, when looking at the description of oxygen dynamics, the reoxidation of anoxic substances other than ammonium and nitrite is ignored. This implicitly assumes that the concentrations of Fe2+, Mn2+, H2S, CH4 are completely removed in the deep parts of the sediment and therefore do not react with oxygen. While such removal processes may occur in certain sediments, it is rare that they completely remove all of these substances. The authors should also list this important (and perhaps unrealistic) assumption in the model assumptions section (on p 10). A suggestion to make the oxygen dynamics more robust would be to explicitly model a lump-sum of anoxic concentrations in the model that are reoxidised with oxygen, and impose a flux of this lump-sum constituent at the lower boundary.
Some minor comments:
The figures were rather difficult to interpret, due to a color scheme that did not provide enough discriminating power. This made it difficult to follow the discussion.
The ammonium deep boundary flux is imposed. How is this flux divided into 14N-NH4 and 15N-NH4?
The model was dynamically run to steady-state. How was steady-state checked?
L 99 sediment reactivity -> organic matter reactivity
L176: manganese , iron -> manganese and iron oxides
L775: [NO3-] instead of 14NO3- + 15NO3- ?
L330 Wagenigen -> Wageningen
L508. A bioturbation coefficient of 1cm2/day seems to be very high.
Table 1. The reaction for anammox produces organic matter; however this is not so for the nitrification reactions, which are also autotrophic.
Table 1a: (1-ksi) -> (1-ksi/1000)
Table C1 shows fractions of NH4 produced based on aerobic mineralisation, denitrification, DNRA and sulphate reduction. This does not seem to be consistent with the text where it is said that this is determined by the organic matter 15N/14N composition. Why not use the stoichiometry of the reactions to estimate the NH4 production from the OM composition?
fNit2 present in table C2 is not in Table A1.
The default parameter values for most parameters can be found in table C1, but not the rates, and the boundary conditions. All parameter values used for the base run should be presented somewhere in the manuscript.
Citation: https://doi.org/10.5194/egusphere-2025-4089-RC3
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,743 | 27 | 11 | 1,781 | 25 | 18 | 9 |
- HTML: 1,743
- PDF: 27
- XML: 11
- Total: 1,781
- Supplement: 25
- BibTeX: 18
- EndNote: 9
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Mazzoli and colleagues present a model of benthic N cycling that includes a description of isotopic signatures. They apply it to a dataset obtained in Lake Lucerne, discuss the dominant processes in the nitrogen cycle, and make the case that the model is broadly applicable.
The combination or measured porewater profiles of N2, NH4 and NO3, together with their isotopic signatures and rates of denitrification, DNRA and anammox with the model is novel. What also sets this paper apart from previous work and many other early diagenetic modeling efforts is that it uses a strong parameter estimation component. They use Bayesian inference to connect prior knowledge with the data, and carry out sensitivity analyses on parameters for which the marginal posterior distributions differed substantially from the prior.
This is a well-written paper, providing a good overview on biogeochemical N models in the introduction. Overall, I really like this paper. The work seems to be technically sound, the analysis related to the model parameters is excellent (see e.g. the first paragraph of the conclusions). I also appreciate the clear demonstration that model fit to the measured profiles is not a strong indicator that the underlying processes are necessarily captured correctly (conclusion line 762ff).
The point that I struggled with most is the very low isotope effect reported for denitrification, which I am a bit skeptical about given this is model-based and lacks direct observational support. However, the model analysis presented in the paper is well thought out, so I wonder about potential assumptions that could lead to such a finding in the model (but potentially not in nature). One such issue might be that the abundance of active cells is not modeled (assumption iii, which is common). This may lead to large variations in cell specific rates and potentially fractionation effects. Are these likely candidates for the “structural limitations” mentioned on Line 398? If so, consider expanding on this in the discussion of your result. For example, rather low fractionation has been reported by Perez-Rodriguez et al. 2017 (https://doi.org/10.1016/j.gca.2017.05.014) and by Kritee et al. 2012 (https://doi.org/10.1016/j.gca.2012.05.020), and Kritee’s figure 1 suggests lower epsilons at low cell specific rates. Nevertheless, the values reported here are much lower, so a broader discussion of existing supporting or contradictory experimental evidence would be welcome (including past work by some of the authors).
Below are a few more comments:
Model formulation
The mineralization of organic matter is modeled depending on intrinsic rate constants for aerobic and anaerobic mineralization reactions, as well as for sulfate reduction and denitrification. Different pathways are modeled depending on the availability of dissolved electron acceptors, as well as inhibition of anaerobic mineralization by O2 and NO3. Notably, the mineralization rate is modeled to not depend on the availability of reactive organic matter.
While this manuscript clearly expands beyond previous work, I suggest to revise the wording on line 77-79 and 765/766. For example, the work by Rooze and Meile also incorporated (some) stepwise process descriptions that emphasize the role of nitrite as intermediate, and uses data (though not porewater profiles) for validation. Also, some existing models deal more comprehensively with the coupling between different elemental cycles (see below), and the applicability of this particular model to other environments really depends on how substantial this coupling is. For example, I did not see how reduced products from anaerobic mineralization reactions are accounted for (other than ammonium). As a consequence, the model will underestimate the use of nitrate and O2 through their reoxidation reaction. It seems that in the particular application for Lake Lucerne, this may not affect the results greatly (compare rates of anaerobic mineralization to those of aerobic mineralization in Figure 3, indicating that there is little anaerobic mineralization). However, these are modeled rates, and not accounting for the use of O2 to oxidize hydrogen sulfide or reduced metals limits the applicability of the model to a broader range of environments.
All state variables considered I believe are solutes, i.e. the code does not track any solids. This shortens the simulation time and makes the MCMC feasible. However, it also implies that bioturbation does not account for solid phase mixing. Thus, it is not surprising that changing Db has largely no effect on the model results (section 4.2, scenario D). As a consequence, I suspect this is largely a result caused by the model structure and unlikely to be true in reality. Furthermore, it is not clear how the value of Dbio was determined. Its current value drops from 1.16e-9 m2/s at z=0 to 4.26e-10 m2/s at 1cm depth. Thus, it essentially only exceeds molecular diffusion in the top cm. However, non-diffusive mixing typically exceeds solute transport by molecular diffusion more substantially, and bioturbation tends to often dominate the movement of solid phases near the sediment-water interface. Thus, this value should be justified for solutes and if this is hard to do, then consider toning down conclusions related to bioturbation. And because there are no solid phases accounted for, the model also ignores the effect of the precipitation of FeS (maybe not a huge pool, but reported present (e.g. Horw Bay; Table 2 in Urban et al. 1999 https://doi.org/10.1016/S0016-7037(98)00306-8), affecting the extent of use of electron acceptors such as nitrate and dissolved oxygen.
The manuscript says that rates of N cycling were measured, but I didn’t see those results. The rate of (total) organic matter mineralization is on the order of 450 uM N/d near the surface, which translate into approximately 3 umol/L/d or 120 nmol/cm3/hr. (Figure 3: I assume these rates are in uM nitrogen/d because different N processes are being compared. Please clarify if these are in mol C). These rate estimates are about 100 times larger than the 1 nmolC/cm3/hr reported in Fiscal et al. 2019 (https://doi.org/10.5194/bg-16-3725-2019) for Lake Lucerne surface sediment. Please clarify the evidence supporting your numbers (measured rates; or fitting profiles using low estimates of transport as discussed above?)
The model description is clear, but I have some questions about the process descriptions
- denitrification and nitrification: these processes are implemented as multi-step processes, with the rates of the steps following the initial step modeled as a fraction of the rate of that first step. This is done “to prevent unrealistic rates” (e.g. L180). But it also predetermines the “sharing” of NO2- as a key intermediate between different processes, and, for example, with the Monod dependency on nitrite, the second step of denitrification could not occur at a rate higher than what is fueled by nitrate. This may not be a huge issue for N dynamics and I consider this to be a reasonable model simplification, but can you discuss the implications in particular for the isotopic signatures?
- the rate expression of reaction 1b (Table 1) follows a Monod dependency on both of the ammonia involved forming the N2O. This contrasts with pretty much all the other reactions (including the O2 dependency in the same reaction), in which the reaction stoichiometry is not accounted for in the rate laws. Is there any experimental evidence that supports this formulation? And maybe more importantly what are the implications of this choice?
- anammox: this process includes not only the NO2 + NH4 —> N2 reaction, but also the production of nitrate from nitrite. If I understand correctly, the parameter fside (Table A1) is therefore representing the 0.3 NO3/NH4 (Line 205). Please clarify; also define [s] and [m] in table A1.
- Table A1: define gamma_Den1, _Den2, _Den3. Please clarify so it is clear why R is 15N/14N, instead of e.g. 15N/(14N+15N).
Model/parameter analysis
The Bayesian inference analysis is well done and very helpful. For example, I was rather skeptical that a, b from the Ji et al. paper can be directly used in the sediment (p.8, 18). However, these parameters are apparently not impacting the results a great deal. This is a great example of the value of assessing the impact of the parameters on the model outcome.
Visualization/presentation of results
Several figures are difficult to read without magnifying them on the screen. For example, in Figure 3 I had a difficult time identifying which line and process belong together. If another color scheme is not feasible, consider putting information identifying the relevant processes into the figure caption.
Figure S1: spell out what the solid (posterior) and dashed (prior) lines represent, and what ess= … means in the titles of each panel
Line 242: in the reactions of the manuscript you refer to Mn2+, but here is it Mn3+.