Ye et al. present in their paper the coupling of the process-based Model of Early Diagenesis in the Upper Sediment (MEDUSA v.2) to an ocean biogeochemistry model (FESOM v.2.1 with reduced model resolution + REcoM v.3 in reduced complexity) as well as first results from an pre-industrial simulation with prognostic atmospheric CO2.

The paper is clearly written and the coupling of a process-based sediment model to an Earth system model is outlined nicely and is an important step, especially when investigating the long-term carbon cycle.

Therefore, I support publication in GMD and hope the authors will find my few comments below helpful and consider their implementation.

General comments:

While I understand that the coupling of model components and work on model code in general can be very time consuming and that the focus of this study is the documentation of this coupling, I still think it would be nice to see a little more results. 

The authors state that they use a reduced complexity version of REcoM3, REcoM3p, that is targeted for paleo simulations. Also, ocean-sediment interactions become especially interesting when looking at long timescales, such as in paleo-simulations.

Therefore, I think it could enrich the manuscript to see some snapshots of the coupled model under, e.g., LGM conditions and update Tables 2 and 3 with the corresponding LGM values.

In such an exercise, carbon isotopes would be of interest as well...

Further. during comparison with observational data, the coarse resolution of the PI mesh is mentioned as a limiting factor (for example l. 246-249, 280-289) and in the conclusions an outlook is given to stay tuned for not only carbon isotopes but also higher spatial resolution. In my eyes, the manuscript could further benefit from including some results of this ongoing effort, if possible, and not save it all for future publications.

Specific comments:

• You could consider illustrating the carbonate chemistry by providing the governing equations for a better overview in the introduction.
• Figure 1: maybe add riverine + dust inputs to figure to close the loop?
• Section 2.4.4 is not fully clear to me: Does the reported performance apply to REcoM3 or REcoM3p and to FESOM2.1 with reduced model resolution?
• Section 3.1: you write that global vertical profiles in the model agree 'rather well' with observations from GLODAPv2. Could you give metrics?
  
  Are there differences between basins? Maybe add profiles for the different basins?
  
  Maybe also add section plots (model, obs, difference) for the main ocean basins.
• Mass conversation in the coupled model: It seems to me that this definitely needs to be addressed for longer (paleo-)simulations!

Technical corrections:

• l. 50: More complex scheme -> More complex schemes
• l. 111: 2) selecting of processes -> 2) the selection of processes
• l. 112: 3) writing the resulting -> 3) the writing of the resulting
• l. 127: were than partitioned -> were then partitioned
• l. 267: are reproduced in the model -> are reproduced to some extent in the model?
• l. 278: The opal belt in the equatorial eastern Pacific is smaller and less pronounced in the model than observed -> not visible to me
• Caption figure 3: Horizontal averages of -> Global horizontal averages of ?
• Figure 6 and 7: label 0 seems misaligned in colorbar
• Caption table 2: Note that the units here are Tmol year-1, not Pg year-1 -> Note the different units in the table
• Table 2: add observational estimates where available

Review of “FESOM2.1-REcoM3-MEDUSA2: an ocean-sea ice-biogeochemistry

model coupled to a sediment model”

by

Ying Ye, Guy Munhoven, Peter Köhler, Martin Butzin, Judith Hauck , Özgür Gürses, and

Christoph Völker

Summary:

Ying Ye and colleagues report on the coupling of the early diagenetic model MEDUSA2 with the ocean biogeochemistry model FESOM2.1-REcoM3. In order to be able to spin-up the model for multiple millennia (i.e., until the sediment-water interface (SWI) in the deep ocean is in steady-state), the authors use a lower horizontal resolution and reduce the complexity of the marine ecosystem model (i.e., simply using one generic zooplankton and detritus class instead of two for each) compared to a very recent model development paper (Gürses et al., 2023). Currently, most global Earth system models poorly represent the coupling between ocean and sediment biogeochemistry. Because the presented setup explicitly addresses the coupling of these domains, it can potentially be a very useful tool – especially for paleo-applications and simulations studying climate and marine biogeochemical feedbacks over multiple thousands of years.

While the model coupling itself represents a substantial contribution to Earth system modelling, the manuscript, unfortunately, lacks a proper evaluation of the performance of the new model setup and has several other weaknesses, omissions, and confusing parts. Not much new model development has been done for the manuscript – as the authors report on the coupling of two existing models. This would be okay if extensive experiments of the new coupled model evaluate its performance properly and show the added value of the new setup. Unfortunately, neither is done here. The authors only perform and show two experiments: one with the previous one-box sediment representation and one where they couple MEDUSA2. Both experiments are run under pre-industrial pCO2 for 2500 years (the new coupled configuration is only run for 1500 years from the previous model setup). Then, the authors compare some features at the end of both runs (Fig. 6, 7). The lesson learned from the results – apart from that patterns and values are slightly different - is unclear to me (e.g., is it a crucial improvement?).

Therefore, I cannot support the publication of this manuscript in Geoscientific Model Development. I hope my general comments will help to improve the useful coupling exercise and the evaluation of the new configuration. I suggest reconsidering the manuscript after major revision.

General comments:

1. The new coupled model is not thoroughly evaluated. I understand that the lower FESOM resolution and the reduced ecosystem complexity in RecoM represent an entirely new configuration. Therefore, a more in-depth evaluation of the model than comparing it simply to global depth profiles of DIC, DIN, Dsi, and ALK (Fig. 3 – the Fig. caption states these are horizontal averages?) and showing a time-series of atmospheric pCO2 (Fig. 4) is necessary. It is not too surprising that the model gives a good match to the globally averaged GLODAP data, considering that it is initialized with it and the model is only run for 2500 years. It is necessary to show that the model – particularly seafloor conditions and the sediment-water interface (SWI) fluxes – is properly spun up. To evaluate if the model is in steady-state, I suggest including time-series plots of global properties such as global mean ocean O2, nutrients, DIC, DSi, ALK; export production, settling and burial fluxes of POC, CaCO3, opal, SWI fluxes of dissolved O2 and nutrients – potentially also concentrations of dissolved species at the seafloor. Maps of NPP and basin-averaged meridional-depth distributions of DIN, O2, and ALK (compared to observations) would further increase the credibility of the model.

2. The manuscript does not show and/or make use of the features that are added by coupling FESOM-REcoM to MEDUSA. MEDUSA2 simulates OM degradation with O2 and NO3. It would be interesting to see a map of the fraction of aerobic OM degradation vs denitrification. This should show a clear difference between the deeper ocean and shallowe sediments with more OM input.

As mentioned in my first comment, I would like to see maps of simulated SWI-fluxes (e.g., of O2, NO3, DIC, ALK). [2.1. As a side-note: What is done in MEDUSA2 when NO3 is exhausted?]

As described in section 2.1.2, MESUSA2 not only simulates a reactive surface sediment layer but also a core layer that records synthetic sediment cores which is a fantastic feature for paleo-applications. It would be very informative to show simulated sediment core layers for different ocean depths (e.g., a shelve vs deep ocean core) for instance during a carbon perturbation experiment.

3. Related to the previous point: The authors want to make REcoM3 ‘fit for paleo-applications (see 2.2.2). Carbon isotopes are of particular interest for paleo-applications and, in my (and also the other reviewer’s) opinion, should be included in this manuscript and not in a (very short) additional publication of pretty much the same authors (Butzin et al., EGUsphere). The reduced complexity configuration of the model here is particularly useful as long spin-ups are necessary to reach an equilibrium for the isotope system. I understand that carbon isotopes are currently being developed in MEDUSA2 and this might not be straight-forward in a vertically resolved diagenetic model. If this feature is not yet available, the sediment coupling of C-isotopes could for instance be simply realised by assuming no fractionation during OM remineralisation etc. in the sediment and calculating:

DIC_13C_swiflux = DIC_flux_OUT / POC_flux_IN * POC_13C_IN

4. The the previous one-layer sediment model box, R_sedbox, is unclear to me. So R_sedbox is not really a reflective boundary but it is also not really a sediment model either – hence, there is no benthic preservation simulated with R_sedbox (lines 59-65). It would be good to clarify how the sediment box calculates the return fluxes differently compared to a reflective boundary condition. I just saw that this is described in the appendix of Gürses et al., (2023) but it would be good to also include it in this manuscript as it is necessary to understand R_sedbox.

5. Related to the previous comment – Section 3.1 and table 1 is confusing. The fluxes given in Table 1 are very confusing – it is not 100% clear if these are settling or burial fluxes. The title of Tabel 1 says sinking fluxes (so settling fluxes onto the seafloor?) but the text refers to “calcite burial fluxes (line 249, 252; but then I thought R_sedbox does not simulate any preservation?).

So I suppose the model estimates are settling fluxes. But some of the observational estimates are clearly burial estimates (e.g., CaCO3 data estimates are burial fluxes, as stated in Cartapanis et al., 2018). Also POC data estimates stated (50 – 2600 PgC kyr-1) probably refer to burial rates and are confusing. First, where does the 50 actually come from? I know Cartapanis et al. state it but Burdige (2007) gives 160 PgC kyr-1 as the lowest value – and these are POC burial fluxes in these publications.

Often cited POC settling rates are 2628 PgC kyr-1 (Burdige, 2007), 2290 PgC kyr-1 (Dunne et al., 2007), or 930 PgC kyr-1 (Muller-Karger et al., 2004).

So if the model estimates are really settling fluxes, as suggested by the title of the table and the text (see, e.g., line 240), then these values are too low and not well distributed between shallow and deeper ocean. I find the argument that the coarse resolution is responsible not convincing. Previous models with similar or even coarser resolution are able to simulate POC (and calcite) settling fluxes and preservation on the shelves much better (e.g., Palastanga et al., 2011; Hülse et al., 2018).

Assuming that model estimates in Table 1 are settling fluxes – one cannot judge how well burial rates are simulated by the model. If the text is correct (in that these are calcite burial rates) then the CaCO3 burial rate in deep sea sediments (~0.3 PgC yr-1) is 2 – 3 times larger than budget estimates (0.1 – 0.15 PgC yr-1). It would be helpful to know the mean surface sediment CaCO3 content (vs observed 34.8 wt%).

In summary: Please make sure the model estimates are compared to the correct observational estimates. Also, it would be very helpful to include export, settling and burial fluxes for POC, calcite and opal in the table. And also please distinguish between settling and burial fluxes for sediments shallower and deeper than 1000m, as done for the calcite data estimates in the table. This will hopefully help to understand what causes some of the mismatches in POC preservation.

6. The POC wt% and the spatial distribution look not very convincing! Large areas show POC wt% > 5 (what are the maximum, mean values in these areas?) where observations show much lower values – mainly at high latitudes. In contrast, other areas were the data shows higher POC wt%, e.g., the major eastern boundary upwelling zones and the Arabian Sea, R_medinit does not simulate any OC preservation (Fig. 5).

[Also why are the observations compared to R_medinit and not to the final results of the coupled model? I know R_medinit is compared to R_coupled in Fig. 8 but I don’t find this comparison very informative.]

I suspect that the simplification to only use one class of detritus (line 78) in the water column might be partially responsible for the poor representation of POC wt% in the sediments (but it is impossible to be sure since no maps of POC settling fluxes are shown). The main reason however might be organic carbon degradation as simulated in MEDUSA2. MEDUSA2 simulates two classes of organic matter to approximate the different C:N stoichiometry of POC, right? What are the degradation rate constants for these classes? Is the more C-heavy class remineralized more slowly?

I would argue that more tuning of parameters (degradation / dissolution rate constants and/or other boundary conditions) are necessary to improve the model-data fit.

E.g., what about sedimentation rate: The terrestrial clay input of 2.5 E-8 mol cm-2 y-1 is spatially uniform. But should this not, as a first estimate, decrease with distance from the continents? This might help with the unrealistic distribution of carbon burial between the shallow and deep ocean (as stated in table 1 and the related text).

7. The manuscript does not include any parameter values. A comprehensive table stating parameter names, values, units, and references is necessary to understand how the model is set-up. The same applies to the riverine (i.e., weathering) inputs stated in Table. Please indicate where the values come from and how they compare to observational estimates.

8. I suspect, that the loss of N (i.e., 0.8% over 1500 years of simulation) will very likely be a problem during longer model runs if not compensated for via N2 fixation and/or weathering input – especially during paleo-applications with larger contributions of denitrification. Is suggest to fix the N-leak.

9. Sedimentary source of iron: The text says (line 82ff.)

“The sedimentary source of iron can be calculated in two ways: 1) in a fixed ratio to degradation of particulate organic nitrogen (PON) in the benthic layer as described in Gürses et al. (2023, Eq. A77 in Appendix A) or 2) in a fixed ratio to the diffusive flux of dissolved inorganic nitrogen (DIN) calculated by MEDUSA2 in coupled simulations.”

Can you please provide a justification for these representations and how well they approximate realistic SWI fluxes of iron. Also how well do they represent Fe fluxes under anoxic conditions – where Fe-cycling behaves very differently. This might be important depending on the paleo-applications the authors have in mind.

Also, A77 does not exist in the Appendix of Gürses et al. (2023).

A few minor comments:

A better motivation & explanation could be included why this model is appropriate for paleo-applications; also ginving potential applications. Ideally this would be compared to existing Earth system modelling approaches, highlighting the benefits of this new model configuration.

Table 2: Why does seafloor deposition (POC + CaCO3) not equal diffusive C flux out of sediment for R_sedbox?

Section 3.3.3

Unclear where the 402 PgC come from an how it compares to the 21 PgC in table 3.

Unclear what should be learned from the last bit of the section, i.e., the discussion of the carbon, alkalinity and silicon inventories not being in steady-state in the coupled run.

Some of the methodology is unclear – I thought “FESOM2.1 was run for 1000 years to spin-up ocean circulation.” (line 201) why does the 2500 year run in Fig. 4 show again ocean circulation stabilisation?

It is also confusing that two experiments are named R_sedbox – one described in 2.4.1 and then one described in 2.4.3 (and shown in Fig. 4?) The description of R_coupled is also not 100% clear. The text states: “(2) a coupled simulation R_coupled was conducted for 1500 model years first using the output from R_medinit as sedimentary input of DIC, Alk and nutrients.” (lines 231-214) I suppose this means, you start with the SWI-exchanges calculated in R_medinit (i.e., for the first 50 years), after that MEDUSA is called every 50 years.

Lines 268 ff. What does this refer to? How is this different to Fig. 5?