Dear reviewers,

We are very thankful for the time invested in reviewing our manuscript and the constructive comments and suggestions. After considering all the comments, we plan to make the following major changes to the manuscript:

1) Add a section with a theoretical derivation for reaction terms before the section with the applications.

2) Add a subsection to compare the formulation of reaction terms for central and non-central moments.

3) Add a table that lists all parameters and explains their meaning.

4) Provide a technical description of the numerical models used for testing the applications in a supplement.

5) Add validation for derived diffusion equations by applying a delta distribution in the supplement.

6) Improve the introduction

7) Elaborate on the use of different distributions in numerical simulations

We will refer to these points in the response and explain why we plan to make these changes. We will also make other minor changes to correct mistakes or improve the presentation.

Response to reviewer 1

Original comment: The authors presented formulation of diffusive transfers of materials with central moments in detail. In its application to OM early diagenesis the authors also introduced production/consumption of central moments focusing on age and OM reactivity though in less detail and less organized way in my opinion. Overall, new approaches are very intriguing in that they may require less assumptions on the continuum and have a potential of more flexibility to track variable properties of interests in materials important for early diagenesis in marine sediments including age and reactivity. My concern is that I had some difficulty in following the formulation of production/consumption of moments and comparison of the features of the new approach to those from others adopting noncentral moments might be insufficient. If the difficult formulation of the production terms is actually a feature of central moments compared to noncentral moments, I would like the authors to provide more detailed explanations/formulations/comparison. Otherwise, I think this paper is suited for GMD.

Answer: We are glad to read that the reviewer finds the new approaches intriguing and thank him/her for his/her constructive remarks. Regarding the formulation of production/consumption terms, he/she raises the concern that a) the derivations are hard to follow and b) that the formulation of reaction/production terms for central moments and non-central moments are not compared and discussed.

To address point ‘a’: we will reorganize the manuscript and add a more formal derivation for the reaction terms in a new section before the applications (see point 1). Many derivations in the section for the applications can then be moved to this section, and this will help to streamline the presentation of equations.

Regarding point ‘b’: The formulation depends on the process. When new material is produced with age = 0, the product of non-central moments and concentration is not affected (i.e., dCμ/dt=0), whereas the central moments are affected, leading to a more complex formulation. For consumption reactions that do not discriminate with respect to age/reactivity, both non-central and central moments are unaffected, and both require a term to account for the changing concentration. When the reaction does discriminate, the formulation for non-central moments is easier. However, the formulation for non-central moments can be succinct and systematic (i.e., automatic) by using a formulation with binomial coefficients. More importantly, in practice it will probably not matter for numerical simulations, as the system of equation for non-central and central moments will depend on the same integral evaluations. Finally, aging does affect the non-central moments (requiring additional terms) but does not affect the central moments. We will elaborate on this in the revised manuscript (see point 2).

Original comment: Definition of variables is sometimes confusing. Multiple uses of X, µ, w, k, j for different variables, for instance.

Answer: In the new manuscript, we will add a table that lists the variables and explains their meaning (see point 3). We will ensure the correct, unambiguous usage of these variables.

Original comment: L163-166. Production of total age in Eq. 25b makes sense to me but is there any way to verify this formulation?

Answer: In the new manuscript, we will give a more formal derivation for reaction terms in a new section (see points 1, 2) and move these derivations from the applications section to this new section.

Original comment: L169. I assume tau is defined as age but did the authors already define the symbol before? If this is the case, then is C necessary before tau_i on the right-hand side of Eq. 26?

Answer: This variable should have been X_i. We apologize for this inconsistency and will improve it in the revised manuscript.

Original comment: Section 3.1. More explanation would be helpful in general. Particularly, derivation of production terms in Eqs. 25b and c. For instance, some examples with specified q value would be helpful, or some more explanation with example when switching tracer from age to reactivity in the next section (Section 3.2)? Also, it may be interesting to have some insight into the difference in formulation of production/consumption terms from previous studies especially those adopting noncentral moments (e.g., Delhez and Deleersnijder, 2002). I assume the central moments may have more complicated formulation for production/consumption compared to non-central counterparts in general? If this is actually a general feature of central moments, this should be described.

Answer: We agree that the presentation should be improved. Please, see the previous reply about putting the formal derivations for reaction terms in a new section (see points 1 and 2). The derivations will be more general but can be directly implemented in the code. Writing binomial coefficients for specific examples (e.g., equation 32) out by hand would take too much space. Therefore, we do not intend to include these steps in the manuscript.

Delhez and Deleersnijder only considered aging but not reaction processes that depend on the distributions of the property of interest (e.g., age). As explained above, we will compare the formulation for central to non-central moments.

Original comment: L187. Need more explanation/description on “particle-based simulation”.

Answer: The particle-based simulation is a simple toy model that is of little interest to the general audience, but we agree it should be described. We will include a technical description of the model in a supplement (see point 4).

Original comment: L187-190. Difference between two experiments was attributed to step size in the Lagrangian simulation. This seems to be easy to check and I think the authors should provide some results on this. Also, if adding production/consumption terms are possible with the code as mentioned by the authors, why not checking the formulation based on those numerical experiments?

Answer: Properly evaluating the effect of changing the step size in the Lagrangian model is not so straightforward since D also depends on the step size. However, we have been very confident in the analytical derivations of the diffusion equations because we have also verified these analytically with a delta distribution. Therefore, the slight mismatch is not due to an error in the Eulerian equations. The proof based on the delta distributions is rather lengthy, but we will include it in the supplement (see point 5).

Original comment: L197-198. Eq. 30b and 30c may want to be explained in more detail with using more generalized equations (something similar to those in Section 3.1 but more generalized as done for diffusion)?

Answer: For these equations, we did not give a derivation. We will explain these equations in more detail in the new version of the manuscript (see points 1 and 2).

Original comment: L222. I thought you do not need lower boundary conditions assuming bioturbation is limited within top ~10 cm.

Answer: A model must always have defined upper and lower boundary conditions. However, the reviewer is correct that since the bioturbation will have dropped to 0, the zero gradient condition at the bottom will not matter for diffusion. For completeness, it is still good to mention the specified condition.

Original comment: L223. More details on 30-G model would be helpful.

Answer: We will add that 30 evenly spaced bins for reactivity represented the entire reactivity spectrum of the deposited organic matter.

Original comment: L232-236. I think the explanation of production of moments is less organized in Section 3 than that for diffusion in Section 2. Could it be possible to make explanations/formulations easier to follow? This might include, e.g., giving more general formulation first and then presenting different cases later (like for diffusion, first Section 2 and then Table 1 and Appendix), giving the final forms of those production terms after the manipulation of equations, and so on.

Answer: We will address this by combining the derivations for applications 2 and 3 in a new section (see points 1 and 2). Although these formulations contain binomial coefficients, they can be immediately included in a script, and it is unnecessary to write them out.

Original comment: L244. As no values of w1-w3 are given, I have no clue what I am supposed to see in Fig. 4A.

Answer: The point was to compare the “real” distribution in black to a reconstruction based on a direct least-squares fit and an indirect fit based on the moments. The revised manuscript will show the parameters used for all examples, and we will make the figure caption more clear.

Original comment: L250. How are the weights determined? Assume they are determined rather arbitrarily after reading Discussion?

Answer: Considering the definition phi = (X-mu)^q, one can expect eps_2 to be of magnitude (eps_1)^2 and eps_3 to be of magnitude (ep2_2)^(3/2). eps_2 and eps_3 are thus scaled to similar magnitude. For eps_1, we divided through a smaller number to create a higher weight artificially. We may change the method in the revised manuscript as we found an easier way to optimize the parameters.

Original comment: L261. More details on “discrete simulations” would be helpful. Confusing because even continuum model does numerically integrate fitted function (L254-256) if I am not mistaken?

Answer: Yes, that is correct. We will explain the discrete simulations in more detail and avoid confusion.

Original comment: L295. Does this difficulty come from numerical diffusion and/or grid size, i.e., numerical errors? Do numerical errors get bigger for higher order of moments?

Answer: For this example, diffusion and grid size do not play a role. In Figure 4, we directly fit a hypothetical distribution (which was produced by interpolating a “discrete” distribution from an arbitrary simulation), treating it as the actual distribution. Hence, how this distribution was formed does not matter. The question is how well it can be reproduced based on the moments (mean, variance, and skewness). The point is that the extremes have much higher weight, as (X_i-mu) will be larger, which will scale exponentially with the order of moment due to the definition phi_q = (X_i-mu)^q. This is problematic since the extremes may not be that important for the overall reactivity. For instance, a tiny amount of extremely old material will only slightly lower the overall reactivity but greatly impact the higher moments. Therefore, simulating a higher number of moments may be unpractical.

Original comment: L301. Or better/more stable solution seeking like with e.g., ensemble Kalman filter?

Answer: This is an approach we have not tested, but that could work.

Original comment: Appendix A. More detailed derivations of Eqs. A10, A11, A12 would be desirable.

Answer: We will add steps to these derivations.

Original comment: L18, 43. “Kuderer et al.” should be replaced with “Kuderer”.

Answer: We will corrrect the mistake.

Response to reviewer 2

Original comment: I have read through the manuscript written by Rooze et al. This manuscript derived a new diffusion-advection-reaction type of partial differential equations based on centralized moment. Compared with Lagrangian frameworks, this Eulerian PDEs can be analytically evaluated and are computationally less expensive. The authors applied this PDEs to simulate organic matter age and reactivity with mixing processing in marine sediments. Overall, this work gives an opportunity to include bioturbation in the early diagenesis model with continuous distributions instead of the Multi-G model. However, the manuscript is not well organized and is difficult to follow, and it needs extensive revision to be accepted for publication.

Answer: We thank the reviewer for his suggestions and comments. The organization of the manuscript will be improved. The main changes will involve adding more formal derivations of production/consumption terms in a new section (see points 1, 2) and removing the corresponding derivations from the section about the applications. Then these derivations will be less scattered throughout the text. We will add a table of the variables with explanations and ensure they are used more consistently.

Original comment: In the introduction, the authors can emphasize the necessity of the continuous distributions with mixing processes, what is the current state of this model? Multi-G models have often been used to simulate the organic matter degradation with mixing processes including bioturbation. The challenges of such continuous model for bioturbation can be reviewed in the introduction.

Answer: Multi-G is a discrete approach often used to simulate mineralization in bioturbated sediments. The drawback of the multi-G approach is the somewhat arbitrary assignment of reactivity classes, which makes comparing different model parameterizations harder/impossible. Eulerian models with continuous OM reactivity description have not accounted for the effect of bioturbation. Freitas et al. (2021) recently compiled parameters used for the continuous gamma distribution in various study sites. However, for strongly bioturbated environments, they applied a discrete multi-G model. This inconsistency highlights a clear problem and poses several questions. As bioturbation typically dominates transport in the upper centimeters of sediment, the approach may be mechanistically flawed, reducing its predictive power in most environments.

We will explain this in the introduction (see point 6). The drawback of the multi-G approach was mentioned in the discussion: “This approach could be attractive to replace the multi-G approach, as it does not require the somewhat arbitrary definition of various reactivity classes (Jørgensen, 1978)” (lines 282-283), but this will also be moved to the introduction.

Original comment: In whole section 2 and Appendix A, all equations were derived without any references except Crank (1956). I don't know which equations are newly derived and which are from literatures. Please check and explain them carefully.

Answer: When we provide derivations without citation, it implies we derived the equations ourselves. When we obtain an equation that other scientists obtained, we specifically mention that (see end of sect. 2.2). For some very general equations, such as basic calculus, the equation for chemical (Fick’s) diffusion, and the conversion of central to zero moments, we did not add a reference.

Original comment: In this manuscript, the authors have used confusing definition. For example, tau and mu are defined as particle age and averaged age (L68 - L69), but later tau=k is defined as reactivity (L192). k=f(tau) was also used (L227). tau was not even defined in int{tau^x C_tau dtau} (L172). Please make all definitions consistency.

Answer: We thank the reviewer for this suggestion. We indeed used tau and X_i by accident for the same purpose. tau=k was used, since this may be more intuitive for chemists, However, we will reconsider this. We will add a table that lists all variables (see point 3) and ensure consistent usage in the new version.

Original comment: I think the age defined in this manuscript is not age, rather transit time. The distribution function in age-structure model was usually expressed as f(tau,k,t) depending on decay reactivity k, age tau and time t. The average reactivity and age, and high-order moment can thus be calculated from f(tau,k,t). At the sediment-water interface, the age of organic matter will not be zero. The reactivity and age will have a distribution (e.g. gamma distribution) at the sediment-water interface. If this distribution is included in the model, could you derive the PDEs for the transit time, age, reactivity based on your model?

Answer: In the context of organic matter reactivity, the word “age” is used. In physical transport models, “transit time” may be more common, but there is not much difference conceptually. We shall mention “transit time” once in the introduction, but we will call it otherwise age for consistency.

We agree that the age of the deposited OM is not zero. In the second application, age is not explicitly modeled, but deposited material has not a single value but a distribution. In application 3, we used equation 37 for freshly deposited OM, which is an age distribution (not a single value). However, we did not specify the parameters used in the test run, which will be added in the new version.

The derived PDEs apply to any distribution. However, mixing and reaction processes will let the distribution evolve into a different shape, which needs to be reconstructed in order to perform the integrations for the reaction terms. The fundamental problem is that the change in distribution shape means it no longer conforms to the distribution type of the deposited matter. Thus, when a gamma distribution is used for the reactivity distribution of deposited OM, the distribution will no longer be a gamma distribution at depth due to mixing. The challenge is to find a distribution that may work well at different depths.

In this manuscript, we do not develop a method to reconstruct distributions resulting from deposited OM with the gamma distribution, which is outside the scope of this manuscript. However, we will elaborate on using different distributions (see point 7) and show examples.

Original comment: Meile and Van Cappellen calculated particulate age and transit time distributions by particle tracking approach in bioturbated marine sediments. Could this model reproduce the results with same setup?

Answer: We have not directly compared our model results with this model. The Lagrangian simulation in the manuscript only serves to validate the diffusion equations (Table 1). In addition, we will add a supplement in the revised manuscript, showing analytical validation by using a delta distribution, which, together with the simple Lagrangian model and mathematical derivations, proves the derivations are correct.

Original comment: The authors emphasized the important of center-moments to derive the PDEs. Please gives the difference in detail between center-moments and non-central moments and how they converted to each other.

Answer: We will elaborate on the differences between using centralized and non-centralized moments in the equations for reactions (see response to reviewer 1). However, a derivation for PDEs for the diffusion of non-central moments was already given by Delhez and Deleersnijder (2002), which we will not repeat. The conversion from central to non-central moments was given in equation 28. We will also add the definition for non-central moments.

Original comment: In application 3.1 a particle tracking simulation has been used to compare current model, but the authors didn't mention how particle tracking simulation works in the text. Similar in application 3.2, it was not mentioned how 30-G model works. The same is in application 3.3 for the discrete simulation with 500 age bins.

Answer: We will elaborate on this in the new version. A description of the simple particle-tracking model will be given in a supplement. We will add text explaining how the bins for the multi-G simulations were designated (see response to reviewer 1).

Original comment: I found one problem in current model is to choose a distribution function to fit the moments. This choice is not unique and it must be numerically evaluated and time-consuming. The multidimensional root-finding procedure may fail. I think that this model lacks generality to model early diagenesis at marine sediment. Maybe the authors can run example with the real data from marine sediment.

Answer: Indeed, a distribution may not be uniquely defined by a limited number of moments, as is discussed in lines 272-275. The largest practical drawback of the approach is that it may be hard to reconstruct a distribution from moments (for which we used the root-finding procedure, eq. 32), but this problem is discussed in the text at length. See, for instance, lines 276-291. For the application in section 3.2, the simulations were stable and ran relatively fast. The simulation was slow and more challenging for the application in section 3.3. However, the last application is also included to show the limitations of the moments-based approach. Applications 1 and 2 could already be used in practice. Application 3 is a proof-of-concept rather than a final product. It facilitates the discussion and shows the limitations of the method. The emphasis of the manuscript is on the mathematical model and not the numerical implementation. However, at the end of the discussion, we give suggestions to improve the method's performance (lines 301-304).

Most early diagenetic models use 1-G, 2-G, or 3-G formulations (see Arndt et al., 2013). The second application demonstrates that our new model is up to the task, as it can reproduce results from a 30-G model and runs reliably and relatively fast. The question regarding data will be addressed in the following item.

Original comment: Although the author listed the application of the Eulerian model under three different conditions, the authors could select a specific site, including organic matter content data and organic matter 14C age data, to further validate the accuracy of the model.

Answer: This paper aims to derive general equations that can be used for reaction-transport models and are not limited to describing OM mineralization in sediments. An application of the model to measured TOC contents, considering good dating proxies, and elaborating on the most suitable distribution for OM reactivity, can better be addressed in a separate paper.

Original comment: The mathematics is rather complex to understand. The meaning of many mathematical symbols is confusing (e.g., Line 59: What does δ x mean?). Therefore, I suggest that the author make a table showing in detail the meaning and value of each mathematical symbol in the text.

Answer: δ_x is the jumping distance of particles. We will add a table with the meaning of all symbols (see point 3).

Original comment: The setting and handling of the upper and lower boundaries of the model is the key to solving the model. However, the description of the upper and lower boundaries in the text is not very specific. I suggest listing more details and formulas to show how the upper and lower boundaries are handled in the solution.

Answer: It is a good point, and we will add this.

Sincerely,

Jurjen Rooze

PS: All cited references can be found in the first submitted manuscript.