Overview

Van der Zant et al. present RADIv2, an updated version of the RADIv1 1D early diagenesis sediment model described by Sulpis et al., (2022). Specifically, they have included improved parameterizations for diffusive boundary layer thickness and pore water dispersion, and added methane cycling. They also added features that allow ensembles to be run, and that allow unknown parameters to be optimized to better fit in situ data. They validate RADIv2 from in situ measurements across a range of different locations, then use RADIv2 outputs to develop a metamodel. This metamodel parameterizes dissolved oxygen, dissolved inorganic carbon, and total alkalinity fluxes in terms of temperature, bottom velocity, bottom-water calcite saturation state, and organic and inorganic carbon fluxes. They suggest that this metamodel could be used to represent more complex sediment dynamics than are currently included in global ocean biogeochemical models (GOBMs), without incurring the significant computational expense of including a fully process based model.

General comments

I think this paper adds some very useful features to RADIv1 and seeks to address a critical issue by demonstrating a pathway towards including more realistic sediment biogeochemical cycling in GOBMs. It’s well written, and I think the authors have done a great job in finding so many datasets to compare RADIv2 output against as I know that useful datasets for calibrating such a model are not readily available. Overall, I think it’s a strong paper. Below are some suggestions primarily focused on validation and the associated discussion which could improve it further.

I think this paper could be strengthened by adding more discussion to the sections validating RADIv2 against in situ observations. The applicability of RADIv2 and the resulting metamodel to a wide range of environments depends on these sections, hence I think they are particularly important.

Comparisons between the modelled and observed fluxes currently focus on whether the modelled fluxes lie within the range of observational values. However, it seems like in several cases the model means lie the observed values. Performing some kind of statistical test to determine whether the model differs significantly from the observations or not would strengthen the case that the model does capture most of what’s happening in the environment, and that the variations attributed to transient processes aren’t as important. If they do differ significantly, it might indicate that these differences are due to more than just these transient variations, or that these variations occur so frequently that they need to be accounted for.

I also think the discussion about the discrepancies between the model output and the observational data could dive deeper into factors that might be affecting these differences on the modelling side. Currently much of this discussion is framed around how the values the model is being compared to are outliers in the wider literature (e.g., the North Sea values for remineralization rates were lower than rates reported from other shelves, the Berelson et al., (2003) values for Monterey Bay are lower than the global-mean C_org remineralization rates for continental shelf environments). I think that if the point of this section is to demonstrate that RADIv2 can capture conditions across a range of different environments then it’s important to focus the discussion on how it compares to the specific environment it’s being tested against rather than against more generic conditions.

Diving more into the factors in the model that could be driving these discrepancies could also offer more insight into ways that the RADIv2 could be improved. In particular, I think that more discussion of why RADIv2 overestimates the POC remineralization rate at each site would be useful, as it seems like these overestimations could be responsible for the modelled O2, TA, and DIC fluxes tending to be higher than the observed values. Could this indicate that the POC fluxes used for the model are too high? How are the remineralization rate constants chosen for each site, and how does changing them affect both the remineralization rates and the partitioning of remineralization between aerobic and anaerobic pathways? There are also other parameters that could affect these fluxes, such as changes in porosity or the temperature dependence of redox reaction rates. I understand that RADIv2 has a lot of dials to turn, and the data necessary to really constrain it doesn’t exist and so these questions are probably impossible to answer definitively. However, given the importance of this section I think more discussion about why it overestimates remineralization would help strengthen this work.

I think it would also be useful to check the metamodel against the experimental data, rather than only comparing it to RADIv2 output. The ultimate goal of the metamodel is to simulate the world, hence I think it’s important to see how it does against observational data in addition to how well it replicates the results of RADIv2.

The explanation of how the metamodel was developed could be made clearer by specifying which results the regressions were performed on. If they were performed on the same outputs used to validate the model, then I think it would be helpful for that to be stated. If that is the case, was there thought given to running RADIv2 over a wider parameter space and performing regressions on those results? One area that could potentially be beneficial to explore would be more tropical, carbonate rich environments. Assessing the performance of a metamodel developed over a wider parameter space against both RADIv2 outputs and the validation sites presented here would provide a better basis for assessing how well the metamodel would fare over the diverse range of environments present within a GOBM.

I think the structure of the paper could be changed slightly to improve clarity. Section 3.1 to 3.5 are all about model validation, so I would suggest renaming section 3 to model validation. I would then add a distinct metamodel section, which could include parts of what is written in section 3, and 3.6

Overall, I think this is a strong piece of work that highlights a path forward and makes valuable contributions towards addressing a large problem in current GOBMs. The suggestions above are intended to help strengthen the validation and discussion in order to further improve the work and maximize its impact. I congratulate the authors for tackling this difficult problem, and look forward to seeing this work develop further.

Specific comments

Line 59: generic is repeated twice

Line 107: How is this incorporated within RADIv2? Manually setting them based on literature values?

Line 112: Does RADIv2 use the same carbonate dissolution/precipitation scheme as RADIv1? If so, I think it’s worth stating here, as RADIv1 is setup for relatively cold temperatures. If not, I think it would be useful to describe the differences given the importance of precipitation and dissolution to sediment biogeochemistry.

Line 256: Move the first bracket to after ITTC et al., so that it reads … ITTC et al., (2011).

Line 296: It could be helpful to be a bit more specific here. Is it the climate change-induced warming itself that could outpace the capacity, or is it the effects of warming? If so, what effects specifically?

Line 336:

Line 348: A citation would be handy for the 12 mmol C m^-2 d^-1 number.

Line 392: The text refers to Luff et al. (2000) while the figure reports values from Epping et al. (2002). Should there be another reference in here?

Line 426: Capitalise the T variables.

Line 457: By keeping track of organic carbon pools with different labilities, won’t RADIv2 capture some of this implicitly? The more labile stuff will be consumed closer to the surface, resulting in a reduction in lability with depth.

Figures & tables

Figure 1: The fitted line in 1a is difficult to make out. It might be worth changing the colour and style so it stands out a bit more. Switching the panels may also improve the flow so that they appear in the order that they’re discussed.

Figure 2:Changing the x axis labels so they’re in alphabetical order would make it easier to find the reference each bar refers to. You could also put all of the references in a single table, either in the text or supplementary, so that you didn’t have to list them in each caption.

If the data’s available and n is large enough that they’re not too noisy, it would be interesting to see the bars as violin plots to better get a sense of the data. Are the long tails in the in situ data caused by just a few extreme points, or do some of these data sets of some kind of bimodality? Are the RADIv2 results relatively normal?

Could the table be presented as a figure as well? Given the number of different values reported it might make it easier to quickly compare different values.

Figures 3-5: Same comment about violin plots, though I think a table’s fine for just a few values. The figure panels for 4 and 5 could potentially be narrowed. The colours in figure 2 correspond to the parameter being plotted. It would be good to keep these colours consistent in the other plots.  

Figure 6: The x labels have RADI2v instead of RADIv2. All the other figures with oxygen fluxes (which I think were into the sediment) have reported them as positive values while here they’re reported as negative values. I think changing the other figures to have negative values would help make the presentation of this paper more consistent, and highlight that the oxygen fluxes are in the opposite direction to the TA and DIC fluxes.

Table 2: Is saturation state calculated from TA and DIC? If so, it might be worth stating it specifically given the larger uncertainties that can occur if you calculate it from other pairs of parameters.

It seems like all the organic carbon is now either fast or slow, whereas RADIv1 had refractory organic carbon. Has that been removed from the model? If so, I think this change should be mentioned somewhere.

Overall Comment

The paper by Van der Zant presented an update to the previously published Radi early diagenetic model. This version of the model improves upon first version by including relevant processes related to the exchange of solutes across the sediment-water interface (SWI) as well as enhancing the model capabilities for simulating biogeochemical processes in highly permeable sediments. A key feature of this model and the overall paper is the showcasing of the adaptability of the model in computing benthic solute fluxes which can be used in coupled benthic pelagic application from a metamodel built from RadiV2. This parameterization of benthic solutes exchange from metamodel is important in the context of better representation of seafloor biogeochemical dynamics in global earth system model. Overall, the paper is well written and concise. The paper take care in explaining the relevant addition of the new development of Radi and showcasing its utilizing in estimating benthic fluxes across different marine environment. However, the paper will do well in at least providing supplementary text on the performance of this new version in capturing the variability of porewater profiles especially in coastal regions where it has previously not been tested. This will help inform users of the strength of the new model as well as the current limitation as highlighted in the latter part of the paper for their application. Overall, this is solid work, and I recommend for publication following some minor comments below.

General Comments

I agree with the first reviewer about the statistical validation of the model derived fluxes in comparison to observations. One other thing I might add as well is the lack of detail in how the five environmental drivers used in the empirical formulation of the metamodel was selected. Was there any stepwise selection screening to diagnose what variables are sensitive to benthic fluxes of DIC, O2, TA? Were the selected variables based on prior literature findings (I can’t find any cross-reference in the text)?

Also, it will be beneficial to comment on the computational performance Radiv2 compared to version 1 as this version has exploited Julia’s differential equations packages. I agree with the comment of improvement in computational efficiency of the newer version of Radi but details on how much and how well this improvement is will be useful to help guide potential users in their expectation as other competing framework are growing in this space (see SedTrace.jl – J Du, 2023).  

Specific comments

L71: How does RadiV2 steady state work. Julia differentialEquations.jl can deal with steady state but depending on the dynamic itself (eg oscillating equilibrium, unstable equilibrium), it can be difficult in finding solution space that is well-posed (Steady State Solvers · DifferentialEquations.jl). The paper state nothing about the numerical scheme employed in this new version of Radi that makes it different from the previous version (which I believe was a first order Euler-based scheme).

L100: Is the dz constant, if so, what interval is it? If it is constant dz over the depth domain, then the resolution of dz is important if RadiV2 is employed in coastal sediment as processes operating in the upper few cm of the sediment drive the computed benthic fluxes. Thus, not properly resolving the underlying diagenetic processes at finer spatial scale will be detrimental to the predicted benthic flux across SWI.

Table 1: My impression of RadiV2, is that it can simulate coastal early diagenetic processes. However, in these settings, the role of anaerobic pathway featuring metals (Iron in particular) coupled to the carbon cycle is strong. This will necessitate modelling secondary redox processes relating to metals (Iron and Manganese, sulfide and its minerals). In iron rich sediment, for example there can be significant routing of reduced sulfide to mineral phase (e.g pyrite) which can be either buried or undergo redox oscillation and ultimately alter alkalinity budget in marine sediment (Canfield, D. E., et al. (2005)). Is there a justification why this is omitted redox pathways are not included in the model?

Section 2.1.2, I think the organization of the paper section can be improved. Generally, the paper can provide more information on the computational cost associated with the model improvement relative to the previous version.

L180 Eqn 6, what does b0 mean? No mention of the term in the text and it appears again in eqn 12 (L215). My assumption is it related to surface bioturbation rate, but I might be wrong.

Same for eqn 8 (L200). The terms in the equation are not clearly defined. What does C2 mean? Concentration in the second layer (or what).

L258 I think the cross-reference equation numbers are wrong. Eqn 10 is related to bioirrigation and has little to do with bottom water viscosity. v should be substituted in eqn 14. Am I right in this assertion? If so, please fix the numbering.

 In the practical application, is U in Eqn 14 derived from the terms (k, v, u*) and calculated implicitly or is it defined by user.  Browsing through the code posed in zenodo, it seems U is imposed by the user or application. I suppose both options are valid.

Table 2: It seems that the depth considered here only validates the model for continental shelf/and adjacent coastal as well as deep Ocean. For delta and shallow coastal area, how will the coverage of the ensample of tested region perform. 

Are the values for pFast/slow as well as kinetic constant literature based? from previous modelling/experimental work done in these sediments or are they arbitrarily chosen?  

Discussion

In the new version of Radi, it was stated that it can perform transient simulation, however in some of these regions there can be significant variability in physical (tide, resuspension/erosion) and biogeochemical forcing (changing boundary condition) that either need to impose or modelled. Did you test some configuration of this new version of the model in simulating these transient dynamics? Perhaps for on some stations. How were the results. Perhaps including some of that in the supplementary text will be useful.

Summary of RadiV2 performance: Here, a short sentence on the computational performance of RadiV2 as compared to the previous version would be useful.

References:

Canfield, Donald E., Erik Kristensen, and Bo Thamdrup. "The iron and manganese cycles." In Advances in marine biology, vol. 48, pp. 269-312. Academic Press, 2005.

Du, Jianghui. "SedTrace 1.0: a Julia-based framework for generating and running reactive-transport models of marine sediment diagenesis specializing in trace elements and isotopes." Geoscientific Model Development 16, no. 20 (2023): 5865-5894.