Submitted as: model description paper 08 Mar 2021
Submitted as: model description paper  08 Mar 2021
OMENSED(RCM) (v1.1): A pseudo reactive continuum representation of organic matter degradation dynamics for OMENSED
 ^{1}BGeosys, Department Geoscience, Environment & Society (DGES), Université Libre de Bruxelles, Brussels, Belgium
 ^{2}Department of Earth Sciences, University of California, Riverside, CA 92521, USA
 ^{1}BGeosys, Department Geoscience, Environment & Society (DGES), Université Libre de Bruxelles, Brussels, Belgium
 ^{2}Department of Earth Sciences, University of California, Riverside, CA 92521, USA
Abstract. OMENSED is a onedimensional, analytical reactiontransport model for early diagenesis in marine sediments. It explicitly resolves organic matter (OM) degradation and associated biogeochemical terminal electron acceptor, reduced species and nutrient dynamics in porous media under steady state conditions. OMENSED has been specifically designed for the coupling to global Earth System Models and the analytical solution of the coupled set of mass conservation equations ensures the computational efficiency required for such a coupling. To find an analytic solution, OMENSED expresses all explicitly resolved biogeochemical processes as a function of OM degradation. The original version of OMENSED contains a relatively simple description of OM degradation based on two reactive OM classes, a socalled 2G model. However, such a simplified approach does not fully account for the widely observed continuous decrease in organic matter reactivity with burial depth/time. The reactive continuum model that accounts for the continuous distribution of organic compounds over the reactive spectrum represents an alternative and more realistic description, but cannot be easily incorporated within the general OMENSED framework. Here, we extend the diagenetic framework of OMENSED with a multiG approximation of the reactive continuum model (RCM) of organic matter degradation by using a finite, but large number of OM fractions, each characterized by a distinct reactivity. The RCM and its multiG approximation are fully constrained by only two free parameters, α and ν, that control the initial distribution of OM compounds over the reactivity spectrum. The new model is not only able to reproduce observed pore water profiles, sedimentwater interface fluxes and redox zonation across a wide range of depositional environments, but also provides a more realistic description of anaerobic degradation pathways. The added functionality extends the applicability of OMENSED to a broader range of environments and timescales, while requiring fewer parameters to simulate a wider spectrum of OM reactivities.
 Preprint
(3955 KB) 
Supplement
(175 KB)  BibTeX
 EndNote
Philip Pika et al.
Status: final response (author comments only)

RC1: 'Comment on gmd20214', Anonymous Referee #1, 06 May 2021
The authors presented a multiG approximation method, in which organic matter (OM) can be divided into a large number of classes with different degradation rate to represent a peudo reactive continuum of OM, as an extension / further development of an existing analytical early diagenetic model (OMENSED) that is originally based on two reactive OM classes (2G model). The approximation is based on two assumptions, namely (1) the distribution of reactivity/degradation rate of OM in marine sediments can be reasonably described by a gamma function, and (2) the vertical OM distribution in sediments is in an equilibrium status (i.e. the temporal gradient of OM is zero at any depth in sediments) so that an analytical solution of the transportreaction equation of OM can be derived. Because the proposed multiG approach is based on the analytical solution, which does not require solving the transportreaction equation dynamically, therefore computational expense is not a hindering factor. This makes the proposed method appealing.
However, I have several major concerns: 1) although the method seems to produce reasonable results, the authors did not provide convincing arguments that the proposed method outperforms the original 2G model; 2) the method is described in an unclear manner; 3) while the method trys to fit observational vertical profiles in sediments, the boundary conditions needed for the modeldata fit at some places do not reflect reality and biophyscial laws; and 4) the precondition for the validity of the approach, namely a zero temporal gradient of OM at any depth in sediments, can hardly be met in a dynamic environment. This makes the approach of limited use for coupling to dynamic models in which sedimentation of OM is variable, which is especially true for continental margins.
Specific comments:
* The method (section 2) is not described clearly.
(a) From eq.6, it is stated that om(k,t) represents the probability density function that determines the amount of bulk OM with a reactivity between k and k+dk at time t. As om (k,t) is a probability density function, the sum of om (k,t) across all k at any speficic t should always be 1. However, this is not satisfied in eq.7, in which om (k,0) is dependent on OM(0). Please clarify this.
(b) I am confused by the definition and use of k. k is supposed to indicate the reactivity of OM, which is a variable. So what is the justification of eq.8 that k is determined by a, v and sediment age? I understand that the latter three parameters at any speficic depth are either predescribed (e.g. v=0.15 in case studies) or derived by modeldata fitting. This means that k is also fixed by these values, which is not variable any more. Further, how is the age of sediment layer at depth z derived? It seems that this quantity is another variable which needs to be solved in the method, in addition to a and v. This contradicts the statement and conclusion that only a and v need to be solved. Please clarify.
(c) Another parameterization of k using eq. 16 clearly violates the original relationship as mentiond in (b). Please justify the validity of the method if a different parameterization of k is used.
* In the case study 3.1, although the free varibales a, v are tuned that the model produces results close to observed sediment profiles, their setting has no mechanistic connection with other environmental variables, e.g. in Tabl 1, it is not clear why z_{bio} is set to 0.01 cm at depth 585m, which means that there is no bioturbation at all, but then why D_{bio} has a nonzero value and how these parameters are related to the setting of a and v? Also it is not clear why a has a very small value (corresponding to very small lifetime of OM, therefore quite labile component) for depth 2213 m. Compared to a very confusing setting in this case study, the setting in 3.2 (Table 2) seems more reasonable and respects reality.
* There is hardly justification for the validity of the approach in global application in section 3.3, as shown in Fig.6. In particular, the part that simulated OPD exceeds 10^{3} mm is not supported by any observation.

AC1: 'Reply on RC1', Philip Pika, 08 Jul 2021
The comment was uploaded in the form of a supplement: https://gmd.copernicus.org/preprints/gmd20214/gmd20214AC1supplement.pdf

AC1: 'Reply on RC1', Philip Pika, 08 Jul 2021

RC2: 'Comment on gmd20214', Bernard Boudreau, 19 May 2021
This paper reports a means by which the OMENSED diagenesis program (Hülse et al., 2018, GMD 11), which uses 2G organic matter (OM) decay kinetics, can be altered to approximate a reactive continuum decay model, instead. The interest in this modification lies in better approximating the observed decrease in OM reactivity with depth in sediments, which is a desirable goal. I do have four comments on the execution of this paper, one of which reflects one that I made on a recent related submission by some of these authors, i.e., Freitas et al. (2020, BG2020435).
1) Equations (4) and (8) in the Introduction and on page 7 are strictly valid only in nonbioturbated sediments, where there is a simple relationship between time and depth, a point that is only briefly mentioned in the present manuscript. These equations are stated and applied without any serious caveat. The correct expression with bioturbation is given in Boudreau and Ruddick (1991), their Equation (49), and it is not in an easily implemented form, as needed here.
No one has adequately explored the effects of applying Equation (8) to a bioturbated sediment. In Table 1 of Boudreau and Ruddick (1991), there are 6 applications of Equation (8) to nonbioturbated or negligibly bioturbated sediments. These generate ð values between 0.05 and 0.3; the Peru Margin sediments in that same table are almost certainly bioturbated and have associated ð values closer to unity (0.8 and 0.91). This dichotomy of results may indicate the effects of improperly dealing with bioturbation, but I really cannot affirm this. I also recognize that the actual calculations are made with a discrete numerical approximation, which removes much of the mathematical challenge, but also assumes a ð characteristic of a nonmixed sediment; this may create unintentional modifications to the calculated “a” value.
Further to this point, om(k,0) in Equation (7) is said to be the initial distribution of reactivities. In a nonmixed sediment, that distribution is exclusively altered by the decay of components with time. In a bioturbated sediment, “older” and “newer” organic matter components are mixed together, so that in leaving the mixed layer, the organic matter reactivity distribution is different than om(k,0) at the SWI. Thus, the mathematics get even more demanding.
Progress in global modelling is needed to meet the challenges of global change; I appreciate that; corners will be cut by necessity. However, explicitly recognizing and stating those shortcuts is not only desirable, bit is essential.
2) The methodology of this paper and of Hülse et al. (2018) is to create a diagenetic module that can be solved analytically and to calculate the resulting integration constants by fast linear matrix methods (see section 2.3.1 of Hülse et al., 2018). I am confused about this second part, as there is insufficient detail in either paper to see how the latter can be accomplished. Take sulfate in its reduction zone. The analytical solution will be a sum of exponentials, two of which will be multiplied by unknown integration constants. Boundary conditions at the top of the zone and at the base of the zone will allow calculation of the two constants, and these resulting equations will indeed be linear in the two integration constants. However, the position of the base of the sulfate reduction zone is a further unknown, a point that was extensively discussed in Boudreau and Westrich (1984, GCA 48, 25032516). A third boundary condition is needed to determine the unknown z_{SO4}.
Hülse et al. (2018) do supply a third equation, and so it should be in the present paper too. What puzzles me is that this third equation, as in Boudreau and Westrich (1984), is a nonlinear function of z_{SO4}. Thus, linear matrix methods cannot be employed to arrive at a solution. What have the authors done to avoid this problem, which seems fundamental to their paper?
3) I fail to agree that the model predictions and observations provided in Figure 6 “generally agree well with observations and capture also the widely observed increase in OPD with water depths (see Fig. 6).” The model in this paper generally overpredicts the OPD, more often by far more than an order of magnitude. The data appear to have an upper limit of 1000 mm with ocean depth, which is not present in the model predictions. The obvious discrepancies are not discussed in the paper. They really need to be. Is the “observed” 1000 mm limit an artifact or real? Why the common overprediction of the OPD? There is something valuable to be learned from negative results!
4) Finally, as much as Equation (15) might prove valuable, it is repeated as Equation (17) with a discussion that is similar. Is this duplication needed or am I missing something?
Signed: Bernard Boudreau
 AC2: 'Reply on RC2', Philip Pika, 08 Jul 2021
Philip Pika et al.
Model code and software
OMENSED(RCM) (v1.1): A pseudo reactive continuum representation of organic matter degradation dynamics for OMENSED Philip Pika, Dominik Hülse, and Sandra Arndt https://doi.org/10.5281/zenodo.4421777
Philip Pika et al.
Viewed
HTML  XML  Total  Supplement  BibTeX  EndNote  

253  79  11  343  15  1  3 
 HTML: 253
 PDF: 79
 XML: 11
 Total: 343
 Supplement: 15
 BibTeX: 1
 EndNote: 3
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1