General comments

The manuscript investigates the effects of benthic fauna on the biogeochemical transformations of major elements in the Baltic Sea via direct (metabolic) and indirect (bioturbation) pathways, by coupling a recently developed model of benthic fauna to an established hydrodynamic-biogeochemical model. Within regional ecosystem models, benthic components are often represented in a simplified way, where their role is largely limited to returning remineralised nutrients back into the water column. In the Baltic Sea, models accounting for indirect impacts of benthic fauna have been used for a while now, incentivised recently by the need to study the effects of bioturbation by invasive species. However, these models did not consider benthic faunal biomass explicitly, and thus did not allow a separation of effects of bioturbation and direct effects of faunal metabolism on benthic-pelagic dynamics. This manuscript bridges this gap and provides a detailed look at the role benthic fauna plays in shaping biogeochemical fluxes in the Baltic Sea. This study builds on the authors’ previous work and contributes to the advancement of benthic system modelling, which has lagged behind its pelagic counterpart in coupled regional models. Although it is an important step in the right direction, I have several criticisms I would like to bring to the attention of the authors. Foremost they are related to the “large-scale” approach taken for this work. It should not serve as a substitute to thorough comparison of model results with data, nor to the detailed presentation and discussion of results (providing quantitative estimates and aligning text narrative with figures). Limitations of the model, unique contributions of this work to the understanding of the Baltic Sea biogeochemistry beyond state-of-the-art knowledge should be appropriately discussed. Below are my detailed comments and questions, that I hope will improve the quality of the manuscript.

Specific comments

Title

1.    The title should clearly indicate that it is a model-based study.

Abstract

1.    It would be informative to add quantitative estimates of changes due to direct and indirect contributions of benthic fauna. “Small proportion of seafloor organic stocks”, “decreases denitrification”, “increases P retention”, “reduction in N fixation” etc all should be accompanied with quantitative estimates (% change). The way information is currently presented, it is hard to see what a unique contribution of this work is to the understanding of biogeochemical impacts of benthic fauna in the Baltic Sea, as the effects of bioturbation via sediment oxygenation on N, P dynamics in the Baltic Sea are generally well known from previous studies.

2.    Line 16: It should be clear from the text that bioturbation affects denitrification, P retention etc indirectly, e.g., increasing oxygen penetration depth and availability. 

3.    Line 20: “chain of indirect effects” is opposed to “direct effects of faunal respiration, excretion and bioturbation”. This create some confusion over which effects are “direct” and which are “indirect” ones.

Introduction

1.    Lines 43-44. “The combined effects of animal bioturbation and metabolism have seldom been studied together.” Several references support this statement, but could the authors elaborate on what are the main reasons for that gap in research, despite recognition of the importance of benthic processes?

2.    Could you provide a more solid justification for your approach? Why do we need to implement modelling – does it address the knowledge gap? Why did biogeochemical models of the Baltic Sea so far did not include benthic fauna?

3.    When specifying reasons for using Baltic Sea as a model area, to what extent also relatively simple benthic community composition play role?

Materials and methods

1.    Line 74: the focus of the study is “on the Baltic Proper and the Gulf of Riga”. Could you please explain your choice of focus areas in more detail? While the focus on Baltic Proper as the largest and deepest basin is straightforward, why Gulf of Riga? It is the only sub-basin to the east of the Danish Straits for which there was no validation data available apart from some literature-based values, so it is not possible to adequately validate the model for this region, and consecutively the confidence in model performance here is lower than in other areas, both in terms of absolute biomasses of benthic fauna and its impact on pelagic biogeochemistry.

2.    Figure 1: As the text often refers to different regions of the Baltic Sea, a legend for all of them should be provided in addition to numbers. Moreover, model results are presented at intervals of 0-30, 30-70, 70-120 (and 120+) meters, so bathymetry on the figure should use the same gradation.

3.    Section 2.2: BALTSEM simulates Baltic Sea as 13 horizontally homogeneous boxes. In the Results and/or Discussion section, it should be elaborated on how this affects model results. Some of the Baltic Sea sub-basins have strong gradients in nutrient distributions, which in turn leads to gradients in productivity and in distributions of benthic fauna. For instance, the eastern part of the Gulf of Finland is heavily influenced by riverine nutrient inputs and have higher primary production rates than its western part. How are these gradients accounted for? 

4.    Lines 128-130: why does the model consider respiration and excretion of benthic fauna to contribute ammonia and phosphate directly to the water column? At least in case of deposit feeders one might expect part of the excreted ammonia to be oxidised directly within sediments and to be released into the water column in form of nitrate, and some of the excreted phosphate to be bound within sediment, the same way as it is considered for microbially remineralised nutrients in the model?

5.    Does the model consider sediment resuspension? It seems to be an important factor mediating organic matter availability in the sediments in coastal seas, especially shallow regions such as Gulf of Riga? Omitting it might have important consequences both for model parameterisation and its results?

6.    Line 134: “degradation” - what type of degradation? At this point, it should be specified in more concrete terms.

7.    Around line 152: how are silicate transformations handled in the benthos?

8.    Line 167: is nitrification considered a sink for oxygen in the model?

9.    Equation 8, page 7: is sequestered phosphorus considered a state variable in the model, as it is not listed among the state variables in Table A1?

10.    Line 176: what is a definition of “severely hypoxic”?

11.    Section 2.3: please provide more detail on model forcing and setup in addition to the references, so the interested reader does not need to look for those details within several previous papers. How was benthic fauna initialised?

12.    Line 184-185. Aggregating the results to means and standard deviations obscures a lot of detail about model capabilities, which cannot be justified by stating that “the purpose of this study was to evaluate large-scale dynamics”. Does the model capture seasonal dynamics of benthic fauna? Are there long-term trends in the model or in the data? There is a general lack of knowledge about benthic processes, so more detailed and varied comparison with data would allow to identify critical gaps and steer discussions which will help to identify directions for improvement. This is especially important for a new model implementation, as presented in the manuscript.

13.    Lines 187-188. This is an important topic for benthic modelling and should be discussed in some more detail. What should observational scientists measure to help constrain the model?

Results

1.    The model run is 1970-2020. Why was validation of pelagic state variables performed for 1970-2015?

2.    Figure 8 shows that there was a noticeable impact of inclusion of benthic fauna/bioturbation on primary productivity and nitrogen fixation, so even if overall relative bias remained almost unchanged, there was an impact on pelagic dynamics. Please provide more detail on how the models with BMM and without BMM compare. Given the impacts of fauna and its activity on nutrient fluxes is the main scope of the paper, a more detailed assessment of changes in pelagic environment is justified.

3.    Line 214: how is “reasonable accuracy” defined in case of benthic variables? As the authors demonstrate, benthic fauna biomass shows high variability, hence large standard deviations, especially when aggregated over long time and large areas, while standard deviations in the modelled fauna are relatively small. Does this high std in data really justify using cost function (Appendix C) as a validation metric?

4.    Line 221-224: see my comment above: why focus on the Gulf of Riga?

5.    Section 3.3, from line 247: this is a key section in the manuscript, as it demonstrates the direct and indirect effects of fauna. Yet it is very short and just skims through the results and through the figures. It should be more detailed, which would also make it easier to relate figures 7-9 to the text. It would also benefit from detailed quantitative information, in particular regarding relative changes.

6.    Line 264: for this, a figure similar to figure 5 could be provided.

7.    Lines 265-271: figure 11 (as well as several other figures of the manuscript) contain multiple subfigures, so it would be very helpful to provide relevant pointers in the text, e.g., reference to Fig 11a on line 268, Fig 11b on line 269 etc, so the text and figures could complement each other. 

Discussion

1.    Line 284: both terms “long-term” and “large-scale” have been used multiple times in the manuscript without proper definition. It creates overall impression of vagueness. Does “long-term” stand for “multidecadal” or “long-term mean”? The words “large-scale” could be omitted altogether without impact on meaning.

2.    Line 285: as it currently stands, it should be “Baltic Sea”, not “coastal sea”, as the manuscript shows the model is not yet applicable in high-salinity regions.

3.    Line 288-289: to support the statement that benthic fauna can alleviate the ‘vicious circle’ of eutrophication, could you show the differences in extent of hypoxic area (e.g. with Figure 8) for simulations with different levels of bioturbation?

4.    Similarly, what about primary production and N-fixation for the two future scenarios (Fig 10)? I think these would also support the discussion on faunal impacts on ‘vicious circle’?

5.     From line 307: as the overall model performance was not improved by adding benthic fauna, and related processes are not well understood, could the authors discuss on what understanding is currently lacking and which model improvements are desired.

6.    Lines 313-316: same as the previous comment, which data is currently lacking? Which mechanisms are not sufficiently understood? The sentence on “the main strength of this study” is too general (“many interlinkages”) – could it be replaced with something more concrete and relevant to the actual modelling work?

7.    Lines 326-334: please keep referencing to figures and sub-figures as discussion requires.

8.    Line 361: is it worth presenting fast running time as an advantage, since the manuscript presents results from only several simulations with sensitivity analysis limited to a single parameter?

Conclusions and outlook

1.    Line 397: “much-studied bioturbation” - I suggest replacing this with “bioturbation, relatively more studied in the modelling context”.

2.    Could the authors provide some outlook on the potential future directions of benthic modelling, in the context of Baltic Sea in particular?

3.    Lines 399-402: These conclusions are far too general. Could it be something emerging from the study rather than vague statement that benthic-pelagic coupling is “modified by multiple drivers, which may change over time”?

4.    Based on your work, has the time arrived for the regional models (of the Baltic Sea) to extend their formulations to explicitly include benthic fauna?

Appendix A, line 817: why is mortality rate chosen to be linear for Limecola balthica, and quadratic for the other two groups?

Technical corrections

1.    Line 89: “extension” - should be “extent”

2.    Line 110: “Ponotporeia” - should be “Pontoporeia”

3.    In equation 1, should “Uci” be “UBFiC” in the denominator?

4.    Page 6: equation 6 is missing.

5.    Line 175 and equation 10: the fitting constants 5-8 are not featuring in the equation.

6.    Figures 7 and 11: in both cases, “phosphorus (c, d)” should be “phosphorus (e, f)”.

7.    Line 319: “suggests” - should be “suggest”

8.    Table 1: what is “total” biomass? Should it be “mean” instead?

9.    Line 712: “and” - should probably be “an”

10.    Line 742: Nitrite and nitrate should be “NO2” and “NO3”, respectively.

11.    Line 784: “and” at the end of the line should be “are”?

12.    Line 794: “feces” - “faeces”.

13.    Line 907: “120-120” should be “70-120”

14.    Figure C2: please use a different y-scale for the Bothnian Bay.

In this paper, a new tool is presented that couples a low-resolution pelagic biogeochemical model with a low resolution benthic biological model for the Baltic sea. The (vertically integrated) benthic state variables are then used to calculate impacts on biogeochemistry using presumed effects of bioturbation and water-column conditions on denitrification and phosphorus dynamics. My main doubts with this paper are connected to the biological focus of the model.

Essentially there exist two schools of modelers: some modelers take a biological approach and ignore or strongly parameterize biogeochemistry. Their models disregard the small-scale vertical gradients of solutes in the sediment and often consider only surface-averaged concentrations of particulate substances (e.g. organic matter). Moreover, their models operate on seasonal time scales, as organisms usually react on these time scales.  Opposed to this are the modelers that tackle sediment dynamics from a biogeochemical perspective and strongly parameterize biology. These modelers take into account the fine-scaled vertical gradients of solids and solutes that are observed in the sediment, and their dynamics includes reactions operating at very different timescales, from very short (< seconds) up to very long time scales (multi-years). In these models, the metabolism of the (higher) organisms is included as “oxic mineralization” of organic matter, while their bioturbation activity is included as a “coefficient”. Thus, these models strongly parameterize the biology, and only explicitly account for the biogeochemistry. As long as the main conclusions of these models are stated in the area of the model focus, there is nothing wrong with any of these approaches. For instance, it is reasonable to assume that a biogeochemical model can rather faithfully reproduce the impacts of certain external conditions on sedimentary nitrogen or phosphorus removal rates, but it is questionable whether such models can also well represent the distribution of the benthic organisms that drive the biogeochemical cycles. Similarly, why would we put a lot of faith in biogeochemical conclusions that come from a model that focusses on biology and parameterizes the biogeochemistry? This is in a nutshell the doubts I have on this paper. While the conclusions seem logical, I am still to be convinced that the tool used to arrive at them is appropriate.

Because of the biological focus, there are quite some assumptions with respect to biogeochemistry that are not dealt with in the manuscript. For instance: the paper talks about the sediment pools of C, N and P, and Si. Biogeochemically one distinguishes between particulate and dissolved pools – here I had to guess that the pools refer only to particles (the ‘food’ of the organisms). Thus, the transient (within season) storage of dissolved components is ignored. Is this a reasonable assumption ? (I could not find any evidence for this). In addition, historical eutrophication in the Baltic may have caused significant storage of dissolved nutrients deep in the sediment (i.e. ammonia, phosphate, sulphide), which are not accounted for in the model. Can these be ignored – what is the effect of ignoring these on long-term simulations?

In addition, the dependencies of the biogeochemical processes on the model variables are so complex that it is very difficult to see how these processes are affected. For instance the formula (5), which essentially describes the dependency of denitrification on water-column oxygen and biota, has 4 “fitting” parameters – to what data have these been fitted? The P-sequestration formula (formula 7) has even 8 “fitting” parameters. On line 187, it is said that it is difficult to constrain the new parameters. Does this mean that these parameters have not been fitted at all – and if they have, on which data? And why would instead running sensitivity analyses by changing the Ebio parameter be a valid alternative? A little more effort in showing that these dependencies are realistic is required. (and where is formula 6?).

I also find the lack of any comparison of model output with biogeochemical sediment data worrisome. On L 313, the authors claim that they cannot “properly validate the simulated sediment stocks or fluxes due to a lack of large–scale data and insufficient understanding of the multitude of mechanisms underlying the biogeochemical transformations and fluxes”.  The first part (lack of data) does not do justice to the multiple biogeochemical studies in the Baltic that have recorded sediment-water exchange fluxes, and measured sediment concentration profiles in great detail. Also, I do not agree with the statement that there is “insufficient understanding” of biogeochemistry. As a quantitative science, biogeochemistry is at least as (and probably much more) advanced as biology!  And even if it were true that we do not understand the biogeochemistry, why would we then trust the simple parameterisations that are used in this manuscript?

In summary, as much as I like the conclusions from this paper, the authors need to try a bit harder to convince that biogeochemistry in the Baltic can be predicted based on presumed effects of biological activity on N and P removal.