Formulation, optimization, and sensitivity of NitrOMZv1.0, a biogeochemical model of the nitrogen cycle in oceanic oxygen minimum zones

. Nitrogen (N) plays a central role in marine biogeo-chemistry by limiting biological productivity in the surface ocean; inﬂuencing the cycles of other nutrients, carbon, and oxygen; and controlling oceanic emissions of nitrous oxide (N 2 O) to the atmosphere. Multiple chemical forms of N are linked together in a dynamic N cycle that is especially active in oxygen minimum zones (OMZs), where high organic matter remineralization and low oxygen concentrations fuel aerobic and anaerobic N transformations. Biogeochemical models used to understand the oceanic N cycle and project its change often employ simple parameterizations of the network of N transformations and omit key intermediary tracers such as nitrite (NO − 2 ) and N 2 O. Here we present a new model of the oceanic N cycle ( Nitr ogen cycling in O xygen M inimum Z ones, or NitrOMZ) that resolves N transformation occurring within OMZs and their sensitivity to environmental drivers. The model is designed to be easily coupled to current ocean biogeochemical models by representing the major forms of N as prognostic tracers and parameterizing their transformations as a function of seawater chemistry and organic matter remineralization, with minimal interference in other elemental cycles. We describe the model rationale, formulation, and numerical implementation in a one-dimensional representation of the water column that reproduces typical OMZ conditions. We further detail the optimization of uncertain model parameters against observations from the eastern tropical South Paciﬁc OMZ and evaluate the model’s ability to reproduce observed proﬁles of N tracers and transformation rates in this region. We conclude by describing the model’s sensitivity to parameter choices and environmental factors and discussing the model’s suitability for ocean biogeochemical studies.


Introduction
Nitrogen (N) limits phytoplankton production over large swathes of the ocean (Moore et al., 2013).Most of the N in the ocean is present as dissolved dinitrogen gas (N 2 ); however, only fixed N, e.g., ammonium (NH + 4 ) and nitrate (NO − 3 ), can be readily utilized by planktonic microorganisms, with the exception of N-fixing diazotrophs (Capone et al., 2008).The inventory and chemical form of N in the ocean are controlled by an active nitrogen cycle, whereby different chemical forms of the element are utilized as substrates for growth by a variety of microorganisms, either to supply building blocks for organic molecules or to fuel metabolism via redox reactions (Capone et al., 2008;Kuypers et al., 2018).As a result, the residence time of fixed N in the ocean is on the order of 3000 years or less, about 1 order of magnitude shorter than for the macronutrient phosphorous (Gruber and Galloway, 2008;Wang et al., 2019).
The ocean's inventory of fixed N is dominated by NO − 3 , and the main N cycle reactions consist of uptake and assimilatory reduction of NO − 3 to NH + 4 (here used interchangeably with ammonia, NH 3 ) and the oxidation of NH + 4 back to NO − 3 following the decomposition of organic matter and nitrification (Fig. 1).Only when the concentration of dissolved oxygen (O 2 ) drops to suboxic or anoxic levels (typically below 5 mmol m −3 ) do additional metabolic pathways involving N become relevant, as observed in the ocean's oxygen minimum zones (OMZs) and low-O 2 sediments (Lam and Kuypers, 2011).These reactions include the three main steps of heterotrophic denitrification, i.e., the oxidation of organic carbon (OrgC) with NO − 3 , nitrite (NO − 2 ), and nitrous oxide (N 2 O), and anammox, the chemolithotrophic oxidation of NH + 4 with NO − 2 .Both denitrification and anammox lead to the production of N 2 and thus remove fixed N from the ocean (Bianchi et al., 2012;DeVries et al., 2012DeVries et al., , 2013)).Ammonia oxidation is another source of N 2 O -a powerful greenhouse gas and a leading agent of ozone destruction in the stratosphere.The number of N 2 O molecules produced per NH 3 oxidized, i.e., the yield of this reaction, increases as O 2 declines (Goreau et al., 1980;Nevison et al., 2003), likely caused by a shift from N 2 O production as a byproduct of hydroxylamine oxidation to nitrifier denitrification (Hooper and Terry, 1979;Wrage et al., 2001;Stein and Yung, 2003).Because of denitrification and enhanced production by ammonia oxidation, OMZs are important sources of N 2 O to the atmosphere (Naqvi et al., 2010;Yang et al., 2020), with the largest emissions observed right above shallow oxygendeficient waters (Arévalo-Martínez et al., 2016).
The emerging picture of the ocean's N cycle is that of a web of inter-dependent transformations that is particularly active in OMZs, where overlapping aerobic and anaerobic re-actions exchange nitrogen metabolites and substrates (Lam and Kuypers, 2011;Kuypers et al., 2018), ultimately controlling fixed nitrogen removal and N 2 O production.While there is evidence that organic matter and O 2 regulate the rates and relative importance of N transformations (Babbin et al., 2014;Dalsgaard et al., 2014), our mechanistic understanding of these environmental controls against the backdrop of oceanic variability remains limited.Ocean biogeochemical models can shed light on the expression of the N cycle reactions in a dynamic environment.These models have included N as a macronutrient since the beginning, representing NO − 3 and NH + 4 assimilation by phytoplankton and subsequent nitrification (Fasham et al., 1990;Sarmiento et al., 1993;Moore et al., 2004).With the advent of more complex earth system models, biogeochemical representations have progressively expanded to include more detailed representations of the N cycle, including N fixation, denitrification, and N 2 O production (Aumont et al., 2015;Séférian et al., 2020;Stock et al., 2020;Long et al., 2021b).
The representation of N transformations in models often relies on crude assumptions that simplify the network of N reactions and their controls to simple empirical parameterizations.For example, models that include N 2 O cycling often rely on parameterizations that link N 2 O production to nitrification or aerobic respiration (Suntharalingam and Sarmiento, 2000;Nevison et al., 2003;Manizza et al., 2012;Jin and Gruber, 2003), overlooking N 2 O sources and sinks by denitrification.These models also conflate anammox and denitrification into a single N 2 production term.Explicit cycling of NO − 2 under low O 2 , with the observed co-occurrence of NO − 2 production from NO − 3 dissimilatory reactions; reduction to N 2 O and N 2 by denitrification and anammox; and reoxidation to NO − 3 are missing (Lam and Kuypers, 2011;Kalvelage et al., 2013;Babbin et al., 2014Babbin et al., , 2015;;Buchwald et al., 2015a;Babbin et al., 2017).
The goal of this paper is to present a new model of the oceanic N cycle designed to be incorporated in current ocean biogeochemical models, with a particular focus on processes occurring within OMZs.We refer to this model as NitrOMZ (Nitrogen cycling in Oxygen Minimum Zones).The model explicitly represents the major forms of N found in seawater as prognostic tracers and parameterizes the transformations that connect them as a function of seawater chemistry.This formulation is informed by recent observations that describe the response of N cycle reactions to environmental controls, in particular the availability of substrates and dissolved O 2 .We detail the implementation of the model in an idealized one-dimensional representation of the water column that allows for comparison to in situ observations, formal optimization, and studies of the model sensitivity to parameter choices and environmental conditions.
The rest of the paper is organized as follows: Sect. 2 discusses the rationale and formulation of the model, Sect. 3 the implementation of the model, Sect. 4 the model optimization against tracer and rate observations, Sect. 5 the performance of the model and its sensitivity to environmental parameters, and Sect.6 the implications and conclusions of the work.

Model rationale
The NitrOMZ model is based on the current understanding of the N cycle in OMZs (Lam and Kuypers, 2011;Kuypers et al., 2018) as mediated by six major species: N 2 , NO − 3 , NO − 2 , N 2 O, NH + 4 , and organic nitrogen (OrgN) in either dissolved or particulate form.We only explicitly model NH + 4 (the dominant dissolved form) and do not distinguish it from NH 3 .We also assume that organic nitrogen is linked to organic carbon by fixed stoichiometry (Anderson and Sarmiento, 1994), although variable stoichiometry can easily be accommodated.
A schematic of the model's tracers and transformation is shown in Fig. 1.Our approach represents a natural progression for current biogeochemical ocean models and takes a "system view" of the N cycle by focusing on the biogeochemistry of N transformation reactions (Lam and Kuypers, 2011) rather than on microbial ecology (Penn et al., 2016;Louca et al., 2016;Zakem et al., 2018;Penn et al., 2019).That is, the model explicitly resolves N chemical tracers and their transformations but not the populations of microbes that are responsible for these reactions.
The underlying assumption is that the occurrence and rates of N transformations are controlled by, and can be predicted from, the physical and chemical conditions of the oceanic environment.Implicitly, the model assumes that diverse populations of microbes are always present in the water column and that their activity (i.e., metabolic rate) is controlled by the abundance of substrates, analogous to chemical reactions, and dissolved O 2 , which inhibits anaerobic reactions (Kalvelage et al., 2011;Babbin et al., 2014;Dalsgaard et al., 2014;Ji et al., 2018a;Sun et al., 2021b).The focus on dissolved N forms and reaction rates bypasses poorly known aspects of microbial population dynamics, which are topics of ongoing research (Louca et al., 2016;Zakem et al., 2018;Penn et al., 2019).
We assume that each reaction is implicitly mediated by specialized microorganism groups, each relying on a distinct metabolism (Lam and Kuypers, 2011;Kuypers et al., 2018).Thus, the model represents a "modular" N cycle, with individual reaction steps (i.e., individual redox reactions) rep-resented separately and connected by the exchange of dissolved substrates (Graf et al., 2014;Kuypers et al., 2018).This premise is grounded on observations of high specialization and streamlined genomes for marine prokaryotes (Giovannoni et al., 2014), including microorganisms carrying genes for N-based metabolic reactions (Ganesh et al., 2015;Kuypers et al., 2018).
These assumptions are sufficient for providing a broad representation of microbial N transformations and their environmental expressions in the ocean, while limiting model complexity and the proliferation of poorly constrained parameters.They are also grounding steps toward models that explicitly represent microbial populations, including their diversity and dynamics in OMZs (Louca et al., 2016;Penn et al., 2016;Zakem et al., 2018;Penn et al., 2019).

Model tracers and processes
The model focuses on microbial processes that take place below the euphotic zone, as driven by the flux of organic matter produced near the surface and exported into the ocean interior by the biological pump (Boyd et al., 2019).We include heterotrophic and chemolithotrophic pathways that are commonly observed in the open ocean and require N species as substrates (Kuypers et al., 2018) (Fig. 1).Additional pathways, for example, involving sulfur or iron, could also be represented following a similar approach.
Heterotrophic reactions resolved by the model (Fig. 1) consist of aerobic organic matter respiration (R rem , pathway 1), which relies on O 2 as the oxidant, and the three main steps of denitrification: dissimilatory NO − 3 reduction to NO − 2 (R den1 , pathway 4), NO − 2 reduction to N 2 O (R den2 , pathway 5), and N 2 O reduction to N 2 (R den3 , pathway 6).Chemolithotrophic processes consist of aerobic oxidation of NH + 4 to both NO − 2 (R no2 ao , pathway 2a) and N 2 O (R n2o ao , pathway 2b, via both hydroxylamine oxidation and nitrifier denitrification); aerobic oxidation of NO − 2 to NO − 3 (R no , pathway 3); and anammox, the anaerobic oxidation of NH + 4 with NO − 2 to produce N 2 gas (R ax , pathway 7).Reactions are parameterized as functions of substrates (i.e., model tracer concentrations) and environmental parameters such as dissolved O 2 and organic matter.Tracers are expressed as concentrations, with units of mmol m −3 .We do not include an explicit representation of nitric oxide, NO, because of the poor understanding of its cycle in the marine environment (Ward and Zafiriou, 1988).NO is thought to be an obligate intermediate or a byproduct of N cycle reactions, including nitrification and denitrification (Schreiber et al., 2012).However, it is a very reactive chemical with extremely low concentrations (on the order of pmol m −3 ) and rapid turnover in seawater (Ward and Zafiriou, 1988).As a consequence, in situ NO observations are limited (Lutterbeck et al., 2018), and rate measurements targeting NO reactions are missing.Implicitly, we assume that NO cycles so rapidly that accumulation and transport by https://doi.org/10.5194/gmd-16-3581-2023 Geosci.Model Dev., 16,2023 the oceanic circulation are negligible and that its dynamics can be folded into the cycle of other N tracers.
There are also several notable processes that are not represented in the current model formulation but could be introduced in future releases.Some of these processes (e.g., dissimilatory NO − 2 reduction to NH + 4 , DNRA) are not thought to be quantitatively relevant in oceanic oxygen minimum zones.Others, while relevant, require further measurements to constrain their significance and response to environmental variability.
Production of N 2 O via NH + 4 oxidation in NitrOMZ is represented as a single O 2 -dependent function designed to model the transition in bacterial metabolisms from predominantly hydroxylamine oxidation to nitrifier denitrification at low O 2 (Hooper and Terry, 1979;Wrage et al., 2001;Stein and Yung, 2003;Nevison et al., 2003).However, growing evidence suggests that ammonia-oxidizing archaea (AOA, which greatly outnumber their bacterial counterparts) can also produce N 2 O via a hybrid mechanism (Santoro et al., 2011;Löscher et al., 2012).Production of N 2 O via AOA appears to be similarly enhanced at low O 2 (Trimmer et al., 2016;Santoro et al., 2021), although evidence from Stieglmeier et al. (2014) argues otherwise.
DNRA, which can be dominant in anoxic sediment, has been sporadically observed in the water column of oxygendeficient zones, where it may provide an additional source of NH + 4 to anammox bacteria (Lam et al., 2009;Lam and Kuypers, 2011;Kraft et al., 2011;Jensen et al., 2011).However, DNRA is commonly undetectable in OMZ waters (Kalvelage et al., 2013;De Brabandere et al., 2014), and its importance to the N cycle of OMZ is still debated (Long et al., 2021a).
Recent tracer incubation studies show substantial and often dominant formation of N 2 O from NO − 3 rather than NO − 2 (Ji et al., 2018b;Frey et al., 2020).This suggests that denitrifying bacteria capable of direct production of N 2 O from NO − 3 reduction (as NO − 2 reduction proceeds entirely within the cell) could be a major source of N 2 O.This idea, which contrasts with the model assumption of a fully modular N cycle, is further supported by isotopic evidence (Casciotti et al., 2018).However, observations needed to constrain the proportion of N 2 O from NO − 3 and NO − 2 and its environmental sensitivity remain limited (Ji et al., 2018b;Frey et al., 2020).
Other work suggests the occurrence of NO − 2 oxidation in apparently O 2 -deficient waters (Buchwald et al., 2015b;Babbin et al., 2020;Sun et al., 2021a).This may involve NO − 2 oxidation coupled to iodate reduction or NO − 2 disproportionation -two poorly characterized processes.It may also reflect the high affinity to O 2 of nitrite-oxidizing bacteria (Bristow et al., 2016) in regions where vanishing O 2 concentrations are maintained by infrequent lateral intrusions (Buchanan et al., 2023).
Finally, the model could easily accommodate missing processes that couple the N cycle with other elemental cycles, in particular carbon and sulfur.These include formation of or-ganic matter by chemolithotrophy; changes in inorganic carbon chemistry (e.g., pH) by anaerobic reactions (Cinay et al., 2022); and additional metabolic pathways such as anaerobic oxidation of sulfide with NO − 3 (Callbeck et al., 2021) and anaerobic oxidation of methane with NO − 2 (Thamdrup et al., 2019), both chemolithotrophic denitrification reactions.

Model equations
Heterotrophic reactions (i.e., organic matter remineralization) are parameterized as a function of the respective oxidants and organic matter concentration and expressed in carbon units per unit volume and time.Heterotrophic reaction rates are assumed to be on the first order in the concentration of organic matter and limited by the oxidant, following a Michaelis-Menten formulation (Johnson and Goody, 2011).Anaerobic reactions are inhibited by the presence of O 2 , based on an exponential limitation term (Dalsgaard et al., 2014).The resulting equation for a general heterotrophic reaction is (1) Here, H indicates the heterotrophic process considered (e.g., dissimilatory reduction of NO − 3 to NO − 2 ), R H the heterotrophic reaction rate (mmol C m −3 s −1 ), k H the specific first-order reaction rate (s −1 ), [X] the concentration of the oxidant (i.e., O 2 , NO − 3 , NO − 2 , or N 2 O), K X H the halfsaturation constant for oxidant uptake (mmol m −3 ), K o 2 H the scale for inhibition of the reaction by O 2 (mmol m −3 ), and POC the concentration of particulate organic matter in units of mmol C m −3 .No O 2 inhibition is applied to aerobic respiration (i.e., K o 2 H can be thought of as arbitrarily large).Chemolithotrophic reactions are proportional to the respective substrates.A maximum reaction rate is modulated by the concentration of oxidants and reductants, following Michelis-Menten dynamics.For anaerobic reactions (here, anammox), an O 2 -dependent inhibition term limits the reactions when O 2 is present.The resulting equation for a general chemolithotrophic reaction is Here, A indicates the chemolithotrophic process considered (e.g., anammox); R A the reaction rate (mmol N m −3 s −1 ); k A the maximum reaction rate when the process is not limited (mmol N m −3 s −1 ); [X] and [Y ] the concentrations of the oxidant and reductant, respectively (e.g., NO − 2 and NH + 4 for anammox); K X A and K Y A the half-saturation constants for oxidant and reductant uptake, respectively (mmol m −3 ); and K o 2 A the scale for inhibition of the reaction by O 2 (mmol m −3 ).For aerobic reactions, K o 2 A is set to infinite, removing O 2 inhibition.
Equations for each of the heterotrophic and chemolithotrophic reactions are presented in Appendix A1 and A2, respectively; parameter names, units, and suggested values from the literature are presented in Table 1.

Model assumptions and parameterizations
In the model, we assume that heterotrophic reactions are firstorder to the concentration of organic matter; thus all organic matter can be utilized by microorganisms without saturation at high concentrations.Because of the low abundance of organic matter in seawater and extensive colonization of particles by heterotrophic bacteria, this is a reasonable first-order assumption.However, see Nguyen et al. (2022) for a discussion of microbial-particle interactions in ocean biogeochemical models and more complex aspects of their dynamics.For simplicity, we represent organic carbon by a single component.This assumption is easily relaxed to include multiple carbon species, for example, separate particulate or dissolved forms.
We do not explicitly model conversion of dissolved CO 2 to organic matter by chemolithotrophy because of the small rates compared to the remineralization of organic matter in the upper ocean.This assumption can also be relaxed in future implementations of the model, allowing a more complete integration between chemolithotrophy and the carbon cycle.
The use of an exponential inhibition term for anaerobic reactions by O 2 is based on the observation that they are limited at O 2 concentrations to a few mmol m −3 or smaller (Dalsgaard et al., 2014;Babbin et al., 2015;Frey et al., 2020).However, coexistence of anaerobic and aerobic reactions at O 2 concentrations of 10-20 mmol m −3 or higher is also observed (Kalvelage et al., 2011), perhaps related to the presence of redox microenvironments within organic particles (Bianchi et al., 2018;Smriga et al., 2021), which are not explicitly considered here.The exponential inhibition formulation has the advantage of being controlled by a single parameter, allows anaerobic reactions at concentrations of finite O 2 , and approximates empirical rates from incubation experiments reasonably well (Dalsgaard et al., 2014).
Parameter values for maximum reaction rates, halfsaturation constants, and O 2 inhibition terms (Eqs. 1 and 2) are informed by analysis of previous work and further optimized against in situ observations of tracers and rates (Sect.4).Table 1 presents a list of the model parameters and measured values based on a review of the literature.Note that these studies are based on shipboard and laboratory incubations that differ in the setup, conditions, and microbial populations tested.Despite these caveats, experimental results provide valuable starting points to further constrain parameter values in the model.

One-dimensional model setup
We implement the model for a one-dimensional water column that includes physical transport by vertical advection and turbulent diffusion (Wyrtki, 1962) and, if required, parameterized lateral transport by horizontal currents and eddies (Gnanadesikan et al., 2013;Bettencourt et al., 2015).The model is configured to represent the typical weak upwelling conditions that characterize open ocean oxygen minimum zones, following previous work (Babbin et al., 2015).
In the one-dimensional framework, the conservation equation for the concentration [C] of a generic dissolved tracer can be written as where w u is the vertical upwelling velocity (m s −1 ) and K v is the vertical turbulent diffusion coefficient (m 2 s −1 , distinct from molecular diffusion, which is much smaller), both of which can be a function of depth.The first and second summations are, respectively, over the N H heterotrophic and N A chemolithotrophic processes that involve the tracer (Eqs. 1 and 2), with r i C,H and r i C,A being the corresponding stoichiometric ratios (Appendix A4).LT represents any parameterized lateral transport process.The explicit equations for each of the model tracers are detailed in Appendix A5.
The lateral transport term LT can be included to parameterize horizontal circulation by advection and diffusion in the one-dimensional framework.Typically, these terms are simplified by a linear restoring to far-field tracer concentration profiles (Babbin et al., 2015), [C] far , with a relaxation timescale τ C (s): For typical open ocean conditions, τ C can be estimated as the minimum of an advective timescale L U and a diffusive timescale, L 2 K H , where L, U , and K H are, respectively, the horizontal spatial scale, the horizontal velocity scale, and the horizontal eddy diffusion.Assuming L on the order of 1000 km, U on the order of 0.01 m s −1 , and K H on the order of 1000 m 2 s −1 results in a timescale τ C = 10 8 s, i.e., on the order of 3 years and in agreement with recent estimates of the residence time of water within the eastern tropical South Pacific (ETSP) (Ji et al., 2015b;Johnston et al., 2014).

Organic matter remineralization
In the one-dimensional model implementation, we represent organic matter (OrgC and OrgN) as a single particulate orhttps://doi.org/10.5194/gmd-16-3581-2023 Geosci.Model Dev., 16, 3581-3609, 2023   ganic carbon (POC) class that sinks through the water column.We assume that this sinking is rapid compared to advection and diffusion, leading to a steady-state distribution of POC that is only controlled by sinking and remineralization (Kriest and Oschlies, 2008).Since remineralization rates are proportional to the concentration of organic matter, the resulting steady-state one-dimensional equation for POC is where w s is the depth-dependent sinking speed of POC in the water column, and k eff,i H (s −1 ) is the effective rate constants for each heterotrophic process, i.e., the maximum rate constants multiplied by the respective substrate limitation and O 2 inhibition terms (Eq.1).
Considering the flux of sinking POC, POC (mmol C m −2 s −1 ), POC = w s • POC. ( 6) Equation ( 5) can be written as or, equivalently, Equation ( 7) can be recast to relate the concentration of POC in the water column to the remineralization of the POC flux with depth: The advantage of Eq. ( 9) is that it allows us to diagnose sinking POC concentrations when the POC flux and remineralization rate constants are known.In the one-dimensional implementation of the model, we parameterize the POC flux following a typical depth-dependent power-law function or Martin curve (Martin et al., 1987;Berelson, 2001;Primeau, 2006): where z 0 is the upper boundary of the model and b the powerlaw or Martin coefficient.A plot of the model POC is shown in Fig. C1.Another advantage of this formulation is that it allows coupling NitrOMZ to more complex parameterizations for the remineralization of organic matter in ocean biogeochemical models, some of which rely on explicit representation of sinking organic particles and some of which only represent sinking organic particle fluxes in the water column (Moore et al., 2004;Sarmiento et al., 2010;Aumont et al., 2015;Stock et al., 2020;Long et al., 2021b).Because NitrOMZ's equation can be cast as a function of prescribed vertical organic matter flux or remineralization profiles, the model can be coupled to existing biogeochemical models with minimal interference in their formulation of organic matter cycles.

Numerical implementation of the one-dimensional model
For the purpose of testing and illustration, we implement Ni-trOMZ in a one-dimensional representation of the water col-umn below the mixed layer, following previous work (Babbin et al., 2015).Model tracers are discretized on a onedimensional vertical grid, with equal spacing z = 10 m, where z is depth.Boundary conditions are set at the top (z 0 ) and bottom grid (z bot ) cells, as Dirichlet (or fixed concentration) boundary conditions, with values taken from observations (Tables B2-B3).The conservation equation for each tracer (following Eq. 3; see Appendix A5 for full equations) is then solved using a forward-in-time, centered-in-space numerical scheme, with a constant vertical grid spacing, and the option for a variable or constant time step.In the baseline simulations (Fig. 2), we adopt a time step of 5 d for the initial 650-year spinup and decrease it to 3 h for the final 2 years of the simulation (years 698 and 699) to increase accuracy.As in Babbin et al. (2015), NitrOMZ does not represent primary production in the surface layer and is instead forced at the uppermost boundary by a flux of sinking POC, POC (z 0 ) = w s (z 0 )•POC (z 0 ), where POC (z 0 ) provides the boundary condition for POC.The flux POC remineralizes in the water column based on a Martin curve profile (Eq.10).At each depth, the steady-state conservation equation for POC (Eq.8) is solved with a forward-in-space method, using a depth-dependent sinking speed w s chosen to produce, together with the maximum aerobic remineralization rate constant, k rem , a POC flux profile matching a Martin curve with exponent b appropriate for the oxygenated ocean (Primeau, 2006;Weber and Bianchi, 2020).To this end, the sinking speed is calculated at each depth as The concentration of POC in the water column is then diagnosed using Eq. ( 9) and used to calculate the heterotrophic remineralization rates R H in Eq. (1) (see Appendix A1).
Under constant forcings and boundary conditions, the model tracers evolve towards steady state ( ∂ [C]  ∂t ≈ 0, Fig. 2) with a timescale τ SS that can be estimated from the advection velocity w u , the turbulent vertical diffusion K v , and the vertical scale H as the minimum between H w u and H 2 K v .For w u on the order of 10 m y −1 , K v on the order of 10 −5 m 2 s −1 , and a vertical scale of 1000 m, the timescale to approach steady state is τ SS = 3 × 10 10 s or about 100 years.
Figure 2 shows an example of model spinup to steady state in NitrOMZ, with parameters taken from an optimal solution discussed in Sect.5.2 and uniform initial tracer concentrations in the water column.At the start of the simulation, high water column O 2 leads the aerobic remineralization (R rem ) to dominate total POC consumption.As the simulation proceeds, an O 2 minimum develops in subsurface waters, reaching suboxic (< 10 mmol O 2 ) concentrations around year 100.NO − 3 reduction rates (R den1 ) are relieved of O 2 inhibition and begin to take up a larger fraction of total POC remineralization, as revealed by the depletion of N

Model optimization strategy
The model contains 23 major parameters that control the N cycle, some of which are relatively well constrained by observations, whereas others are poorly known and can plausibly span a broad range of values (Table 1).In the model, these parameters approximate complex or poorly known aspects of microbial physiology, metabolism, and ecology and are thus intrinsically uncertain.In order to select a set of parameters that produces a realistic representation of the N cycle in OMZ, we adopt a "metaheuristic" approach based on application of an optimization algorithm, following an established strategy in ocean biogeochemistry (Schartau and Oschlies, 2003;Ward et al., 2010;Kriest et al., 2017).
To conduct this optimization, we compile available tracer and biogeochemical rate observations for the ETSP OMZ from a July 2013 cruise aboard the R/V Nathaniel B. Palmer, for which abundant trace and rate measurements are available (Fig. 5) (Ji et al., 2015b;Peng et al., 2016;Babbin et al., 2017Babbin et al., , 2020)), as well as from other cruises in the region (Kalvelage et al., 2013).The observations are then used to define a cost function based on normalized squared deviations between model profiles and observations.The cost function is minimized by applying a covariance matrix adaptation evolutionary strategy algorithm (CMA-ES; discussed in Sect.4.1), which finds a local optimal solution in the model's multi-dimensional parameter landscape.
The optimization is characterized by large dimensionality, strong non-linearity, a significant computational cost (requiring several 10 000 s model runs to converge), and inherent flexibility in the formulation of the cost function (Schartau and Oschlies, 2003;Kriest et al., 2017).Thus, instead of seeking a single global optimal solution, we generate an ensemble of optimal solutions that provide equally acceptable representations of OMZ processes based on the cost function.
To this end, we apply the optimization multiple times, varying the formulation of the cost function slightly and assigning a random error to the observations for each optimization (Table B4).As a result, we produce a set of equally plausible optimal solutions that we further evaluate to select a final parameter set based on additional comparisons with observations, which we use for further analysis.

Optimization algorithm
The CMA-ES is a stochastic, population-based algorithm that seeks to minimize an objective cost function (Hansen et al., 2009).The CMA-ES falls within the broader class of evolutionary optimization algorithms, where the search for an optimal solution proceeds by an iterative improvement of a population of parameters, with each iteration including a stochastic "evolutionary" element, in loose analogy with biological processes of mutation, recombination, and selection (illustrated in Fig. 3).In contrast with typical evolutionary computation algorithms such as genetic algorithms, in the CMA-ES the mutation and recombination operations are substituted by sampling from a multivariate normal distribution in which parameters (the covariance matrix) are deterministically updated based on previous iteration steps (Hansen, 2006).
The CMA-ES has been shown to be more efficient (i.e., requiring fewer objective function evaluations), accurate (i.e., able to approximate the global optimum when it is known to exist), and robust (i.e., not overly sensitive to the initial choice of parameters) compared to other optimization algorithms, when applied to multi-dimensional, non-linear optimization problems (Hansen et al., 2009;Hansen, 2023).These properties make it suitable for optimization of ocean biogeochemical models (Kriest et al., 2017).A detailed description of the algorithm procedure can be found in Hansen (2023); an overview of the main steps of the algorithm and its application to ocean biogeochemistry are presented in Kriest et al. (2017).

Optimization implementation
As an illustration of NitrOMZ, we perform a series of optimizations against ETSP OMZ observations.For this configuration, we set a constant upwelling velocity (w up ) but impose a variable vertical diffusion (K v ) profile, with lower diffusion in upper stratified layers, and a transition to higher diffusion in deeper layers (Fischer et al., 2013) (Fig. C1, left panel).This is a simplifying assumption that allows us to control the vertical scale for advective-diffusive transport (given by the ratio between vertical diffusivity and upwelling velocity, K v w up ), without requiring vertical divergence terms in the conservation equation for tracers associated with variable vertical velocities.Since this simulation targets the core of the OMZ, generally characterized by sluggish horizontal circulation (Karstensen et al., 2008), we turn off far-field tracer restoring.This simplifies analysis of model balances between transport and reaction rates, while resulting in realistic tracer distributions.The top and bottom boundary conditions are listed in Table B3 and are extracted from observations.
As a first step, we select parameters that control aerobic remineralization processes (R rem ) and lead to a realistic vertical O 2 profile relative to ETSP observations, including the vertical position and thickness of oxygen-deficient waters  (O 2 < 5 mmol m −3 ) (Fig. 5).These consist of the vertical diffusion and upwelling magnitude, the Martin curve coefficient (b), and the upper-ocean POC flux ( top poc ), based on values consistent with observations (Table B2 and Fig. C1).For simplicity, we also set the maximum aerobic remineralization rate (k rem ) and the O 2 half-saturation constants for NH + 4 and NO − 2 oxidation (K o 2 ao and K o 2 no , respectively) to reported values in the literature (see Table 1).We then employ the CMA-ES algorithm in NitrOMZ to optimize the remaining 20 parameters that control heterotrophic and chemolithotrophic reactions in Fig. 1, using the range of parameter values listed in Table B1.
To optimize more uncertain parameters that control the anaerobic N cycle, we then conduct four sets of optimizations, with cost functions devised to match desired characteristics of tracer and rate profiles in the ETSP OMZ.Briefly, the cost function is calculated as the mean square of the difference between observations and model output profiles for a series of variables that include tracers and N transformation rates (listed in Table B4).Before each optimization, a random error of up to 20 % is assigned to each observation to increase the variability in observational constraints and improve the robustness of the optimization ensemble by preventing it from always converging in the neighborhood https://doi.org/10.5194/gmd-16-3581-2023 Geosci.Model Dev., 16, 3581-3609, 2023 of a specific local minimum controlled by non-relevant features of the observations.Three additional constraints are imposed to improve the fit to observations for N cycle processes occurring within the core of the OMZ.First, all rates are weighted equally, whereas different weights are assigned to each tracer, giving higher weight to N 2 O and NO − 2 , which are central to the anaerobic N cycle.Because of possible influence from horizontal advection in observations, discrepancies exist between modeled and observed NO − 3 and PO 3− 4 .To compensate for this, we also assign lower weights to NO − 3 and PO 3− 4 and higher weight to N * .Second, a depthdependent weighting scheme is included to emphasize the match to observations in the OMZ interior.This vertical weight is shaped as a Gaussian curve centered at the core of the observed OMZ, where the bulk of anaerobic transformations targeted by our model occurs so that values within the core of the OMZ are weighted up to twice as much as values outside the OMZ.Finally, N cycle transformation rates are shifted vertically to match their depth relative to the oxycline (here defined as O 2 = 1 mmol m −3 ) in both model and observations and rescaled by a factor proportional to observed vs. modeled POC flux in the upper ocean.The only difference between the four sets of optimization is the relative weights assigned to each tracer, listed in Table B4.In total, we obtain 382 optimized parameter sets for further analysis.

Optimization results
The distributions of the parameter values from the 382 sets of optimizations (see Sect. 4.2 and Table B4) are shown in Fig. C2.Rather than always converging to the same set of parameters, the optimization shows some variability for specific parameters.This reflects the stochastic nature of the CMA-ES algorithm, the inclusion of random variations in the observations, and the highly non-linear nature of the optimization problem, which may allow for non-unique optimal solutions.Optimized maximum rates (such as k ao , k no , k den1 , and k den3 ) and exponential O 2 inhibition parameters for step-wise denitrification (K o 2 den2 and K o 2 den3 ) reveal more variability than half-saturation concentration coefficients (K terms), which often settle to the minimum or maximum allowed value (Table B1).
Pairwise correlations in Fig. 4 reveal several parameter pairs which exhibit strong relationships, reflecting the fact that in a significantly non-linear optimization, similar results can be obtained by trade-offs between different parameters and processes.Notably, the exponential O 2 inhibition constants for NO − 2 and N 2 O reduction (K o2 den2 and K o2 den3 , respectively) are strongly correlated with each other (R = 0.73) and with other parameters controlling the denitrification steps.These include positive correlations with the maximum rate parameters for NO − 3 and NO − 2 reduction (k den1 and k den2 , respectively) and negative correlations with the halfsaturation constants for NO − 2 and N 2 O reduction (K no2 den2 and K n2o den3 , respectively).These correlations suggests tight couplings between modeled denitrification steps, wherein high-/low maximum denitrification rates can be compensated by lower/higher half-saturation coefficients, respectively.
Considering the variability in the optimal parameter sets and the complexity of the cost function, which depends on observations for multiple variables at different depths, the resulting N cycle profiles show similar features across all optimal solutions (Fig. 5, top panels; see also Fig. C3 for macronutrient profiles).When compared to observations, the majority of parameter sets are able to skillfully model (1) the vertical distribution of O 2 , including the oxygen-deficient layer between roughly 100 to 400 m; (2) the subsurface maximum in NO − 2 ; (3) the rapid attenuation of NH + 4 with depth; and (4) the subsurface minimum in N * .N cycle transformation rates also show similar consistency in their vertical profiles, albeit with more notable discrepancies with observations, possibly reflecting the higher variability and more complex nature of these measurements.Lower rates than observed may also reflect the fact that incubation experiments provide potential rates rather than in situ rates.In general, the yield of N 2 O from NH + 4 oxidation (R n 2 o ao ) is O(100) times less than for NO − 2 (R no 2 ao ), following Eq.( A8) and (A9), consistent with observations (Ji et al., 2015a(Ji et al., , 2018a;;Santoro et al., 2021).The step-wise denitrification rates (R den1 , R den2 , and R den3 ) show remarkably similar vertical profiles, with higher NO − 3 reduction rates (R den1 ) and nearly identical magnitudes between R den2 and R den3 .Anammox (R ax ) shows a similar profile as denitrification, albeit with enhanced local maxima near the upper-and loweroxycline depths surrounding the OMZ core, consistent with observations (Kalvelage et al., 2013).Several robust features emerge from the optimized parameter solutions, suggesting underlying mechanisms that need to be captured for a faithful representation of the OMZ N cycle.In particular, the differences in the exponential O 2 inhibition parameters for denitrification, shown in Fig. 6 (left panel), reveal the existence of progressively lower O 2 tolerance for step-wise denitrification (K o2 den3 < K o2 den2 < K o2 den1 ) from all optimized parameter sets.As a result, denitrification can stop at either N 2 O or NO − 2 as O 2 increases above anoxic levels, leading to "incomplete" denitrification (Babbin et al., 2015).
Within the anoxic core of the OMZ (∼ 100 to 350 m depth), O 2 is low enough in all optimizations to allow each of the steps to proceed unimpeded (Fig. 5).The large differences between NO − 3 and NO − 2 reduction (R den1 − R den2 , middle panel of Fig. 6) allow accumulation of a characteristic subsurface peak in NO − 2 near the OMZ core.Conversely, N 2 O produced via NO − 2 reduction (R den2 ) is quickly consumed via N 2 O reduction (R den3 ), leading to a pronounced N 2 O deficit near the OMZ core.The progressive O 2 inhibition of the three steps of denitrification results in a decoupling between these reactions that is particularly evident in the oxycline layers above and below the OMZ, where N 2 O accumulation dominates, as N 2 O reduction (i.e., consumption) is more strongly inhibited by O 2 than NO − 2 reduction (i.e., N 2 O production, right panel of Fig. 6).Thus, the O 2 range defined by K o 2 den2 and K o 2 den3 can be thought of as a N 2 O production "window" that allows net N 2 O accumulation in the water column (Babbin et al., 2015).This O 2 -driven decoupling of anaerobic reactions is consistent with the observed sequential inhibition of N 2 O and N 2 production in incubation experiments (Dalsgaard et al., 2014), although we find O 2 inhibition thresholds that are somewhat higher than suggested by those experimental studies.Conversely, other studies have suggested much higher O 2 inhibition thresholds for anaerobic processes, on the order of several mmol m −3 (Kalvelage et al., 2011;Ji et al., 2018a).
The vertical profile of the step-wise denitrification rates (R den1 , R den2 , and R den3 ) shows remarkable agreement across solutions, with only a small subset of parameter sets that behave as outliers (Fig. 5).As a consequence, the fraction of POC remineralized by each heterotrophic reaction remains consistent across optimizations (Fig. 7, top panels).Near the base of the euphotic zone, at around 30 m depth, aerobic remineralization (R rem ) far exceeds denitrification, reflecting O 2 inhibition of the latter.However, as O 2 decreases to suboxic levels around 100 m depth, NO − 3 reduction becomes the dominant remineralization pathway (up to 60 % of total remineralization).As O 2 drops further within the OMZ core (∼ 100 to 350 m depth), NO − 2 and N 2 O re-duction rapidly take up the remaining fraction (∼ 25 % and 15 %, respectively), albeit with more variability than near the euphotic zone.Below the OMZ, as the water column reverts to oxic conditions, aerobic remineralization dominates, and by 500 m depth, all solutions show essentially no denitrification.
The processes responsible for fixed N loss (anammox, NO − 2 reduction, and N 2 O production from NH + 4 oxidation) are also consistent across optimizations (Fig. 7, bottom panels).Within oxygenated waters, N 2 O production from NH + 4 oxidation (R n2o ao ) is by far the dominant fixed N loss term, as all other sources are inhibited by O 2 .Anammox (R ax ) becomes the dominant term within the upper and lower oxycline due to increased availability of both NO − 2 (from denitrification and nitrification) and NH + 4 (from the decomposition of sinking POC), consistent with observations (Babbin et al., 2020).In the anoxic OMZ core, relief from O 2 inhibition allows NO − 2 reduction to outcompete anammox for NO − 2 and contributes up to 60 % of the total N loss, with anammox making up the remaining 40 % (also see Fig. 5).This is somewhat higher than expected from purely stoichiometric constraints (Koeve and Kähler, 2010;Bianchi et al., 2014), likely reflecting vertical transport of NO − 2 and NH + 4 , co-occurrence of aerobic and anaerobic processes, and the higher O 2 threshold for anammox inhibition in oxygenated waters.The resulting profile of total N loss thus reveals subsurface maxima predominantly driven by anammox, with denitrification leading total OMZ losses.

Selected solution for the eastern tropical South Pacific
Among tracers, N 2 O profiles show significant variability between optimizations.While all optimizations generate two peaks in N 2 O surrounding the oxygen-deficient core, only a subset is able to reproduce the observed magnitude of the secondary peak at the lower oxycline (roughly 500 m depth; see Fig. 5).This subset forms a "cluster"' of optimizations that share common features that facilitate the formation of a realistic deep N 2 O peak, including higher O 2 inhibition thresholds (between 1.0 and 2.0 mmol m −3 for NO − 2 reduction and between 0.5 and 1.0 mmol m −3 for N 2 O reduction) and a wider O 2 window where net N 2 O production is favored (between 0.5 and 1.0 mmol m −3 width).Additionally, while most optimizations are able to reproduce the OMZ peak in NO − 2 , significant variability in its magnitude exists.Given the central roles of N 2 O and NO − 2 in both nitrification and denitrification pathways (Fig. 1) and the importance of oceanic N 2 O emissions to the atmosphere, we assign high priority to optimizations that reproduce realistic features in the distribution of these tracers, in particular a higher magnitude for the secondary N 2 O maximum.To this end, we select a parameter set (hereafter Opt sel ), which results in N 2 O and NO − 2 profiles closer to observations (bold red curves in Fig. 5, with parameter values reported in Table B1).We use https://doi.org/10.5194/gmd-16-3581-2023 Geosci.Model Dev., 16, 3581-3609, 2023  this Opt sel parameter set for further analysis of the model sensitivity.
Compared to the other parameter sets, Opt sel is characterized by weaker maximum NH + 4 and NO − 2 oxidation rates (k ao and k no , respectively) and smaller half-saturation constants for reductant uptake (K nh 4 ao and K no 2 no , respectively) (Fig. C2).In surface oxygenated waters, this results in relatively higher NH + 4 and NO − 2 (Fig. 5).In contrast, maximum denitrifica-tion rates (k den1 , k den2 , and k den3 ) are close to the median values from all optimizations.Rates of NO − 2 and N 2 O reduction (R den2 and R den3 , respectively) are generally larger than other solutions, in particular near the lower oxycline (Fig. 5).This increases POC consumption within this depth range via denitrification compared to other solutions (Fig. 7).As a consequence, the residual between the NO − 3 and NO − 2 reduction (R den1 − R den2 ; see Fig. 6) leads to higher NO − 2 accumula- tion at these depths, providing the necessary NO − 2 substrate to fuel either NO − 2 reduction (i.e., N 2 O production) or anammox.Since the parameterization scheme in Opt sel also results in reduced NO − 2 oxidation (R no ) and anammox (R ax ) rates (see Fig. 5), likely because of higher anammox halfsaturation constants for substrate uptake (K nh 4 ax and K no 2 ax ), more NO − 2 is available for reduction by denitrification, leading to a surplus in production (R den2 ) relative to consumption (R den3 ) and high concentrations of N 2 O at the lower oxycline.

Sensitivities to model parameters
As shown in Sect.5.1 and Fig. 4, strong correlations exist between parameter pairs in the optimization ensemble.Since Opt sel demonstrates good comparisons with ETSP tracer and rate observations, we perform a series of sensitivity tests around parameters (P ) most responsible for controlling specific features (F ) of the tracer distributions.These include concentrations of NH + 4 and NO − 2 at 50 m depth, the peak NO − 2 concentration in the OMZ, the N 2 O concentrations at the primary and secondary N 2 O maxima, and the minimum in the OMZ NO − 3 deficit (i.e., N * ).Additionally, we evaluate which parameters govern total N loss, including the fractional contribution of anammox; the partitioning of POC consumption via NO − 3 , NO − 2 , and N 2 O reduction; and total N 2 O production and air-sea flux (here approximated by the vertical transport at the upper model boundary).To this end, we calculated the sensitivity coefficient (φ ij ) for each P and F pairing by evaluating the impact of varying each Opt sel P value by ±5 % of its range in Table B1 and recording the resulting relative change in the F : The results demonstrate high sensitivity to changes in the maximum rates for all reactions (Fig. 8).Specifically, higher maximum rates correlate negatively with the concentrations of their substrates and positively with the concentrations of their products.For example, increasing k den1 results in an increase in OMZ NO − 2 and a decrease in OMZ N * .Similarly, increasing k den2 decreases OMZ NO − to the atmosphere.These impacts are further modulated by the half-saturation and O 2 inhibition constants.
Figures 9 and 10 further summarize the sensitivities to the maximum denitrification rates and their inhibition by O 2 , detailing the resulting changes to O 2 , N 2 O, NO − 2 , and N * profiles.As expected, changes in maximum rates affect reaction substrates and products in opposite ways.For example, a positive perturbation of k den1 (top panels) stimulates NO − 3 reduction, causing an increase in OMZ NO − 2 and a decrease in N * as expected.Similarly, a positive perturbation of k den2 increases N 2 O and decreases NO − 2 nearly everywhere.However, these sensitivities also have specific depthdependent signatures.While changes in NO − 2 are more pronounced within the OMZ core, in particular the upper section, changes in N 2 O are stronger at the upper and lower oxyclines, i.e., within the N 2 O production window defined by K o 2 den2 and K o 2 den3 (see Sect. 5.1).Notably, by increasing k den1 (top panels in Fig. 9) or k den2 (middle panels) from Opt sel values, the vertical extent of oxygen-deficient waters is reduced as a result of increased POC consumption via denitrification (not shown).This enhances aerobic remineralization and nitrification below the OMZ, providing an enhanced source of NO − 3 that partly offsets the OMZ losses seen via k den1 enhancement.This may indicate a potential negative feedback: if denitrification is locally enhanced (i.e., via increased competition for POC by denitrifying heterotrophs), a resulting reduction in the vertical extent of the OMZ would inhibit further N loss.
Figures 8 and 10 highlight significant sensitivities to the O 2 inhibition constants, which control O 2 -dependent modulation of the maximum reaction rates.These effects are particularly evident at the boundaries of the OMZ.For example, an increase in K o 2 den2 allows for more NO − 2 reduction at higher O 2 , leading to a slight depletion in OMZ NO − 2 and, as a consequence, an increase in suboxic N 2 O concentrations (Fig. 10, middle panels), consistent with observations of these processes in the Peruvian oxygen-deficient zone (Frey et al., 2020).In a similar manner, an increase in K o 2 den3 leads to more N 2 O reduction, reducing the magnitude of both the primary and secondary N 2 O peaks, while leaving other OMZ tracers (NO − 2 , N * ) relatively unaffected.

Sensitivities to environmental variables
The main features of the OMZ simulated by the model are strongly dependent on environmental parameters such as upwelling and mixing; organic matter fluxes; and the model boundary conditions, including mixed-layer depth and O 2 concentrations.Critically, these parameters are likely to vary over time under the effects of natural climate variability (e.g., Deutsch et al., 2011) and anthropogenic climate change (Bopp et al., 2013).While each of these parameters control OMZ tracer profiles and N cycle reactions in complex ways, the main responses can be ascribed to changes in the position, thickness, and strength of the anoxic OMZ layer.Perturbations that replenish O 2 above the thresholds for anoxic processes -such as those predicted under climate warming scenarios (Busecke et al., 2022) -thus have cascading impacts on anaerobic N cycle intermediates, such as NO − 2 and N 2 O, and on the fixed N removal and NO − 3 deficit of the oxygen-deficient zone.
Figure 11 shows the sensitivity of the optimal solution Opt sel to the magnitudes of vertical upwelling (w up ) and turbulent diffusion (K v ).Increasing w up results in higher O 2 supply from below the OMZ, leading to increasing O 2 concentrations, and an upward shift and thinning of the anoxic layers.At high upwelling, the anoxic layer is effectively wiped out and is replaced by a suboxic layer.Similar results are obtained with higher K v values, with an increase in diffusive O 2 supply from both above and below the OMZ, resulting in a progressive shrinking of the anoxic layer.As this layer vanishes, anaerobic processes cease, drastically reducing the concentration of NO − 2 and the N deficit in the OMZ core.Notably, as the OMZ reaches the brink of anoxia, i.e., as the minimum O 2 concentration falls within the N 2 O production window, the upper and lower N 2 O maxima merge into a single N 2 O spike with particularly high N 2 O concentrations, reflecting the largest imbalance between production and consumption.

Discussion and conclusions
We developed a model of the N cycle in low O 2 waters and optimized it to reproduce observations from the ETSP OMZ.The model is able to simulate the distribution of multiple N cycle tracers, including NO − 2 and N 2 O, and their transforma-tion rates, capturing the underlying dynamics and environmental sensitivity of the underlying reactions (Fig. 5).In general, the model reproduces observed tracer concentration profiles more accurately than transformation rates.Mismatches with transformation rates may point to processes that need improvement in the model but also underscore limitations in rate measurements, which rely on shipboard incubation experiments that are usually more uncertain and limited than tracer measurements and may not perfectly reflect in situ conditions.However, by matching observed reaction rates to a reasonable degree, the model approximates the complex dynamics of the system in a way that allows it to reproduce tracer distributions.Co-located tracer and rate measurements for multiple processes are thus an effective way to constrain The optimization indicates that multiple parameter sets can produce equally good fits to tracer and rate profiles (Fig. 5).This is expected given the non-linear nature of the model and limitations in the observations.Even when rate measurements are used to constrain the model, as done here, an ensemble of equally good solutions is thus possible.This optimized ensemble shows that significant variability and trade-offs can exist between specific parameters (Fig. 4), suggesting that compensation between different processes can lead to similar profiles of tracers and transformation rates.Refinements to the criteria used to optimize the model, i.e., additional constraints in the definition of the cost function, could allow us to further narrow down plausible sets of parameters.For example, to evaluate the model sensitivity (Figs.8-10), we select a parameter set from our optimization ensemble that better captures the magnitude of the secondary N 2 O maximum, while reproducing other observed features equally well.While we adopt a relatively simple cost function definition, additional constraints such as this one could be explicitly built into its formulation and weighted more heavily to revise model parameters.
A better characterization of environmental sensitivities to substrate concentrations (e.g., half-saturation constant for substrate uptake) and O 2 sensitivities would also help parameter selection, for example, by narrowing down the prior and posterior range of values for these and other variables (e.g., maximum reaction rates).To this end, rate measurements under a range of O 2 and substrate concentrations are especially helpful.Similarly, simultaneous optimization of the model to reproduce observations across multiple regions of an OMZ characterized by different conditions, e.g., the core and the boundaries, or across different OMZ and oceanographic regimes would likely result in more robust optimizations.
Despite the variability in parameter values, analysis of the optimal ensemble reveals emerging features that appear robust across multiple optimizations and that compare well with observations.For example, the sensitivity of denitrification processes to O 2 shows systematic variations, with weaker O 2 inhibition for NO − 3 reduction and stronger for N 2 O reduction (Fig. 6).Accordingly, NO  with tracer incubation experiments (Dalsgaard et al., 2014).However, we note that the specific value of these O 2 sensitivities is far from well-established, with some experiments showing smaller thresholds than those found in our optimization (Dalsgaard et al., 2014) and others finding similar or larger thresholds (Ji et al., 2018a).In the model, the sequential sensitivity of denitrification steps to O 2 supports an O 2dependent window for N 2 O production, which allows accumulation of N 2 O at the margins of the OMZ core.This, and other systematic relationships between parameters and features of the solutions, as revealed by a sensitivity analysis (Figs. 8-10), sheds light on specific balances in the N cycle and can be exploited as a powerful tool to fine-tune the model, both in the one-dimensional setup used here and in more complex and resource-intensive three-dimensional implementations where a formal optimization would be unfeasible (McCoy et al., 2022).
Because the model is based on a mechanistic representation of N transformations, it is suitable for investigating the response of the N cycle to environmental variability and other perturbations .For example, the model could be used to investigate the effects of eddy variability near the boundaries of OMZs or the effects of OMZ expansion and change under global warming.With these goals in mind, the model is designed to be coupled to the biogeochemical component of the current generation of earth system models, enabling accurate simulation of NO − 2 and N 2 O dynamics, with minimal interference in the representation of the cycles of oxygen, nutrients, carbon, and organic matter.
Because the model reflects an evolving understanding of the N cycle, its assumptions should be re-evaluated as new N transformation processes and aspects of microbial dynamics are uncovered.The model is built around two major simplifications: the modularity of the N cycle and the representation of microbial metabolisms as bulk chemical reactions that avoid explicitly tracking diverse microbial populations.Both are approximate views of the N cycle.For example, recent evidence suggests that microorganisms with the ability to carry out intracellular reduction of NO − 3 to NO − 2 and NO − 2 to N 2 O may dominate production of N 2 O in oxygen-deficient waters (Ji et al., 2018a;Frey et al., 2020), although the sensitivity of this process to environmental factors is still being uncovered.
Our bulk approach assumes that metabolic reaction rates are proportional to substrates following a Michaelis-Menten dependency.However, in reality, reaction rates also depend on the abundance of microorganisms present in the water column.If microorganism biomass is assumed to be propor-tional to substrates, then a higher-order dependency of reaction rates may be more appropriate, as adopted by some biogeochemical models (e.g., Paulot et al., 2020).A different dependence on substrates, in turn, may affect the variability of reaction rates with depth and the model sensitivity to parameters such as maximum reaction rates.
Indeed, previous modeling studies have pointed out the value of explicitly resolving the biomass of microbial populations (Penn et al., 2016;Zakem et al., 2020).This, in turn, enables a more direct comparison of model results with molecular observations (Louca et al., 2016) and would favor the emergence of complex feedbacks between microbes and their substrates driven by resource competition and oceanic circulation (Penn et al., 2019).However, explicitly simulating microbial biomass requires a number of additional parameters that remain poorly constrained and adds computational burden that may not always improve the realism of biogeochemical simulations (Galbraith et al., 2015).Our model provides a valuable framework for continuing the exploration of these ideas in both idealized and realistic settings (McCoy et al., 2022).
Appendix A: NitrOMZ equations A1 Heterotrophic rate equations A2 Chemolithotrophic rate equations Production of N 2 O via the nitrification pathway in NitrOMZ (pathway 2b in Fig. 1) is modeled as a byproduct of R ao with enhanced yields at lower O 2 concentrations.The partitioning between N 2 O and NO − 2 production from R ao is calculated using the function proposed by Nevison et al. (2003), which was derived by fitting measured N 2 O and NO − 2 yields (Y n 2 o ao and Y no 2 ao , respectively) to oxygen concentrations (Goreau et al., 1980) and re-fit by multiple observations in the eastern tropical North and South Pacific OMZ (Ji et al., 2015a(Ji et al., , 2018a;;Santoro et al., 2021): Nitrification-derived NO − 2 and N 2 O production rates (R no 2 ao and R n 2 o ao , respectively; pathways 2a and 2b in Fig. 1) are therefore represented as

A4 Stoichiometry
The stoichiometry of heterotrophic redox reactions is based on an electron balance and follows the procedure outlined in Paulmier et al. (2009), under the assumption that the composition of organic matter (POC) follows the average oceanic ratios from Anderson and Sarmiento (1994)

Figure 2 .
Figure 2. Example of spinup of the model.(top) Temporal evolution of O 2 , N 2 O, NO − 2 , N * , and NH + 4 from initial ETSP boundary conditions at year 0 to the final model solution at year 700 using the selected parameter set (Opt sel ) discussed in Sect.5.2.Dashed black curves highlight the 1 and 10 mmol O 2 m −3 contours.(bottom) The same as (top) but for the heterotrophic rates of aerobic respiration (R rem ), NO − 3 reduction (R den1 ), NO − 2 reduction (R den2 ), and N 2 O reduction (R den3 ).The chemolithotrophic anammox rate (R ax ) is also shown in the far-right panel.

Figure 3 .
Figure 3. Flowchart of the CMA-ES optimization algorithm used to constrain uncertain model parameters.

Figure 4 .
Figure 4. Pairwise correlations between model parameters for model solutions optimized for the ETSP OMZ.See Table B1 for a list and a description of the model parameters.Correlation is shown as the Pearson correlation coefficient, with dots representing p values < 0.01.

Figure 5 .
Figure 5. Results from the optimized ensemble of model solutions.(top) Tracer (O 2 , N 2 O, NO − 2 , NH + 4 , N * , and N 2 ) profiles from all 382 optimized ETSP parameter sets.The bold red curves show the selected parameter set (Opt sel ) discussed in Sect.5.2.Observations used to define the optimization cost function are shown as circles in each panel.Macronutrient profiles (NO − 3 and PO 3− 4 ) are shown in Fig. C3.(bottom) The same as in (top) but for reaction rate profiles of N 2 O and NO − 2 production from NH + 4 oxidation (R n2o ao and R no2 ao , respectively); NO − 3 , NO − 2 , and N 2 O reduction (R den1 , R den2 , and R den3 , respectively); and anammox (R ax ).

Figure 6 .
Figure 6.Progressive O 2 inhibition of denitrification steps.(a) Histogram showing the distribution for all optimized solutions of the difference in the O 2 inhibition constant for NO − 3 and NO − 2 reduction (K o2 den1 and K o2 den2 , in dark gray) and NO − 2 and N 2 O reduction (K o2 den2 and K o2 den3 , in light gray).The small red markers denote the values from Opt sel .(b) Rate differences between NO − 3 and NO − 2 reduction (R den1 and R den2 ).Shading represents the 10/90 and 25/75 percentile at each vertical level from the 382 analyzed parameter sets.The bold red curves denote Opt sel results.(c) The same as (b) but for the difference in NO − 2 and N 2 O reduction rates (R den2 and R den3 ).

Figure 7 .
Figure 7. Contribution of different reactions to organic matter remineralization and fixed N loss.(top) Fraction of total POC remineralized by each heterotrophic rate (R rem , R den1 , R den2 , and R den3 ).Shading represents the 10/90 and 25/75 percentile at each vertical level from the 382 analyzed parameter sets.The bold red curves denote the selected parameter set (Opt sel ) discussed in Sect.5.2.(bottom) The same as (top) but for the fraction of total fixed N loss (via production of N 2 and N 2 O) from anammox (R ax ), NO − 2 reduction (R den2 ), and N 2 O production from NH + 4 oxidation (R n2o ao ).The total fixed N loss is also shown.Note the different vertical axes for the bottom panels.

Figure 8 .
Figure 8. Sensitivity coefficient (φ ij = P i F j • ∂F j ∂P i ) for Opt sel parameters (P i ) and features (F j ) of the model solution.Here, each parameter is varied by ±5 % of their respective CMA-ES-allowed ranges in Table B1 to evaluate the relative impact on each feature of the model solution.Concentrations of NH + 4 and NO − 2 at 50 m depth (−50 m) are used as proxies of near-surface values.

Figure 9 .
Figure 9. Model sensitivity to parameter values.Panels show O 2 , N 2 O, NO − 2 , and N * for the Opt sel parameter set after varying the maximum NO − 3 , NO − 2 , and N 2 O reduction rate parameters (k den1 , k den2 , and k den3 ) by ±5 % of their Opt sel value in TableB1.

Figure 10 .
Figure 10.Model sensitivity to parameter values.Panels show changes to N 2 O, NO − 2 , and N * for the Opt sel parameter set after varying the O 2 inhibition constants for NO − 3 , NO − 2 , and N 2 O reduction (K o2 den1 , K o2 den2 , and K o2 den3 ) by ±5 % of their Opt sel value in TableB1.Background gray shadings show O 2 concentrations, with horizontal lines highlighting O 2 = 1 mmol m −3 (dotted lines) and O 2 = 10 mmol m −3 (dashed lines).
upward O 2 supply by upwelling (Fig.12, bottom panel).Increasing bottom O 2 progressively decreases the thickness of the OMZ, shifting it upwards and eventually eroding the anoxic layer.Conversely, decreasing bottom O 2 leads to a downward expansion of the OMZ and an intensification of anoxic conditions and the resulting anaerobic reactions.

Figure 11 .
Figure 11.Model sensitivity to physical drivers.(top) Sensitivity of the Opt sel optimized solution to the constant vertical upwelling velocity (w up ).(bottom) Sensitivity to the vertical turbulent diffusion coefficient (K v ).The bold black curves indicate original Opt sel values, which are also indicated in their respective color bars.

Figure 12 .
Figure 12.Model sensitivity to biogeochemical drivers.The same as in Fig. 11 but for surface POC flux ( top poc ) and O 2 concentration at the lower model boundary.In the top panels, more negative values of top poc correspond to an increasing sinking POC flux.

Table 1 .
Summary of the main NitrOMZ parameters, with any reported values from the literature (NA -not available).
Bianchi et al.: A biogeochemical model of the ocean nitrogen cycle (Fig. 2).With newly available NO − 2 substrate and low-O 2 conditions, NO − 2 reduction (R den2 ) begins, resulting first in a subsurface spike in N 2 O.With further decrease in O 2 concentrations, N 2 O is reduced to N 2 , leading to a layer of low N 2 O concentrations within the OMZ that persists to the end of the simulation.Anammox (R ax ) is similarly relieved of O 2 inhibition as the O 2 minimum is established, reaching maximum values near the upper oxycline, reflecting a relatively high supply of both NO − 2 and NH + 4 .

Table B1 .
NitrOMZ nitrogen cycle parameters and CMA-ES optimization ranges.NA -not available.

Table B2 .
ETSP configuration for optimization routines.NA -not available.