Review of Barré et al. 2023b : Implementation and assessment of a model including mixotrophs and carbonate cycle (Eco3M_MIX_CarbOx v1.0) – Part 2. – Julien Palmiéri.

Barré et al’s study make use of a complex biogeochemical model – developed and extensively used in the Mediterranean sea – introducing a new component to it in the first part, and analysing in great detail the carbonate chemistry of the bay of Marseille.

To do so, they use the model in a dimensionless version, representing a 1m³ surface box of the bay.

The analysis in this study are well driven, investigating all the different mechanisms impacting the air-sea CO₂ fluxes through a whole year, including the impact of fresh water inclusion from the nearby Rhône river, and summer upwellings.

Although the analysis are well done, and are interesting, the whole study is mined by the choice of the dimensionless mode.

* First, the reason of the choice is not given. In the introduction, it’s made mention of the need for high resolution model for coastal regional study, and the next sentence announce the use of a dimensionless (0D hereafter) configuration. I understand it is greatly needed for developments like adding the mixotrophs in Eco3M. It’s a huge task that require such lightweight configuration to test, verify that nothing is broken, and make sure the fluxes between each element of the model are reasonable. So this choice is easily explained for the part 1, but less for this one.

* A second point is that this 0D brings more questions than answers. Because it is a surface box, the model does not represent advection and mixing. The physics variable/forcing come from observations and hence include annual cycle and external forcing, including specific phenomenons like summer upwellings or the Rhône waters passing by. But what about the nutrients ? What are the external forces driving the biology of the model ? This is not explained until the discussion, what is extremely frustrating, as we don’t really understand what the model sees and feels or not, until the very end, when the author reveals some of the experiment limitations.

* Still about the 0D, what happens to the POM ? Do they sink ? Are they removed from the surface box ? Or do they float there and are slowly remineralized (as if the bay is mixed enough to keep the particles around) ? This is important as it has an impact on TA and DIC and all other nutrients concentration.

* Somehow it looks like (and I am sorry to say that, but I am sure you agree with me) the work you’ve done here (changing from autochtonous TA formulation – what is what you ideally want to use – to the abiotic, allochthonous formulation) is a way  to fix a problem due to  the configuration choice, that is not done for this kind of study. Your conclusion (you need to switch from 0D to 3 or at least 1D) should have been one of Lajaunie-Salla et al. 2021’s study.

These are my main concerns. I do recommend the publication of the paper as it is a nice study. The authors investigate the different mechanisms controlling the CO₂ fluxes, and manage to get surprisingly reasonable results, despite the use of the 0D. On top of that I would think it gives a good example, in the first part paper, of when to use the 0D and when not to, in this one.

Before publication I would require the author to better explain the choice and implications of the 0D in the method section, so that the reader can really understand the experiments and the results. How are the nutrients managed (initialized with annual average value)? Are total N, P, SI, Fe, Alk supposed to be conserved within the box ? Or are they allowed to fluctuate with some external sources and sinks from/to outside the box, apart from the air-sea CO₂ flux? what happens to all sinking materials ? Even if they stay in the box, we need to know. I ask for major revision, just to be sure this part is improved. For the rest i cannot ask you to re-do everything in 1 or 3D, this will most probably be your next paper anyway.

----------------------------------------------

Also, I would ask for some modifications in the text.

* The English might need some rewording. I am not an english native, so I cannot help much for that, but I would recommend a second read. For example, you make an extensive use of the word “yielding”. It is a nice word, but you should replace some of them with relevant synonyms.

* Page2 line 25 : add some “-” or () : “model – consistent with observations – predicted….”

* P3-L65 : You need to add something here about the reason for 0D.

* P5-table1: You only fill the time resolution information for the wind. Does that mean they’re all the same ? It seems from the text that some data are daily.  You should feel them all, or tell in the table description why the other data have no time resolution information.

* P6-L124 to 126. “ In addition to …”. what you say there sounds obvious, but BGC-model not including mixotroph represent reasonable TA and DIC. Do you have a Reference paper for this statement ?

* P7-equation1 : You should specify that all terms are define in the appendix A.

* P7- equ2 and 3 : You can specify the unit at the end of the equation, and remove the following sentence.

* P7-equ4 : I might be wrong but, shouldn’t the “photo” terms be more like uptake terms ? Phyto absorbs more DIC than the only ones used for the photosynthesis. Isn’t this equation missing the remineralization terms as a source of DIC ?

* P8-equ5 and L179-181 : I don’t understand why you define Aera being negative when the CO₂ flux is toward the sea. In Equ 4, ∂DIC/∂T increases with Aera being positive, what means CO₂ flux toward the sea, and pCO_2,sw > pCO_2,atm. There’s a discrepancy here you might want to correct.

* P11-L237 : “the first three terms of Eq.(10)”, I think you refer to Eq.11, not 10.

* P14-Table3 : just to mention, comparing pH is tricky. Comparing pH change or bias in pH unit can be misleading. Best practice is to compare H⁺ concentration. See Kwiatkowski and Orr, 2018 (https://www.nature.com/articles/s41558-017-0054-0).

* P15 -Fig 5 : the e,f,g and h panels are not useful. There is no additional information, and it’s not even zoomed-in. Instead, I would remove them, make the picture slightly bigger, and highlight the SUP like you do in Fig. 6 with a shading or something similar.

* P16-Fig 6 : The panel d is quite difficult to look out, it can be quite difficult to differentiate the different blue lines (especially nTA and nDIC have very similar colours). Plus, most of the time the curves in this panel are between -100 to +100 µatm, while the y-axe goes from -600 to +600 µatm. Apart from the big events, it’s quite difficult to see what’s happening there. Maybe take the whole page for this picture ?

* P17-L365 to 370 : You forgot to refer to Fig. 6e somewhere in this section.

* P19-L445 : “we could (not cloud) provide”.

* P20 – L477-8 : “While we only considered TA inputs” (only in the allochthonous formulation, I guess), “Rhône river intrusion can also bring nutrient”. This is never explained till now. I already said it in the first part of the review, but you have to be clear about this. The reader cannot fully understand your results otherwise. The model biology only feels the environment changes/variations through the physical forcing only (T, S and light). The biology reacts to the Rhône water only because it is fresher, or to the upwelling because it is colder, but not because of the associated nutrient changes (that do not occur). It is important to tell it because the biology can react in the opposite way than otherwise expected, and explain that because it is a 0D model you probably don’t have much choice (as I understand it). Knowing that, I am surprised by the DIC variations Fig. 3, that are surprisingly good.

* P21-L482-4 : same remark than just above.

Scope of the manuscript, major comments and recommendation

This manuscript is part of a series of papers (Lajaunie-Salla et al, 2021) and manuscripts (Barre et al, 2023, two manuscripts, submitted)

that deal with modelling the ecosystem and biogeochemical fluxes in the Bay of Marseille (BoM). The companion paper (Barre et al, part I)

is centered around the addition of a mixotrophic plankton class into the model that had previously been presented in Lajaunie-Salla et al,

2021, and the consequences thereof. The present manuscript builds on this ecosystem model and studies how well it can capture variations of

the carbonate system and air-sea carbon fluxes. Specifically it asks whether the replacement of a prognostic equation for total alkalinity

(TA) by a diagnostic relationship with salinity can lead to a better representation of the carbonate system. The diagnostic relationship is

chosen to consist of two linear relationships, one with positive S-TA slope (valid for salinities > 27.8), to account for dilution of TA by

freshwater input, and one with negative slope (salinity<37.8), to take into account the high alkalinity of fresh water from the nearby Rhone

river mouth.

The main results of the manuscript are:

• Biological fluxes (uptake of nutrients, nitrification etc.) alone are unable to explain the amplitude of the observed TA variations
• Replacement of the prognostic by a diagnostic equation for TA yields a range of variability in TA that is similar to the observed one
• with this the model is somewhat better able to reproduce some patterns in the variability of pCO2 over the year
• but the model fails to reproduce the observational strong annual air-to-sea CO2 flux in the Bay of Marseille

Although these results indeed show some progress over the previous studies, the paper is well written and the documentation of the model

is detailed and comprehensible, I cannot recommend publication of the manuscript in GMD. This has several reasons:

• Firstly, I would say that it does not fit the scope of GMD, which is there to present new developments in models. While the companion paper, with its presentation of a mixotroph compartment, meets this criterium, the main new thing in this manuscript is a diagnostic relation between TA and salinity, which may improve the results, but conceptually is a fairly small step and has been used in many different models so far. From this side I would rather recommend publication in a different journal, where the focus is more on the considered system itself, i.e. a model for the BoM.
• The switch to a salinity-TA relationship is motivated by the desire to represent the episodic intrusion of freshwater from the nearby Rhone into the BoM, and also the influence of evaporation and precipitation. My first question here is: If these freshwater fluxes affect the TA balance so strongly, should they not also influence DIC and nutrients as well? That this may lead to biases is discussed in lines 477 to 484; but given the extremely high DIC concentration in Rhone water quoted on line 482, I wonder whether this inconsistency may not invalidate the main results.
• This leads me to a more conceptual difficulty with the aproach. The model concept is that of an arbitrary one cubic metre volume at the surface of the bay, and that the model just represents fluxes within  this volume. Spatial fluxes are excluded (except for CO2 flux, more  on that below). This only allows either to model a variable as  purely forced from what is happening inside the box, or to prescribe  it, e.g. as a function of salinity. For a proper modelling of how  external fluxes (e.g. in mol/s) change concentrations (mol/m^3/s)  inside the modeled region, one would have to define the volume that  is affected by these fluxes. A reasonable choice might be to model a  column of water within the mixed layer, as was done in many  zero-dimensional models, e.g. Fasham et al, 1990 or Hurtt and  Armstrong 1999. That would allow a consistent treatment of the  effects of mixing on TA, DIC, nutrients.
• This difficulty becomes especially clear when the authors discuss  the possible reasons for their low net annual air-sea flux of CO2,  which is in contrast to observation-based estimates. Here they state  that "aeration is is simulated by applying Eq (5) to 1 m^3 of  surface water at the SOLEMIO station, which tends to overestimate  the effect of aeration processes on DIC..." (line 534 ff). Indeed:  if the control volume is that shallow, it will be lead to a too fast  approach of DIC towards equilibrium, and hence an underestimate of  fluxes.
• And finally, while the diagnostic TA leads to an improvement in  model results, as evidenced by decreases in %BIAS, RMSD and a 'cost  function' presented in Table 3, the improvements are fairly  modest. Indeed I would be interested in knowing whether a model with  TA prescribed constant at the average of observations would not have  fared at least similarly good as the two presented model cases. I do  think that the results, with their careful discussion of the  influence of upwelling events, the role of temperature, salinity and  ecosystem activity are interesting and can be published, but I'd say   rather in a journal that focusses on coastal processes, not in a   journal with a focus on methodological improvements.

In summary, I think the manuscript can be published after major revision, but I doubt whether GMD is the right journal for it (while it clearly is for the companion paper).

Major comments

Figure 1 is identical to the one on the companion paper and definitively isn't needed should this paper be published in GMD.

State equation for TA, Eq. (1): The terms in the equation are not properly defined. The definition of the terms is given in the Appendix (Table A1), but the table is not referenced here.

The two linear S-TA relations, presented on page 7, which are valid below and above a salinity threshold of 37.8 are discontinuous at S=37.8. This should lead to sudden jumps in the TA value if this threshold is crossed. Are there any effects of this discontinuity visible in the results?

Page 8, line 175: In principle the model equations would not change had you chosen to assume the effected layer to be deeper, except that then the flux would then be distributed over a larger volume. Why not take at least H as the annual average mixed layer depth> Taking it as 1m is equivalent to speeding up the gas exchange by a factor H_real, the real affected layer.

Page 8 and Appendix B, pH and pCO2 calculation: It is good to see that the pH scale differences are taken properly into account, and fugacity has been calculated correctly. But much of this is fairly standard, e.g. the iterative calculation of pH, described in Figure B1. This could be left away.

Page 9, Figure 3: The quality of the Figure is awful. But also it does not convey much information, I would leave it away.

Page 10, lines 218-220: It is not clear to me how the salinity-normalized nTA and nDIC are exactly defined, by a linear correlation with salinity with zweo intercept? If so, why do that if the observed S-TA relation in the oceanographic region is different?

page 11, definition of the statistical indicators: while the definition of RSMD and %BIAS is rather clear, that of the cost function is less clear: Typically, a cost function aggregates model-data-disaggreement for different variables, possibly with different units, into a single scalar variable (Stow et al, 2009). But what exactly the variables are that enter the CF, and how the different variables are nondimenionalized and aggregated into one CF should be properly defined.

page 11, interpretation of statistical indicators: Whether a CF<1 is considered very good, would probably depend on the definition of CF, and cannot be stated as generally as on line 255-256. If the individual cost function terms e.g. consist of the squared model-data difference scaled by the variance in the individual variables, and are then added together, the expected height of the CF would depend on how many different variables are finally added together.  Also, I don't think one can generally say (line 252-253) that a %BIAS<10% is excellent; I would think that depends on the ratio of natural variability to the mean of the variable in question. For TA, with a high background value, a 10% BIAS is rather large.

Figure 4, page 12: The time-series of the difference between the model runs (right panel) does not convey much new information, I would remove them.

Also, I have a question to the data (crosses in Figure 4): to me it is not clear whether all four carbon system variables were measured independently, or whether e.g. DIC and TA were measured, and pH and pCO2 calculated from them. If they were measured independently, how consistent are they with respect top each other, given the used set of carbon system equations?

Table 3, Page 14: If the variance of the observed TA and DIC values is on the order of 20 micromol/kg (note, units should be given in the table), then I'd say a RMSD of about the same order of magnitude is not an excellent agreement. It is not terrible either, though. A similar remark holds for %BIAS.

Figure 5, page 15: The subpanels on the right are simply a cutout of the panels on the left for the summer period. What is the purpose of this duplicated information?

Page 21, Lines 506 ff: Would including the DIC and nutrient input from upwelling improve the model-data agreement, or the converse?

Page 22, Line 530ff: Can one give a conjecture why the model overestimates pCO2 during winter? 

Technical comments:

Line 302: form -> from

Line 483: the overall (decrease) of

References:

Fasham, M.J.R., Ducklow, H.W., McKelvie, S.M., 1990. A nitrogen-based model of plankton dynamics in the oceanic mixed layer. Journal of Marine Research 48, 591-639.

Hurtt, G.C., Armstrong, R.A., 1999. A pelagic ecosystem model calibrated with BATS and OWSI data. Deep-Sea Research I 46, 27-61.

Stow, C.A., Jolliff, J.K., McGillicuddy Jr., D.J., Doney, S.C., Allen, J.I., Friedrichs, M.A.M., Rose, K.A., Wallhead, P., 2009. Skill assessment for coupled biological/physical models of marine systems. Journal of Marine Systems 76, 4–15. https://doi.org/10.1016/j.jmarsys.2008.03.011