Review of “FESDIA (v1.0): Exploring temporal variations of sediment biogeochemistry under the influence of flood events using numerical modelling” by Stanley I. Nmor et al.

--- General comments ---

The authors provided a model (FESDIA) that facilitates examination of the effect of a flood event on early diagenesis based on the published model OMEXDIA. The data-model comparison was carefully made and relaxation to pre-flood profile was examined as an application example. While the validity of the model looks sufficiently examined, I cannot say the manuscript is carefully written, and more clearance may be desired as a model development paper. Below are my comments that I hope can be of some use to the authors.

(1) Model description

It would be helpful if the authors can provide a short overview of model development relative to the previous and already published works/models. Current manuscript referred to these works, but what is exactly new is kind of obscure. Also, brief summary of the capacity/features of previous models from which the current model has been developed will be useful, including programing language, governing equations, and algorithms utilized for numerical solutions (finite difference/volume/element method for equation differencing and Newton iteration if adopted for solution seeking etc.). The current manuscript provides some of above information, but some information is still missing. And it is not 100% clear which part is new to the current model and which part is not new.

(2) Definition of relaxation time

Eq. 22 describes the change rate of species concentration within sediment profile but does not necessarily define deviation from the pre-flood profile. Also, current manuscript lacks figures that directly compare the solute profile development relative to the pre-event profile. I think Figs. 4-8 need to be improved so that profile development relative to the pre-event profile is more visually obvious.

(3) Description of model limitations and future development

Relevant to the above point, but limitations of model should be discussed more. There are several simplifying assumptions in the model but its influences on e.g., model validation, comparison with observation and estimation for relaxation time are not discussed. For instance, the authors assume that burial rate/porosity does not change with the flooding, but it is not discussed whether this assumption is defendable or close to what we observe. Model validation or comparison with the observed data is essentially based on solute profiles, which likely resulted from a good fit of TOC and may be achievable under different assumptions (those that allow changes of burial rate, porosity, bio-mixing and irrigation etc.). If this is the case, the relaxation time is likely quite different under different assumptions. While most of manuscript discussed how relaxation time is calculated under the specific assumption adopted for this study, it is not discussed how the relaxation time is affected by adopted assumptions.

--- Specific comments ---

In model description, it may be better if you say what programing language you are using earlier on (even in abstract).

L134. Three OM fractions? Thought the authors are using two.

L138. Froelich et al. (1979) rather than Froelich (1988) according to Bethke et al. (2011, AJS 311, 183)?

L166. What does ‘a coupled reaction formulation’ mean?

Eqs. 6, 7. It does not make any sense to use/define ‘maximum rates’ when one is not using Monod or Michaelis-Menten type of equation.

Eq. 7. What is the definition of 𝑟𝐻2𝑆𝑜𝑥𝑖𝑑?

Eqs. 8, 9. What is the assumption behind the formulation of these equations? For instance, how do you obtain Eq. 9 for aqueous NH4+ with accounting for adsorbed NH4+? Can you track NH4+ adsorbed onto solid species along with OM and other solids or do you have to simplify that adsorption is depth-independent and/or time-independent? This can be important if solid materials with unoccupied exchange sites are flooded to sediment depocenter in a short period of time. If such case is possible, one would expect a large sink of NH4+ to the exchange sites? Related to this, do you model PO4 adsorption onto Fe hydroxides or you do not have to do this?  

L205. Does porosity ‘decay’?

L209. According to Eq. 10, the authors seem to assume intraphase biodiffusion (Meysman et al., 2005, GCA 69, 3601). The statement here mentioning an interphase biodiffusion is inconsistent with Eq. 10.

Eq. 14. Irrigation term is not found in the governing equation. Is it included as a reaction term?

Section 2.2.5 & Section 2.2.6. More details are desirable as adding grid for implementing a deposition event must be an important addition to the previous modeling framework. For instance, how you define Zpert, e.g., number of grids and their geometry etc.  Some examples, not only schematics may also be useful.

L276. ‘maximum in the spring and minimum in fall and winter’. This line does not make sense to me.

Eq. 21. What is the units of TOC? Also, how do you derive Eq. 21? Is this simply analytical solution of the governing equation? In any case, it would be helpful if the authors can provide the procedure to obtain Eq. 21 somewhere.

L321. ‘were utilizes the R programming language’. Correct English?

Section 2.2.9. It would be helpful if the authors can make a short description of what part of “deSolve package” they used, not only directing the reader to the R-forge webpage. More specifically, how the authors numerically solve the governing equations, apart from “method-of-lines” methods? Use of any finite difference/volume/element method? How is the time-integration of governing equations made (time-implicitly or -explicitly)?

L331. What is a “slow” stationary state?

Eq. 22. Not quite sure this is a legit mathematical expression. phi(tau) < threshold is what I thought is consistent with what the authors described.

L334. “threshold (i.e given by the median over the entire time duration).” Do you mean that the run is finished when difference becomes less than the median value throughout the simulation and then tau is defined as the model time required for this?

L340. I probably do not fully understand the ensemble of simulations here to estimate the uncertainty in tau. What parameter do you randomly re-sampled exactly? Median of the reference run through time? If so, the runs for determining the uncertainty in tau is conducted until (randomly-chosen) prescribed median is crossed? But this does not necessitate re-running of the model as the boundary conditions are not changed?

L383. “a thickness scale of 1 cm to 30 cm in 5 cm increments”. This line does not make sense to me. What exactly did you use for thickness in sensitivity analysis?

Section 2.2.11.2. If the tested values are not too many, it would be better to list exact values you used for sensitivity analysis.

L507. Please specify what “RiOmar” stands for.

L621. above à below?

--- Technical comments ---

Table 1. What does unequal mark on Fe(OH)3 mean? Is this typo? At least notion should be consistent with that in main text.

L139. Eq. 3 à Eq. 2?

L165. Typo in the second line of Eq. 5.

L168. Right parenthesis in the last line of Eq. 6 is missing. 

L212. Where à where?

L214. specify à specified?

L230. i.e à i.e.?

L232. occur à occurs?

L246. Figure. 2 à Fig. 2 or Figure 2?

L318. “method-on-lines” à “method-of-lines”?

L409. dissolved DIC à DIC?

L415. as thus à as follows?

L418. Table. 3 à Table 3

L433. “Solid” should not be superscript

L434. “Solid” should not be superscript

L535. Improve à improved

L581. introduction à introduction of or introducing

L660. thickness. à thickness

L694. design à designed

Review of FESDIA (v1.0): Exploring temporal variations of sediment biogeochemistry under the influence of flood events using numerical modelling by Stanley I. Nmor et al.

The consequences of floods on the chemistry of solids and porewater in sediments of the Rhône River delta have been reported by previous investigators. Here, the authors build up on previous developments of the OMEXDIA model, updated recently by De Borger et al. (2021, https://doi.org/10.5194/gmd-15-2105-2022) to represent the effects of bottom trawling as a process that removes an upper layer of sediment. The present study aims at simulating a conceptually opposite process: (1) a flood that adds an upper layer of sediment. Another new addition is (2) the presence of new metabolic pathways involving S and Fe, and (3) the introduction of a an enrichment factor, denoted alpha, that aims to account for the fact that deposited sediments can be either depleted or enriched in reactive carbon depending on the flood type and on the nature of sediments transported. Besides these three aspects, the model does not seem to introduce any update compared to earlier versions. The text also suffers from a lack of clarity and structure, inconsistency in the names of variables, and between parameter values reported in the text and in tables. Even though the parts about relaxation time and carbon enrichment were novel and interesting, particularly in the context of the Rhône delta flood case study, the manuscript is overall hard to follow. I have a few general comments and a list of specific comments for improvement.

General comments:

• Rassmann et al. (2016 https://doi.org/10.5194/bg-13-5379-2016, 2020 https://doi.org/10.5194/bg-17-13-2020) described Rhône river delta sediments rich in calcium carbonates, and reported signs of the various reactions associated with calcium carbonates occurring in those sediments. If the focus here is Rhône river delta sediments, why not including any calcium carbonate species in the current model?
• Section 2.2.6: this seems to be the most important section in terms of model development, but it is also the part I had the hardest time to follow. There are a lot of new terms introduced here and they are not well defined. The second paragraph starts to explain how post-flood organic carbon contents are derived, then mentions the solutes, then goes back to post-flood organic carbon content with Eq.(15)… I suggest putting more effort clarifying this section, defining terms with precise and consistent words, in order not to confuse the reader. Specifically, what is the carbon enrichment factor (confac) exactly, and how does it differ from the proportionality constant (pfast)? What is Cflood/org and how does it differ from TOC (both are present in Eq. (15)? Is confac tuned for each simulation or is it constant? Is pfast tuned for each simulation or is it constant?
• It is stated L187-188 that “For dynamic simulation, w can change as a function of time but in most cases we assumed a constant value.” In which cases exactly was w changing? Changing w in all cases seem like a necessity given that the novelty of the model is to simulate events in which the flux of deposited material (thus w) is strongly changing with time. How can a constant w be appropriate to simulate a flood? w also changes with sediment depth, because of chemical reactions occurring within the sediment (see Munhoven, 2021 https://doi.org/10.5194/gmd-14-3603-2021). Can the authors either better justify their choice of a non-changing w or update that in the model simulations?
• Section 2.2.7: in most O2 and pH microprofiles from the Rhône delta presented in Rassmann et al. (2016) we can see the influence of a diffusive boundary layer. Please discuss and justify the absence of diffusive boundary layer control on solutes as an upper boundary condition, or update the upper boundary condition accordingly to include this, as other models do in a simple manner (Boudreau et al., 1996 https://doi.org/10.1016/0098-3004(95)00115-8; Munhoven, 2021 https://doi.org/10.5194/gmd-14-3603-2021; Sulpis et al., 2022 https://doi.org/10.5194/gmd-15-2105-2022).
• There are a lot of inconsistencies between number reported in the text and those in the tables (w, NC ratio, rslow, bottom boundary conditions). Please update and be consistent.

Specific comments:

Shouldn’t “Rhone” be spelled “Rhône”, even in English language?

Abstract

L21-24: Here the enrichment factor alpha is mentioned but not clearly defined. This is confusing. Please update.

Introduction

L35: The use of the acronym RiOmar is not really needed, since only used once after. In general, avoid unnecessary acronyms.

L36: Although more commonly used, POC is also an unnecessary acronym here, since only used once after.

L36-39: The sentence is unclear. “because it serves as a sink for particulate organic carbon and nutrients as well as an intense site of carbon and nutrient”: what is the “it” referring to?

L40: I am not convinced that all the cited models have time-dependent capabilities, unlike several other, more recent models, published in this journal that explicitly do. Please update the list.

L43: “massive episodic events” could refer to lots of processes, please be more specific.

L47-50: Sentence unclear. “Attempts to use mathematical models to understand perturbation-induced events on early diagenetic processes have resulted in a variety of approaches that incorporate this type of local phenomenon.”: what is “this type of local phenomenon” referring to?

L48-50: “As an example, previous research in deep-sea systems suggests that megafaunal perturbation can cause a 35% increase in silicic flux when compared to steady-state estimates (Rabouille and Gaillard, 1990)” this is interesting but this level of precision seems unnecessary, what is the relevance for this study? Besides, what is a “silicic flux”? In which direction is the mentioned flux going?

L50: What is the “redox boundary”?

L52-53: What does the “redistribution of solid-phase manganese with multiple peaks” mean?

L62: “porewater species like oxygen (O2) can be restored after a few months”: it is unclear. Do you mean that porewater concentrations can be restored to their pre-flood levels?

L66: what does “short-lived species” mean?

L67: DIC is a component, not a species.

Materials and methods

L93: “the organic matter delivered reflects the Rhone River inputs (Lansard et al., 2008; Cathalot et al., 2013)”, in terms of what? Composition? Reactivity?

Fig.1: I assume that the dashed and solid gray linings shown on the map depict bathymetry; it would be useful to precise it in a caption/legend

L107: what does “mode of behavior” mean?

L140: how exactly do “the reactivities decrease with depth” in the present model? From Table S1, it seems that the reactivities are constant.

L140: The sentence formulation is awkward: it is the degradation that would “cease”, not its rate. Saying this also slightly exaggerated, degradation rates become indeed very small deep below the sediment-water interface but they are never really equal to zero (e.g. Bradley et al., 2020 https://doi.org/10.1126/sciadv.aba0697).

Eq. (7): What are FeSpro and H2Soxid? Are they different from the FeS and the H2S produced by the reactions shown in Eq. (5)? Why is one rate a capital R and the other a lower case r?

Eq. (9): What is the value of kads and can you give some information on this aspect of the model?

L216-220: How is irrigation implemented into the model, i.e., where does it appear in Eqs. (8 & 9)? 

L236: What is a “time run”?

Eq. (15) Please precise here that TOC_old­ is the TOC concentration at the old sediment-water interface

L248: it would be good to have more information on confac (alpha): here it is tuned. How should it be used in future applications? Always to the same value? Does its value depends on type and magnitude of flood?

Fig.2: change “reactive Corg” for the notation “C (superscript)fast (subscript)org” for consistency

Table 2: the value for rslow is 0.0 d^-1, but in the text it is indicated as 0.0031 d^-1. Please clarify that

Section 2.2.7: what about bottom boundary conditions? Is the concentration really set to 0 for all species, as indicated in Table S1? That would seem unjustified.

L289: By sedimentation rate do you mean solid burial velocity? Porewater burial velocity? Both should be different because porosity is not constant with depth.

L289: Is w 0.027 or 0.03 cm per day? Be consistent between the text and tables.

L283: Why a different NC ratio for both organic matter fractions? How were the values of 0.14 and 0.1 obtained? Why are these values different from those shown in Table 2?

L306: Can you provide details (i.e., show the formula) on how are equations 8-10 integrated?

L306: Please provide guidance on what dt values should users set depending on the simulation

L313: First time the “mix” perturbation is mentioned. What is that?

Eq.(22) Why is it summed over the total number of grid points? Any perturbation following a flood should be the highest near the sediment-water interface, so wouldn’t using data coming from deeper in the sediment to compute the relaxation time dilute the true signal and induce additional uncertainties?

Results

Fig4: what is the alpha value for the slow organic carbon fraction?

L409-412: Can the authors interpret the mismatch between modelled and observed SO4, DIC and NH4 values at depth? Wouldn’t that argue for overestimated organic carbon reactivities at depth?

Table 3: How can there be an oxygen flux to the sediment that is ten times smaller than a DIC flux from the sediment? Wouldn’t a value closer to one be expected?

Fig.9: what is “degradable OM”? does that mean that the alpha value is the same for both fast- and slow-decay organic carbon? If so, precise it.

L524: “mixing events” are again mentioned as something the model is able to simulate, but they are not described earlier, so it is unclear what they are.

Section 4.4: to add to this discussion, and in reference to the mention earlier in the manuscript of a “perturbed trajectory frequently arbitrarily divided into a fast, transient phase and a slow, asymptotic stage”: should we instead think about relaxation time as the time necessary for most of, rather than all, changes to occur, similar to the concept of half-life in radioactivity?

Reference

[Ait Ballagh et al., 2021] is missing from the list