The authors presented a multi-G 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 (OMEN-SED) that is originally based on two reactive OM classes (2-G 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 transport-reaction equation of OM can be derived. Because the proposed multi-G approach is based on the analytical solution, which does not require solving the transport-reaction 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 2-G 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 model-data fit at some places do not reflect reality and bio-physcial 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 pre-described (e.g. v=0.15 in case studies) or derived by model-data 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 non-zero 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³ mm is not supported by any observation.

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). 

1) Equations (4) and (8) in the Introduction and on page 7 are strictly valid only in non-bioturbated 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?

Signed: Bernard Boudreau