Review of Matsumoto et al.: MESMO 3: Flexible phytoplankton stoichiometry and refractory DOM

The manuscript of Matsumoto et al. describes the third version of the MESMO model, a model of intermediate complexity dedicated to studies on ocean biogeochemistry. Major model developments include three phytoplankton functional types, a flexible phytoplankton stoichiometry based on ambient nutrient, temperature and light conditions, as well as an implementation of a refractory DOM pool, in addition to the existing semi-labile pool, also including stoichiometry of C, N, P and Fe. The model is well and comprehensively described and the substantial developments after version 2 fully justify a new version description. The paper is well structured by starting with the developments since version 1, the newly implemented processes and a discussion of the model results. It will be a useful and complete reference for the future use of this model version, and fits perfectly within the scope of the journal. Therefore, I recommend publication after the comments below are addressed:

General comments:

• I think the overall improvement resulting from the changes in the new model version could be made more clear in a quantitative way. This is qualitatively described in paragraph l. 484ff, and on a ‘global-mean’ basis in Table 3. Adding more parameters that are difficult to constrain (such as e.g. the partitioning between DOMsl and DOMr, or uptake stoichiometry) should somehow “pay off” in terms of better agreement with observations or a more flexible representation of processes needed to model biogeochemistry in a changing ocean scenario. It is described that e.g. the nutrient pattern is better in MESMO3 than MESMO2, can you quantify this (with e.g. an RMSE, or a scatter plot/taylor diagram)? Or in Tab. 3, export C:N:P ratio seems to be less close to observations for the new version – what is the cause? Is there a better representation of the spatial pattern instead?
• In combination with the above comment: Could you describe more in detail what your calibration method was (in l. 152f you state that there was a calibration procedure, but you only state the absence of it)? Did you adjust parameters to fit a global mean like in Tab. 3 or a global mean profile of e.g. concentrations, or a specific dataset?
• It would be helpful if some runtime information is given for the new model in comparison to the old one. EMICs are great tools for “cost” efficient calculations, so it would be helpful for future users of the model to be able to calculate runtime.

Specific comments:

• 224: “NPP is produced” – maybe revise, e.g. “all NPP is immediately routed to…”
• 227: do you mean “master nutrient variable”?
• 246: If I understand correctly, your “deep POC split” pathway represents solution of DOM from POM (by e.g. exoenzymes or the like) – is there a reason why this happens only at depth and not at the surface? Or is the relative role of this process in producing DOM larger at the subsurface compared to the surface and therefore it is neglected at the surface?
• 320: Does photodegradation transfer DOMr to DOMsl in your model, or does it just represent a sink for DOMr?
•  463: Can you provide the units for this equation? I'm probably missing sth here, but why is density a factor in the gas-exchange equation  (sticking to SI dimensions L=length, T=time, N=chemical amount/mol, M=mass unit):
  •  [N/(T*L²)] ≠ [M*N/(T*L⁵)] = [L/T] * [M/L³] *   [N/L³] *               [-]
  •  flux                                          k         density     delta conc.         ice fraction
• 489: Based on the changes in the code, can you pinpoint the reason for the improved nutrient fields? A sentence here would be helpful (similar to line 493, where you state why the oxygen minimum depth improved – I find this very helpful).
• 594: Can you give reasons for why the DOM concentration at the surface is so high compared to observations? I was wondering whether this has to do with the way the model in Hansell, 2012 is formulated, that led to estimating the turnover times: There they scale the DOM production with the square root of NPP, this potentially may lead to lower surface concentrations and – by using their turnover time fitted to this model set-up – cause the overestimated concentrations in MESMO3. (eq. 1 in Hansell, D. A., C. A. Carlson, and R. Schlitzer (2012), Net removal of major marine dissolved organic carbon fractions in the subsurface ocean, Global Biogeochem. Cycles, 26, GB1016, doi:10.1029/2011GB004069)

Figures

• Figure 1: The illustration of the sink processes of DOMr might be misleading: does the background decay happen all the time or only at the surface? If it is a true background decay happening everywhere, it might be helpful to write sth like “T_photo + T_background” instead of only “T_photo”. Or does the background decay only happen at a specific depth level – then this would need to be justified why there is no first-order loss process at the surface?
• Figure 3: The caption beginning with “the effect of…” suggests that an effect is shown, i.e. the difference between a simulation with or without RNPG. Please either show an effect, or adjust the caption accordingly.
• Figure 4: Panel b would benefit from a latitudinal mean from observations for comparison, maybe from the dataset you already cite: Martiny, A. C., Pham, C. T., Primeau, F. W., Vrugt, J. A., Moore, J. K., Levin, S. A., & Lomas, M. W. (2013). Strong latitudinal patterns in the elemental ratios of marine plankton and organic matter. Nature Geoscience, 6(4), 279-283.
• Figure 7: Please define DOCt again in the caption

I agree with Reviewer #1’s analysis that this manuscript represents a meaningful contribution to the current and potential users of the MESMO model and is well-suited in scope for GMD. I too recommend publication after my comments below are addressed. Many of my comments ask for better documentation and/or clarity of the methods used while I also offer some recommendations on the future calibration of the optional DOM cycling module.

Equation 3, Line 132-135. It is not immediately apparent how this equation is analogous to the Q10 = 2 relationship. For example, the value of FT for a temperature of 288 K is 0.973 and for a 10 degree higher value of T = 298 K the value of FT is 0.974 or only 0.1% higher.

Line 141 – 144. It is unclear how Eq 1 varies as a function of zml from the statements as written. Does this mean that nutrient uptake can be confined to model level 1 if zml < zc?

Table 2b typo: [NO3]0 is listed as [PO4]0

Equation 5: It is unclear how this equation produces an estimate of P:C (or C:P) that approaches the observed range of ocean C:P stoichiometry of marine phytoplankton and associated nutrient uptake to produce said plankton stoichiometry. For example: if I use Eq. 5 and the parameters from Table 2 to compute the P:C of uptake for the cyanobacteria in a low latitude oligotrophic setting with in situ [PO4] = 0.1 µM, [NO3] = 0.1 µM, T = 298K, and I = 150 W/m2 I get: 6.3 x 0.614 x 0.411 x 0.827 x 0.962 = 1.266 for P:C or 1/P:C = 0.79 for the C:P ratio. Redfield C:P is 106 or 117 to 1 and cyanobacteria should likely exhibit a C:P > 117:1. What am I missing? Is P:C expressed in per mille (‰)? There isn’t any mention of this in the text or Table 2. If that is the case I compute the same example P:C = 0.001266 for C:P = 790:1. Is there any evidence from the literature for such a P-depleted uptake stoichiometry for small phytoplankton? The maximum C:P for cyanos that I am more familiar with maxes out around 350 or 400:1. Figure 4d shows the cyano C:P and its max is ~280:1. Better clarification of how Eq. 5 (and Eq. 6) in conjunction with the parameter values in Table 2 are used to compute PFT C:N:P stoichiometry is needed.

Equation 19 and the discussion on the variability of fDOM (Line 234-243). I refer the authors to Roshan & DeVries 2017 in Nature Communications who report on the latitudinal variability of this parameter from an artificial neural network model of DOC production and export. The variability is within the range ~0.01 – 0.7 which differs from the 0.28 – 0.96 values used in this study. I do not mean to infer that the one study’s range is more correct than the other. I do think that the manner in which the authors vary fDOM as a function of temperature in this study to be a clever formulation and as it captures the likely environmental factors that cause fDOM to vary related to the differences in nutrient levels, light levels, and phytoplankton community structure which co-vary in the ocean.

Equation 31 and value of Kr in Table 2. CaCO3 is more soluble at colder temperatures unlike most solids due to the positive Gibbs free energy change at standard conditions for the dissolution reaction. The parameter Kr is listed as positive 0.69 in Table 2. Thus when computing the temperature scaling term for CaCO3 remin rate exp(Kr*T) you will get a faster remin rate at higher temperatures for CaCO3, not at colder temps. Should the value of Kr be negative 0.69? This would give a faster CaCO3 remin rate at colder temps. CaCO3 solubility is also pH dependent, through the changes of [CO32-] as a function of pH, but without a full carbonate system formulation in MESMO I think it is ok to simply vary CaCO3 remin as a function of temp since pH and temp co-vary in the ocean. The handling of the opal remin rate appears to obey the correct temperature dependence for opal with a higher opal remin rate at warmer temperatures in Eq. 33.

Line 482. Global average fDOM of 0.71 seems problematic. Other estimates that have been constrained by the available DOC concentrations across the global ocean are much lower, e.g. 20% Hansell et al., 2009; 25% Letscher et al., 2015; 17% Romera-Castillo et al., 2016; 20% Roshan & DeVries, 2017. I recommend that the authors revisit the prescribed range on fDOM in future calibration efforts of the DOM module. To be clear, in my opinion I do not believe that a fully calibrated DOM module for MESMO 3 to be a prerequisite for publication of the current manuscript.

Line 486-489 and Figure S1. What happened to the EqPac HNLC in MESMO 3? It’s gone. This suggests there is not enough Fe limitation in that HNLC to keep NO3 and PO4 in the surface or that the new variable P:C and N:C formulations allow for too much NO3 and PO4 drawdown in this HNLC, although this latter mechanism should operate in the Southern Ocean HNLC too which is retained in MESMO 3. This particular deficiency of MESMO 3, it’s lack of an EqPac HNLC should be identified in the text with a possible explanation given by the authors. Also can the authors provide a more formal statistical comparison of MESMO NO3 and PO4 to World Ocean Atlas, e.g. point by point correlation and overall mean bias.

Figure 2. Eukaryotes are limited by silica nearly everywhere in the global ocean except the Southern ocean and subarctic N Pac. How does this compare to the literature? My thoughts are that silica limitation is not as widespread in the real ocean. Perhaps the RNPG formulation for Euk growth is contributing to the high levels of Si limitation? Cyanobacteria are limited by nitrogen almost everywhere except the Southern Ocean. This is likely a result of the biased low surface NO3 concentrations simulated by MESMO 3. Most field estimates suggest the Sargasso Sea is P-limited to all PFTs with this feature only simulated for the diazotrophs in MESMO 3. I recommend a discussion of the PFT nutrient limitation patterns in comparison to the broad patterns identified by the extant literature on the topic, pointing out where MESMO 3 does a better job and where it is deficient.

Figure 4 PFT C:P stoichiometry looks great and largely matches the global patterns and variability seen in field observations.

Figures 7, 8, and 9. Some comments about tuning MESMO DOC. MESMO3 DOC has a positive mean bias when the authors used literature values for the DOC production and remineralization fluxes. The authors then made adjustments to the production flux (from 1% of NPP produces DOCr to 0.2%). I wish to point out to the authors that they can also make adjustments to the remineralization fluxes, i.e. increase the decay rate constants to attempt a better calibration of simulated DOC vs. observed. Neither the DOMr production or consumption fluxes are well constrained observationally at present. The MESMO 3 deep ocean DOCt gradient captures well the thermohaline circulation gradient. The simulation of DOCsl in the upper ocean is reasonable as compared to the neural network data product of Roshan & DeVries. The authors choice of tying the DOMsl production flux to temperature is a good first step to capture the behavior seen in the DOC observational data and the neural network product that DOC production peaks in the subtropical gyres where temperatures are also high. On line 647-650 the authors state that they may try to tie the fDOM parameter negatively to nutrient concentrations to try and achieve a higher DOMsl production flux in the subtropical gyres. In my expert opinion I believe that tying fDOM to both temperature and inversely to nutrients will likely capture the behavior they wish to reproduce. I refer the authors to a review of this topic on ‘extracellular release’ by phytoplankton within Chapter 3 of “Biogeochemistry of Marine Dissolved Organic Matter” 2^nd Edition. Edited by D.A. Hansell & C.A. Carlson published by Academic Press. Rather than temperature, to be more physiologically consistent, the authors may wish to tie fDOM to PAR rather than temperature (although both co-vary in the ocean). This could be left to a future study that focuses on calibrating the DOM module in MESMO 3.

Figure S6. Why not also show the water column depth integrated denitrification flux alongside the N2 fixation flux?