Summary and general appraisal

The manuscript titled “Formulation, optimization and sensitivity of NitrOMZv1.0, a biogeochemical model of the nitrogen cycle in oceanic oxygen minimum zones” by Bianchi et al. presented a new model of the marine nitrogen (N) cycle in and around oxygen minimum zones, which can be incorporated into current ocean biogeochemical models. This model includes the major forms of nitrogen as prognostic tracers and parameterizes their transformations as a function of organic matter availability, substrate concentrations, and the concentration of dissolved O₂. The model is formulated in an idealized 1D representation that allows for formal optimization against in situ observations and sensitivity tests to model parameters and environmental conditions.

The authors found that multiple optimal parameter sets can reproduce observed tracer and rate profiles, given the nonlinearity of the optimization problem, and trade-offs between different model parameters. Despite this variability, the authors found some consistent conclusions:

• Systematic relationships exist between parameters and features of the solutions, with cascading impacts of model parameters and environmental conditions on tracer concentrations in and around the oxygen deficient region;
• Different steps of the denitrification pathway had different O₂ sensitivities, with increasing inhibitions for the different steps, creating an O₂-dependent production window for N₂

The model presented in this manuscript represents a significant improvement over other dynamical N cycle models in terms of the complexity and accuracy of its N cycle representation. This is especially important in oxygen-deficient regions such as the Eastern Tropical South Pacific, where multiple aerobic and anaerobic N cycle metabolisms may overlap. In this regard, the manuscript represents an important and valuable contribution to N cycle modeling.

The strengths of this paper lie in the development of the model itself, the thorough (and very useful) explanation of model equations and parameters, and the explanation of the numerical trade-offs between processes that lead to multiple optimal solutions. In addition, the clever use of an evolutionary optimization algorithm allows for the optimization of the complex and nonlinear cost function. Where the paper falls short is in its discussion of model results. There are interesting features of the solutions and sensitivity tests that bear more discussion than what the authors have provided. For example, the O₂ range defined by K^O2_den2 and K^O2_den3 allows for N₂O production at higher levels of O₂ than previously thought (Dalsgaard et al. 2014). This is an interesting result and should be further discussed. Likewise, increasing the rate of nitrite reduction to N₂O leads to an increase in N* — it is not immediately obvious why this would be the case. These and other results should be discussed in more detail, both with regards to model mechanics and in terms of their biogeochemical relevance.

My other major concern is that the authors do not discuss processes potentially missing from the model until the very last paragraph of the manuscript. In my opinion, this discussion should be moved up to Section 2, “Nitrogen cycle model formulation.” To this list the authors should add dissimilatory nitrate reduction to ammonium (DNRA) (Lam and Kuypers 2011; Kraft et al. 2011), hybrid N₂O production by ammonia-oxidizing archaea (Stieglmeier et al. 2014; Kozlowski et al. 2016; Trimmer et al. 2016), and N₂O production directly from nitrate vie denitrification, i.e. nitrate reduction to nitrite and N₂O entirely within the cell (Ji et al. 2015, 2018). The latter has been shown to occur at rates orders of magnitude higher than nitrite reduction to N₂O in the eastern tropical South Pacific and eastern tropical North Pacific oxygen minimum zones (Ji et al. 2018; Frey et al. 2020) and is probably the most important omission from the model’s N cycle formulation.

It was also somewhat unclear which measurements were used to constrain the model. It seems like some of the model solutions were constrained by rate measurements and others were not; this should be clarified in the text. Additionally, eqns. (5) – (10) in Section 3.1 seem to follow directly from the equations in the previous section. The authors could consider combining these two sections for clarity.

Finally, the authors could spend more time discussing what is potentially lost by modeling bulk reaction rates instead of individual microbial populations. In Section 2.1, “Model Rationale,” the authors begin to address this, but one or two sentences about the trade-offs between the two approaches would be helpful to readers not intimately familiar with either type of modeling.

Specific Comments

Line 33: Given that bacterial ammonia oxidizers oxidize ammonia, not ammonium, it might be more accurate to use the term "ammonia oxidation" instead of “ammonium oxidation.” It remains unclear whether ammonia oxidizing archaea can oxidize ammonium as well as ammonia.

Line 35: Note that this increase in yield is often attributed to a shift from N₂O production as a byproduct of hydroxylamine oxidation to nitrifier-denitrification (Wrage et al. 2001; Stein and Yung 2003).

Lines 88-89: The authors state the advantages of the "system view" approach. But what kinds of information are lost — what are the disadvantages — of modeling chemical reactions instead of microbial populations?

Line 91: An emerging view of N₂O cycling indicates that nitrate reduction to N₂O with no exchange of nitrite with an external pool is a major source of N₂O in and around OMZs (Ji et al. 2018; Casciotti et al. 2018). This important exception to the modularity of the N cycle should be noted.

Line 109: Explain that pathway 2b implicitly represents both N₂O production as a byproduct of hydroxylamine oxidation (Hooper and Terry 1979) and nitrifier-denitrification (Wrage et al. 2001; Stein and Yung 2003).

Line 128: Explain that K_H is a rate constant first-order to POC.

 Line 136/Eqn. (2): It would be helpful to explain why this is better than using a second-order rate constant and making the reaction second order to [X] and [Y].

Line 145: (are first order to the concentration of organic matter)

Line 156: Also cite Frey et al. (2020).

Line 159: Many of these “anaerobic” reactions are shown in this manuscript to occur at non-zero oxygen concentrations. Would it then be more accurate to refer to these reactions as "suboxic" instead of "anaerobic"?

Line 174: Explain that k_v is different from molecular diffusion, which is much slower.

Lines 185-186: What is the residence time of water in the ETSP, and how do the chosen relaxation timescales compare to this residence time?

Line 217/Section 3.3: Consider merging these two sections to make the equations easier to follow.

Line 279: Is the ability of CMA-ES to find a global optimum relevant for a problem that has multiple optimal solutions?

Line 287: Why chose a constant upwelling velocity but a variable vertical diffusion?

Lines 288-289: Explain why targeting the core of the ODZ justifies turning off far-field restoring.

Line 293: Is "anoxic" the right word here if we're talking about water that still has dissolved oxygen (albeit in very small concentrations)? I would suggest "oxygen deficient" as an alternative.

Lines 294-295: Is the mixed layer depth also specified in the model? Because I would expect this to affect the oxygen profile as well.

Lines 308-309: Doesn't this indicate that far-field restoring should be turned on?

Lines 309-310: How does weighting N* instead of nitrate and phosphate individually work around the discrepancies between modeled and observed nitrate and phosphate?

Lines 311-312: What about centering the vertical weight curve around the oxic-anoxic interface, where we would expect the greatest overlap between aerobic and suboxic processes?

Lines 324-325: Is it a problem that these parameters are converging at the boundary?

Line 316: Are there the biogeochemical implications for these trade-offs?

Lines 334-335: This seems somewhat circular. They should show similar profiles, right? If they're being fit to profiles of data that have the same shape?

Lines 338-339: Could some of these discrepancies be attributed to the fact that the rate data are potential rates, especially in the cases where the model produces smaller rates than the measured ones?

Line 339: "Indirect" could potentially be confusing here since rate measurements with ¹⁵N tracer are often referred to as "direct rate measurements."

Lines 358-359: It seems as though this window allows for N₂O production at higher levels of O₂ than 3 µM, as suggested by Dalsgaard et al. (2014) - this is an interesting result.

Line 378-379: What does it mean that most of the optimizations were not able to reproduce the observed N₂O profile, and that the one that does appears to be an edge case with respect to many of the transformation rates in Figure 5 and the proportions of total N loss in Figure 7? Is this evidence for processes missing in the model?

Line 385-386: Is this indicative of an ammonia-oxidizing archaea-dominated regime, as opposed to ammonia-oxidizing bacteria?

Lines 384-395: Again, these trade-offs are discussed primarily in terms of the mechanics of the model, but could they also represent biogeochemical trade-offs? For example, competition between anammox bacteria and nitrite reducing denitrifiers for nitrite?

Line 394: Specify that these half saturation constants are for anammox — it could be confusing since you mention a different set of half saturation constants above.

Lines 421-422: This relates to similar results in Frey et al. (2020).

Line 425/Section 5.4: I know it's outside of the scope of the paper to quantify climate impacts, but is it possible to tie the results in this section to predicted climate-related impacts in the ocean, such as deoxygenation, warming, and increased stratification?

Lines 426-427: What about mixed layer depth?

Lines 429-430: … and so too do the thresholds, themselves.

Line 440: …and yet a small increase in the magnitude of the near-surface N₂O maximum when K_v is increased - why?

Line 452/Section 6: The discussion feels rather short relative to the richness of the results presented above. I would have liked to see how the results (particularly regarding trade-offs and sensitivities to model parameters and environmental variables) relate to recent experimental results.

Line 460: "Capturing the correct underlying dynamics" seems like an overly strong assertion - the underlying dynamics of the model are able to capture observed tracer distributions, but this does not rule out the potential for other formulations of the N cycle to capture tracer distributions equally well.

Line 461: The model contains both aerobic and anaerobic processes. Change this to "the N cycle in and around anaerobic environments."

Lines 463-464: So some of the model solutions are constrained by rate measurements and some are not?

Lines 495-505: Add to this list additional N cycle pathways, such as DNRA, hybrid N₂O production, and N₂O production directly from nitrate.

Figures and tables

Table 1: Explain why some values are N/A. Are these parameters that are optimized for in the model?

Figure 7: It would help to show a panel with the total rate of N loss.

Figure 8: Clarify what "-50m" means.

Figure 9: Why would increasing k_den2 actually increase N*?

Figure 9: There's an interesting asymmetry in the responses of N₂O to changing k_den2 and k_den3. Why is this?

Figure 12: Clarify in the caption that more negative numbers correspond to an increasing POC flux.

Technical corrections

Line 219: typo: "discretized"

Line 298: Typo: "algorithm"

Line 300: Typo: "optimize"

Line 305: Grammar: "from always converging"

References

Casciotti, K. L., M. Forbes, J. Vedamati, B. D. Peters, T. S. Martin, and C. W. Mordy. 2018. Nitrous oxide cycling in the Eastern Tropical South Pacific as inferred from isotopic and isotopomeric data. Deep Sea Res. Part II Top. Stud. Oceanogr. 156: 155–167. doi:10.1016/j.dsr2.2018.07.014

Dalsgaard, T., F. J. Stewart, B. Thamdrup, L. D. Brabandere, N. P. Revsbech, O. Ulloa, D. E. Canfield, and E. F. DeLong. 2014. Oxygen at Nanomolar Levels Reversibly Suppresses Process Rates and Gene Expression in Anammox and Denitrification in the Oxygen Minimum Zone off Northern Chile. mBio 5: e01966-14. doi:10.1128/mBio.01966-14

Frey, C., H. W. Bange, E. P. Achterberg, and others. 2020. Regulation of nitrous oxide production in low-oxygen waters off the coast of Peru. Biogeosciences 17: 2263–2287. doi:doi.org/10.5194/bg-17-2263-2020

Hooper, A. B., and K. R. Terry. 1979. Hydroxylamine oxidoreductase of Nitrosomonas. Production of nitric oxide from hydroxylamine. Biochim. Biophys. Acta 571: 12–20. doi:10.1016/0005-2744(79)90220-1

Ji, Q., A. R. Babbin, A. Jayakumar, S. Oleynik, and B. B. Ward. 2015. Nitrous oxide production by nitrification and denitrification in the Eastern Tropical South Pacific oxygen minimum zone. Geophys. Res. Lett. 42: 10,755-10,764. doi:10.1002/2015GL066853

Ji, Q., E. Buitenhuis, P. Suntharalingam, J. L. Sarmiento, and B. B. Ward. 2018. Global Nitrous Oxide Production Determined by Oxygen Sensitivity of Nitrification and Denitrification. Glob. Biogeochem. Cycles 32: 1790–1802. doi:10.1029/2018GB005887

Kozlowski, J. A., M. Stieglmeier, C. Schleper, M. G. Klotz, and L. Y. Stein. 2016. Pathways and key intermediates required for obligate aerobic ammonia-dependent chemolithotrophy in bacteria and Thaumarchaeota. ISME J. 10: 1836–1845. doi:10.1038/ismej.2016.2

Kraft, B., M. Strous, and H. E. Tegetmeyer. 2011. Microbial nitrate respiration – Genes, enzymes and environmental distribution. J. Biotechnol. 155: 104–117. doi:10.1016/j.jbiotec.2010.12.025

Lam, P., and M. M. M. Kuypers. 2011. Microbial Nitrogen Cycling Processes in Oxygen Minimum Zones. Annu. Rev. Mar. Sci. 3: 317–345. doi:10.1146/annurev-marine-120709-142814

Stein, L. Y., and Y. L. Yung. 2003. Production, Isotopic Composition, and Atmospheric Fate of Biologically Produced Nitrous Oxide. Annu. Rev. Earth Planet. Sci. 31: 329–356. doi:10.1146/annurev.earth.31.110502.080901

Stieglmeier, M., M. Mooshammer, B. Kitzler, W. Wanek, S. Zechmeister-Boltenstern, A. Richter, and C. Schleper. 2014. Aerobic nitrous oxide production through N-nitrosating hybrid formation in ammonia-oxidizing archaea. ISME J. 8: 1135–1146. doi:10.1038/ismej.2013.220

Trimmer, M., P.-M. Chronopoulou, S. T. Maanoja, R. C. Upstill-Goddard, V. Kitidis, and K. J. Purdy. 2016. Nitrous oxide as a function of oxygen and archaeal gene abundance in the North Pacific. Nat. Commun. 7: 13451–13451. doi:10.1038/ncomms13451

Wrage, N., G. L. Velthof, M. L. Van Beusichem, and O. Oenema. 2001. Role of nitrifier denitrification in the production of nitrous oxide. Soil Biol. Biochem. 33: 1723–1732. doi:10.1016/S0038-0717(01)00096-7

Apologies for the late review.

Bianchi et al. present a highly coherent and well written account of a newly developed motel formulation of nitrogen cycling in oxygen deficient zones of the ocean. Their formulation consists of fuctional reactions within the nitrogen cycle, parameterized as a function of the  available carbon pool (POC), michaelis menten saturation sensitivity to substrate pools, and sensitivity to molecular oxygen. Optimizations of a 1D simulation reveals high coherence among runs, wherein salient biogeochemical features of oxygen deficient waters are adroitly represented, inclidng nitrite and nitrous oxide.

The model parameterization is an important step forward in aiding the represenationt of nitrogen in the global ocean, and may enable simulations of the response of different environmental forcings to the produciton of N2O, an important greenhouse gas.

I have no criticisms of the work per se, but do have questions that the authors may want to clarify in their work:

(a) What was the stoichiometric representation of nitrification and anammox? I may have missed it, but did not see it in teh text.

(b) How do the authors justify their respiratory quotients for respiration and denitrification, which are arguably outside empirical limits (those presented in the supplements)?

(c) How sensitive are the simulations to said respiratory quotients?

(d) One limitation  the model partameterizations that I see  is that maximum rates need to be prescribed ("k" values). In reality, maximum rates will depend on the abundance of the functional group of organisms in the water column. This likely explains the high sensitivity of the simulations to prescribed k values, which, in reality, will be depth dependent (e.g, more nitrifiers directly at teh base of the euphotic zone than elsewhere in the water column). While the authors allude to this, I think it merits more discussion.

(e) I'm not sure I understood the sentence at line 340.

On a minor note, the reference list ascribes all findings in OMZs to a select number of recent papers and reviews. I urge the authors to acknowledging original work where pertinent.