Comment on gmd-2021-4

This paper reports a means by which the OMEN-SED diagenesis program (Hülse et al., 2018, GMD 11), which uses 2-G 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, BG-2020-435).

This paper reports a means by which the OMEN-SED diagenesis program (Hülse et al., 2018, GMD 11), which uses 2-G 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, BG-2020. 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 non-bioturbated 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 non-mixed 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 non-mixed 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, 2503-2516. 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?