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 biogeochemistry by limiting biological productivity in the surface ocean; influencing the cycles of other nutrients, carbon, and oxygen; and controlling oceanic emissions of nitrous oxide (N₂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${}_{\mathrm{2}}^{-}$) and N₂O. Here we present a new model of the oceanic N cycle (Nitrogen cycling in Oxygen Minimum Zones, 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 Pacific OMZ and evaluate the model's ability to reproduce observed profiles 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.

1 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₂); however, only fixed N, e.g., ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) and nitrate (${\mathrm{NO}}_{\mathrm{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).

Figure 1 Schematic of the main N cycle tracers and reactions represented by NitrOMZ. Tracers are shown in bold, ordered by the oxidation state of N, and consist of organic nitrogen (OrgN), ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$), nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), nitrous oxide (N₂O), and dinitrogen (N₂). N transformation reactions are shown as arrows connecting reactants and products. Green arrows represent nutrient uptake and nitrogen fixation rates (not explicitly represented in the model, which focuses on subsurface reactions). The black arrow corresponds to the (1, R_rem) release of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ by organic matter remineralization (i.e., ammonification). Blue arrows denote aerobic transformations: (2a, $R_{\mathrm{ao}}^{\mathrm{no}\mathrm{2}}$) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to ${\mathrm{NO}}_{\mathrm{2}}^{-}$, (2b, $R_{\mathrm{ao}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$) ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation to N₂O, and (3, R_no) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation to ${\mathrm{NO}}_{\mathrm{3}}^{-}$. Red arrows represent anaerobic processes: (4, R_den1) ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$, (5, R_den2) ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction to N₂O, (6, R_den3) N₂O reduction to N₂, and (7, R_ax) anaerobic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation (anammox) with ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to N₂.

The ocean's inventory of fixed N is dominated by ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and the main N cycle reactions consist of uptake and assimilatory reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (here used interchangeably with ammonia, NH₃) and the oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ back to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ following the decomposition of organic matter and nitrification (Fig. 1). Only when the concentration of dissolved oxygen (O₂) drops to suboxic or anoxic levels (typically below 5 mmol m⁻³) do additional metabolic pathways involving N become relevant, as observed in the ocean's oxygen minimum zones (OMZs) and low-O₂ sediments (Lam and Kuypers, 2011). These reactions include the three main steps of heterotrophic denitrification, i.e., the oxidation of organic carbon (OrgC) with ${\mathrm{NO}}_{\mathrm{3}}^{-}$, nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), and nitrous oxide (N₂O), and anammox, the chemolithotrophic oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ with ${\mathrm{NO}}_{\mathrm{2}}^{-}$. Both denitrification and anammox lead to the production of N₂ and thus remove fixed N from the ocean (Bianchi et al., 2012; DeVries et al., 2012, 2013). Ammonia oxidation is another source of N₂O – a powerful greenhouse gas and a leading agent of ozone destruction in the stratosphere. The number of N₂O molecules produced per NH₃ oxidized, i.e., the yield of this reaction, increases as O₂ declines (Goreau et al., 1980; Nevison et al., 2003), likely caused by a shift from N₂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₂O to the atmosphere (Naqvi et al., 2010; Yang et al., 2020), with the largest emissions observed right above shallow oxygen-deficient 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 reactions exchange nitrogen metabolites and substrates (Lam and Kuypers, 2011; Kuypers et al., 2018), ultimately controlling fixed nitrogen removal and N₂O production. While there is evidence that organic matter and O₂ 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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NH}}_{\mathrm{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₂O production (Aumont et al., 2015; Séférian et al., 2020; Stock et al., 2020; Long et al., 2021b).

The ultimate goals of these models are multifold and include improving predictability of oceanic N₂O emissions (Suntharalingam et al., 2012; Martinez-Rey et al., 2015; Battaglia and Joos, 2018; Buitenhuis et al., 2018; Ji et al., 2018a), providing a more realistic representation of the redox state of seawater (Louca et al., 2016), or resolving aspects of microbial dynamics underlying the oceanic N cycle (Penn et al., 2016; Zakem et al., 2018; Penn et al., 2019).

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₂O cycling often rely on parameterizations that link N₂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₂O sources and sinks by denitrification. These models also conflate anammox and denitrification into a single N₂ production term. Explicit cycling of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ under low O₂, with the observed co-occurrence of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ dissimilatory reactions; reduction to N₂O and N₂ by denitrification and anammox; and reoxidation to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ are missing (Lam and Kuypers, 2011; Kalvelage et al., 2013; Babbin et al., 2014, 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₂. 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.

2 Nitrogen cycle model formulation

2.1 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₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, N₂O, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, and organic nitrogen (OrgN) in either dissolved or particulate form. We only explicitly model ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (the dominant dissolved form) and do not distinguish it from NH₃. 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₂, 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) represented 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).

2.2 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₂ as the oxidant, and the three main steps of denitrification: dissimilatory ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (R_den1, pathway 4), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction to N₂O (R_den2, pathway 5), and N₂O reduction to N₂ (R_den3, pathway 6). Chemolithotrophic processes consist of aerobic oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ to both ${\mathrm{NO}}_{\mathrm{2}}^{-}$ ($R_{\mathrm{ao}}^{\mathrm{no}\mathrm{2}}$, pathway 2a) and N₂O ($R_{\mathrm{ao}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$, pathway 2b, via both hydroxylamine oxidation and nitrifier denitrification); aerobic oxidation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (R_no, pathway 3); and anammox, the anaerobic oxidation of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ with ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to produce N₂ gas (R_ax, pathway 7). Reactions are parameterized as functions of substrates (i.e., model tracer concentrations) and environmental parameters such as dissolved O₂ and organic matter. Tracers are expressed as concentrations, with units of mmol m⁻³.

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⁻³) 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 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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction to ${\mathrm{NH}}_{\mathrm{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₂O via ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation in NitrOMZ is represented as a single O₂-dependent function designed to model the transition in bacterial metabolisms from predominantly hydroxylamine oxidation to nitrifier denitrification at low O₂ (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₂O via a hybrid mechanism (Santoro et al., 2011; Löscher et al., 2012). Production of N₂O via AOA appears to be similarly enhanced at low O₂ (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 oxygen-deficient zones, where it may provide an additional source of ${\mathrm{NH}}_{\mathrm{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₂O from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ rather than ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Ji et al., 2018b; Frey et al., 2020). This suggests that denitrifying bacteria capable of direct production of N₂O from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (as ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction proceeds entirely within the cell) could be a major source of N₂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₂O from ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and its environmental sensitivity remain limited (Ji et al., 2018b; Frey et al., 2020).

Other work suggests the occurrence of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation in apparently O₂-deficient waters (Buchwald et al., 2015b; Babbin et al., 2020; Sun et al., 2021a). This may involve ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation coupled to iodate reduction or ${\mathrm{NO}}_{\mathrm{2}}^{-}$ disproportionation – two poorly characterized processes. It may also reflect the high affinity to O₂ of nitrite-oxidizing bacteria (Bristow et al., 2016) in regions where vanishing O₂ 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 organic 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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Callbeck et al., 2021) and anaerobic oxidation of methane with ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Thamdrup et al., 2019), both chemolithotrophic denitrification reactions.

2.3 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₂, based on an exponential limitation term (Dalsgaard et al., 2014). The resulting equation for a general heterotrophic reaction is

$$\begin{array}{}\text{(1)}& R_{\mathrm{H}}=k_{\mathrm{H}}\cdot {\displaystyle \frac{\left[X\right]}{\left[X\right]+K_{\mathrm{H}}^X}}\cdot e^{-\frac{{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{H}}^{\mathrm{o}\mathrm{2}}}}\cdot \mathrm{POC}.\end{array}$$

Here, H indicates the heterotrophic process considered (e.g., dissimilatory reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to ${\mathrm{NO}}_{\mathrm{2}}^{-}$), R_H the heterotrophic reaction rate (mmol C m⁻³ s⁻¹), k_H the specific first-order reaction rate (s⁻¹), [X] the concentration of the oxidant (i.e., O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, or N₂O), $K_{\mathrm{H}}^X$ the half-saturation constant for oxidant uptake (mmol m⁻³), $K_{\mathrm{H}}^{{\mathrm{o}}_{\mathrm{2}}}$ the scale for inhibition of the reaction by O₂ (mmol m⁻³), and POC the concentration of particulate organic matter in units of mmol C m⁻³. No O₂ inhibition is applied to aerobic respiration (i.e., $K_{\mathrm{H}}^{{\mathrm{o}}_{\mathrm{2}}}$ 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₂-dependent inhibition term limits the reactions when O₂ is present. The resulting equation for a general chemolithotrophic reaction is

$$\begin{array}{}\text{(2)}& R_{\mathrm{A}}=k_{\mathrm{A}}\cdot {\displaystyle \frac{\left[X\right]}{\left[X\right]+K_{\mathrm{A}}^X}}\cdot {\displaystyle \frac{\left[Y\right]}{\left[Y\right]+K_{\mathrm{A}}^Y}}\cdot e^{-\frac{{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{A}}^{{\mathrm{o}}_{\mathrm{2}}}}}\cdot \end{array}$$

Here, A indicates the chemolithotrophic process considered (e.g., anammox); R_A the reaction rate (mmol N m⁻³ s⁻¹); k_A the maximum reaction rate when the process is not limited (mmol N m⁻³ s⁻¹); [X] and [Y] the concentrations of the oxidant and reductant, respectively (e.g., ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ for anammox); $K_{\mathrm{A}}^X$ and $K_{\mathrm{A}}^Y$ the half-saturation constants for oxidant and reductant uptake, respectively (mmol m⁻³); and $K_{\mathrm{A}}^{{\mathrm{o}}_{\mathrm{2}}}$ the scale for inhibition of the reaction by O₂ (mmol m⁻³). For aerobic reactions, $K_{\mathrm{A}}^{{\mathrm{o}}_{\mathrm{2}}}$ is set to infinite, removing O₂ 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.

2.4 Model assumptions and parameterizations

In the model, we assume that heterotrophic reactions are first-order 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₂ 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₂ is based on the observation that they are limited at O₂ concentrations to a few mmol m⁻³ or smaller (Dalsgaard et al., 2014; Babbin et al., 2015; Frey et al., 2020). However, coexistence of anaerobic and aerobic reactions at O₂ concentrations of 10–20 mmol m⁻³ 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₂, and approximates empirical rates from incubation experiments reasonably well (Dalsgaard et al., 2014).

Parameter values for maximum reaction rates, half-saturation constants, and O₂ 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.

Babbin et al. (2015)Bristow et al. (2016)Peng et al. (2016)Ward (2008)Bristow et al. (2016)Sun et al. (2017)Sun et al. (2017)Bristow et al. (2016)Babbin et al. (2015)Martens-Habbena et al. (2009)Bristow et al. (2016)Peng et al. (2016)Sun et al. (2017)Bristow et al. (2016)Bristow et al. (2016)Dalsgaard et al. (2014)Ji et al. (2018a)Dalsgaard et al. (2014)Dalsgaard et al. (2014)Nevison et al. (2003)Santoro et al. (2021)Ji et al. (2018a)Nevison et al. (2003)Santoro et al. (2021)Ji et al. (2018a)

Table 1 Summary of the main NitrOMZ parameters, with any reported values from the literature (NA – not available).

3 Model implementation

3.1 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

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{\displaystyle \frac{\partial \left[C\right]}{\partial t}}& =-{\displaystyle \frac{\partial \left(w_{\mathrm{u}}\cdot \left[C\right]\right)}{\partial z}}+{\displaystyle \frac{\partial }{\partial z}}{K}_{\mathrm{v}}{\displaystyle \frac{\left(\partial \left[C\right]\right)}{\partial z}}\\ & +\sum _{i=\mathrm{1}}^{N_{\mathrm{H}}}\left(r_{C,\mathrm{H}}^i\cdot R_{\mathrm{H}}^i\right)+\sum _{i=\mathrm{1}}^{N_{\mathrm{A}}}\left(r_{C,\mathrm{A}}^i\cdot R_{\mathrm{A}}^i\right)+\mathrm{LT},\end{array}\end{array}$$

where w_u is the vertical upwelling velocity (m s⁻¹) and K_v is the vertical turbulent diffusion coefficient (m² s⁻¹, 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_{C,\mathrm{H}}^i$ and $r_{C,\mathrm{A}}^i$ 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):

$$\begin{array}{}\text{(4)}& \mathrm{LT}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_C}}\cdot \left(\left[C\right]-[C{]}^{\mathrm{far}}\right).\end{array}$$

For typical open ocean conditions, τ_C can be estimated as the minimum of an advective timescale $\frac{L}{U}$ and a diffusive timescale, $\frac{L^{\mathrm{2}}}{K_{\mathrm{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⁻¹, and K_H on the order of 1000 m² s⁻¹ results in a timescale τ_C=10⁸ 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).

3.2 Organic matter remineralization

In the one-dimensional model implementation, we represent organic matter (OrgC and OrgN) as a single particulate organic 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

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial \left(w_{\mathrm{s}}\cdot \mathrm{POC}\right)}{\partial z}}=-\sum _{i=\mathrm{1}}^{N_{\mathrm{H}}}{R}_{\mathrm{H}}^i=\sum _{i=\mathrm{1}}^{N_{\mathrm{H}}}\left(k_{\mathrm{H}}^{\mathrm{eff},i}\right)\cdot \mathrm{POC},\end{array}$$

where w_s is the depth-dependent sinking speed of POC in the water column, and $k_{\mathrm{H}}^{\mathrm{eff},i}$ (s⁻¹) is the effective rate constants for each heterotrophic process, i.e., the maximum rate constants multiplied by the respective substrate limitation and O₂ inhibition terms (Eq. 1).

Considering the flux of sinking POC, Φ_POC (mmol C m⁻² s⁻¹),

$$\begin{array}{}\text{(6)}& {\mathrm{\Phi }}_{\mathrm{POC}}=w_{\mathrm{s}}\cdot \mathrm{POC}.\end{array}$$

Equation (5) can be written as

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial {\mathrm{\Phi }}_{\mathrm{POC}}}{\partial z}}=-\sum _{i=\mathrm{1}}^{N_{\mathrm{H}}}{R}_{\mathrm{H}}^i=-\sum _{i=\mathrm{1}}^{N_{\mathrm{H}}}\left(k_{\mathrm{H}}^{\mathrm{eff},i}\right)\cdot \mathrm{POC}\end{array}$$

or, equivalently,

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\partial {\mathrm{\Phi }}_{\mathrm{POC}}}{\partial z}}=-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{H}}}\left(k_{\mathrm{H}}^{\mathrm{eff},i}\right)}{w_{\mathrm{s}}}}\cdot {\mathrm{\Phi }}_{\mathrm{POC}}.\end{array}$$

Equation (7) can be recast to relate the concentration of POC in the water column to the remineralization of the POC flux with depth:

$$\begin{array}{}\text{(9)}& \mathrm{POC}=-{\displaystyle \frac{\mathrm{1}}{{\sum }_{i=\mathrm{1}}^{N_{\mathrm{H}}}\left(k_{\mathrm{H}}^{\mathrm{eff},i}\right)}}\cdot {\displaystyle \frac{\partial {\mathrm{\Phi }}_{\mathrm{POC}}}{\partial z}}.\end{array}$$

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):

$$\begin{array}{}\text{(10)}& {\mathrm{\Phi }}_{\mathrm{POC}}={\mathrm{\Phi }}_{\mathrm{POC}}\left(z_{\mathrm{0}}\right)\cdot {\left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)}^{-b},\end{array}$$

where z₀ is the upper boundary of the model and b the power-law 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.

3.3 Numerical implementation of the one-dimensional model

For the purpose of testing and illustration, we implement NitrOMZ in a one-dimensional representation of the water column below the mixed layer, following previous work (Babbin et al., 2015). Model tracers are discretized on a one-dimensional vertical grid, with equal spacing Δz=10 m, where z is depth. Boundary conditions are set at the top (z₀) 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, ${\mathrm{\Phi }}_{\mathrm{POC}}\left(z_{\mathrm{0}}\right)=w_s\left(z_{\mathrm{0}}\right)\cdot \mathrm{POC}\left(z_{\mathrm{0}}\right)$, where POC(z₀) 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

$$\begin{array}{}\text{(11)}& w_{\mathrm{s}}={\displaystyle \frac{k_{\mathrm{Rem}}\cdot z}{b}}.\end{array}$$

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 ($\frac{\partial \left[C\right]}{\partial t}\approx \mathrm{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 $\frac{H}{w_{\mathrm{u}}}$ and $\frac{H^{\mathrm{2}}}{K_{\mathrm{v}}}$. For w_u on the order of 10 m y⁻¹, K_v on the order of 10⁻⁵ m² s⁻¹, and a vertical scale of 1000 m, the timescale to approach steady state is ${\mathit{\tau }}_{\mathrm{SS}}=\mathrm{3}\times {\mathrm{10}}^{\mathrm{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₂ leads the aerobic remineralization (R_rem) to dominate total POC consumption. As the simulation proceeds, an O₂ minimum develops in subsurface waters, reaching suboxic (<10 mmol O₂) concentrations around year 100. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction rates (R_den1) are relieved of O₂ inhibition and begin to take up a larger fraction of total POC remineralization, as revealed by the depletion of N^*, signaling ${\mathrm{NO}}_{\mathrm{3}}^{-}$ consumption in the water column. Reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ also leads to a subsurface peak in ${\mathrm{NO}}_{\mathrm{2}}^{-}$ within the O₂ minimum (Fig. 2). With newly available ${\mathrm{NO}}_{\mathrm{2}}^{-}$ substrate and low-O₂ conditions, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den2) begins, resulting first in a subsurface spike in N₂O. With further decrease in O₂ concentrations, N₂O is reduced to N₂, leading to a layer of low N₂O concentrations within the OMZ that persists to the end of the simulation. Anammox (R_ax) is similarly relieved of O₂ inhibition as the O₂ minimum is established, reaching maximum values near the upper oxycline, reflecting a relatively high supply of both ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$.

Figure 2 Example of spinup of the model. (top) Temporal evolution of O₂, N₂O, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, N^*, and ${\mathrm{NH}}_{\mathrm{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₂ m⁻³ contours. (bottom) The same as (top) but for the heterotrophic rates of aerobic respiration (R_rem), ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (R_den1), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den2), and N₂O reduction (R_den3). The chemolithotrophic anammox rate (R_ax) is also shown in the far-right panel.

4 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., 2017, 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.

4.1 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).

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

4.2 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, $\frac{K_{\mathrm{v}}}{w_{\mathrm{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₂ profile relative to ETSP observations, including the vertical position and thickness of oxygen-deficient waters (O₂<5 mmol m⁻³) (Fig. 5). These consist of the vertical diffusion and upwelling magnitude, the Martin curve coefficient (b), and the upper-ocean POC flux (${\mathrm{\Phi }}_{\mathrm{poc}}^{\mathrm{top}}$), 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₂ half-saturation constants for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation ($K_{\mathrm{ao}}^{{\mathrm{o}}_{\mathrm{2}}}$ and $K_{\mathrm{no}}^{{\mathrm{o}}_{\mathrm{2}}}$, 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 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₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$, which are central to the anaerobic N cycle. Because of possible influence from horizontal advection in observations, discrepancies exist between modeled and observed ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$. To compensate for this, we also assign lower weights to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ and higher weight to N^*. Second, a depth-dependent 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₂=1 mmol m⁻³) 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.

5 Results and sensitivity

5.1 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₂ inhibition parameters for step-wise denitrification ($K_{\mathrm{den}\mathrm{2}}^{{\mathrm{o}}_{\mathrm{2}}}$ and $K_{\mathrm{den}\mathrm{3}}^{{\mathrm{o}}_{\mathrm{2}}}$) 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₂ inhibition constants for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O reduction ($K_{\mathrm{den}\mathrm{2}}^{\mathrm{o}\mathrm{2}}$ and $K_{\mathrm{den}\mathrm{3}}^{\mathrm{o}\mathrm{2}}$, 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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (k_den1 and k_den2, respectively) and negative correlations with the half-saturation constants for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O reduction ($K_{\mathrm{den}\mathrm{2}}^{\mathrm{no}\mathrm{2}}$ and $K_{\mathrm{den}\mathrm{3}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$, 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.

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.

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₂, including the oxygen-deficient layer between roughly 100 to 400 m; (2) the subsurface maximum in ${\mathrm{NO}}_{\mathrm{2}}^{-}$; (3) the rapid attenuation of ${\mathrm{NH}}_{\mathrm{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₂O from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation ($R_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}$) is 𝒪(100) times less than for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ ($R_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}$), following Eq. (A8) and (A9), consistent with observations (Ji et al., 2015a, 2018a; Santoro et al., 2021). The step-wise denitrification rates (R_den1, R_den2, and R_den3) show remarkably similar vertical profiles, with higher ${\mathrm{NO}}_{\mathrm{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 lower-oxycline depths surrounding the OMZ core, consistent with observations (Kalvelage et al., 2013).

Figure 5 Results from the optimized ensemble of model solutions. (top) Tracer (O₂, N₂O, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NH}}_{\mathrm{4}}^{+}$, N^*, and N₂) 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 (${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) are shown in Fig. C3. (bottom) The same as in (top) but for reaction rate profiles of N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation ($R_{\mathrm{ao}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$ and $R_{\mathrm{ao}}^{\mathrm{no}\mathrm{2}}$, respectively); ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N₂O reduction (R_den1, R_den2, and R_den3, respectively); and anammox (R_ax).

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₂ inhibition parameters for denitrification, shown in Fig. 6 (left panel), reveal the existence of progressively lower O₂ tolerance for step-wise denitrification ($K_{\mathrm{den}\mathrm{3}}^{\mathrm{o}\mathrm{2}}<K_{\mathrm{den}\mathrm{2}}^{\mathrm{o}\mathrm{2}}<K_{\mathrm{den}\mathrm{1}}^{\mathrm{o}\mathrm{2}}$) from all optimized parameter sets. As a result, denitrification can stop at either N₂O or ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as O₂ 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₂ is low enough in all optimizations to allow each of the steps to proceed unimpeded (Fig. 5). The large differences between ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den1−R_den2, middle panel of Fig. 6) allow accumulation of a characteristic subsurface peak in ${\mathrm{NO}}_{\mathrm{2}}^{-}$ near the OMZ core. Conversely, N₂O produced via ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den2) is quickly consumed via N₂O reduction (R_den3), leading to a pronounced N₂O deficit near the OMZ core. The progressive O₂ 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₂O accumulation dominates, as N₂O reduction (i.e., consumption) is more strongly inhibited by O₂ than ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (i.e., N₂O production, right panel of Fig. 6). Thus, the O₂ range defined by $K_{\mathrm{den}\mathrm{2}}^{{\mathrm{o}}_{\mathrm{2}}}$ and $K_{\mathrm{den}\mathrm{3}}^{{\mathrm{o}}_{\mathrm{2}}}$ can be thought of as a N₂O production “window” that allows net N₂O accumulation in the water column (Babbin et al., 2015). This O₂-driven decoupling of anaerobic reactions is consistent with the observed sequential inhibition of N₂O and N₂ production in incubation experiments (Dalsgaard et al., 2014), although we find O₂ inhibition thresholds that are somewhat higher than suggested by those experimental studies. Conversely, other studies have suggested much higher O₂ inhibition thresholds for anaerobic processes, on the order of several mmol m⁻³ (Kalvelage et al., 2011; Ji et al., 2018a).

Figure 6 Progressive O₂ inhibition of denitrification steps. (a) Histogram showing the distribution for all optimized solutions of the difference in the O₂ inhibition constant for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction ($K_{\mathrm{den}\mathrm{1}}^{\mathrm{o}\mathrm{2}}$ and $K_{\mathrm{den}\mathrm{2}}^{\mathrm{o}\mathrm{2}}$, in dark gray) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O reduction ($K_{\mathrm{den}\mathrm{2}}^{\mathrm{o}\mathrm{2}}$ and $K_{\mathrm{den}\mathrm{3}}^{\mathrm{o}\mathrm{2}}$, in light gray). The small red markers denote the values from Opt_sel. (b) Rate differences between ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O reduction rates (R_den2 and R_den3).

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₂ inhibition of the latter. However, as O₂ decreases to suboxic levels around 100 m depth, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction becomes the dominant remineralization pathway (up to 60 % of total remineralization). As O₂ drops further within the OMZ core (∼100 to 350 m depth), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O reduction 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, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction, and N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation) are also consistent across optimizations (Fig. 7, bottom panels). Within oxygenated waters, N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation ($R_{\mathrm{ao}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$) is by far the dominant fixed N loss term, as all other sources are inhibited by O₂. Anammox (R_ax) becomes the dominant term within the upper and lower oxycline due to increased availability of both ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (from denitrification and nitrification) and ${\mathrm{NH}}_{\mathrm{4}}^{+}$ (from the decomposition of sinking POC), consistent with observations (Babbin et al., 2020). In the anoxic OMZ core, relief from O₂ inhibition allows ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction to outcompete anammox for ${\mathrm{NO}}_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$, co-occurrence of aerobic and anaerobic processes, and the higher O₂ 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.

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₂ and N₂O) from anammox (R_ax), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den2), and N₂O production from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation ($R_{\mathrm{ao}}^{\mathrm{n}\mathrm{2}\mathrm{o}}$). The total fixed N loss is also shown. Note the different vertical axes for the bottom panels.

5.2 Selected solution for the eastern tropical South Pacific

Among tracers, N₂O profiles show significant variability between optimizations. While all optimizations generate two peaks in N₂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₂O peak, including higher O₂ inhibition thresholds (between 1.0 and 2.0 mmol m⁻³ for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction and between 0.5 and 1.0 mmol m⁻³ for N₂O reduction) and a wider O₂ window where net N₂O production is favored (between 0.5 and 1.0 mmol m⁻³ width). Additionally, while most optimizations are able to reproduce the OMZ peak in ${\mathrm{NO}}_{\mathrm{2}}^{-}$, significant variability in its magnitude exists. Given the central roles of N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in both nitrification and denitrification pathways (Fig. 1) and the importance of oceanic N₂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₂O maximum. To this end, we select a parameter set (hereafter Opt_sel), which results in N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ profiles closer to observations (bold red curves in Fig. 5, with parameter values reported in Table B1). We use 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 ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation rates (k_ao and k_no, respectively) and smaller half-saturation constants for reductant uptake ($K_{\mathrm{ao}}^{{\mathrm{nh}}_{\mathrm{4}}}$ and $K_{\mathrm{no}}^{{\mathrm{no}}_{\mathrm{2}}}$, respectively) (Fig. C2). In surface oxygenated waters, this results in relatively higher ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Fig. 5). In contrast, maximum denitrification rates (k_den1, k_den2, and k_den3) are close to the median values from all optimizations. Rates of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (R_den1−R_den2; see Fig. 6) leads to higher ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulation at these depths, providing the necessary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ substrate to fuel either ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (i.e., N₂O production) or anammox. Since the parameterization scheme in Opt_sel also results in reduced ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (R_no) and anammox (R_ax) rates (see Fig. 5), likely because of higher anammox half-saturation constants for substrate uptake ($K_{\mathrm{ax}}^{{\mathrm{nh}}_{\mathrm{4}}}$ and $K_{\mathrm{ax}}^{{\mathrm{no}}_{\mathrm{2}}}$), more ${\mathrm{NO}}_{\mathrm{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₂O at the lower oxycline.

5.3 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 ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ at 50 m depth, the peak ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration in the OMZ, the N₂O concentrations at the primary and secondary N₂O maxima, and the minimum in the OMZ ${\mathrm{NO}}_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N₂O reduction; and total N₂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:

$$\begin{array}{}\text{(12)}& {\mathit{\varphi }}_{ij}={\displaystyle \frac{P_i}{F_j}}\cdot {\displaystyle \frac{\partial F_j}{\partial P_i}}.\end{array}$$

Figure 8 Sensitivity coefficient (${\mathit{\varphi }}_{ij}=\frac{P_i}{F_j}\cdot \frac{\partial F_j}{\partial 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 ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ at 50 m depth (−50 m) are used as proxies of near-surface values.

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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and a decrease in OMZ N^*. Similarly, increasing k_den2 decreases OMZ ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and increases N₂O concentrations in the upper and lower oxycline and its flux to the atmosphere. These impacts are further modulated by the half-saturation and O₂ inhibition constants.

Figures 9 and 10 further summarize the sensitivities to the maximum denitrification rates and their inhibition by O₂, detailing the resulting changes to O₂, N₂O, ${\mathrm{NO}}_{\mathrm{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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction, causing an increase in OMZ ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and a decrease in N^* as expected. Similarly, a positive perturbation of k_den2 increases N₂O and decreases ${\mathrm{NO}}_{\mathrm{2}}^{-}$ nearly everywhere. However, these sensitivities also have specific depth-dependent signatures. While changes in ${\mathrm{NO}}_{\mathrm{2}}^{-}$ are more pronounced within the OMZ core, in particular the upper section, changes in N₂O are stronger at the upper and lower oxyclines, i.e., within the N₂O production window defined by $K_{\mathrm{den}\mathrm{2}}^{{\mathrm{o}}_{\mathrm{2}}}$ and $K_{\mathrm{den}\mathrm{3}}^{{\mathrm{o}}_{\mathrm{2}}}$ (see Sect. 5.1).

Figure 9 Model sensitivity to parameter values. Panels show O₂, N₂O, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N^* for the Opt_sel parameter set after varying the maximum ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N₂O reduction rate parameters (k_den1, k_den2, and k_den3) by ±5 % of their Opt_sel value in Table B1.

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 ${\mathrm{NO}}_{\mathrm{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₂ inhibition constants, which control O₂-dependent modulation of the maximum reaction rates. These effects are particularly evident at the boundaries of the OMZ. For example, an increase in $K_{\mathrm{den}\mathrm{2}}^{{\mathrm{o}}_{\mathrm{2}}}$ allows for more ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction at higher O₂, leading to a slight depletion in OMZ ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and, as a consequence, an increase in suboxic N₂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_{\mathrm{den}\mathrm{3}}^{{\mathrm{o}}_{\mathrm{2}}}$ leads to more N₂O reduction, reducing the magnitude of both the primary and secondary N₂O peaks, while leaving other OMZ tracers (${\mathrm{NO}}_{\mathrm{2}}^{-}$, N^*) relatively unaffected.

Figure 10 Model sensitivity to parameter values. Panels show changes to N₂O, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N^* for the Opt_sel parameter set after varying the O₂ inhibition constants for ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N₂O reduction ($K_{\mathrm{den}\mathrm{1}}^{\mathrm{o}\mathrm{2}}$, $K_{\mathrm{den}\mathrm{2}}^{\mathrm{o}\mathrm{2}}$, and $K_{\mathrm{den}\mathrm{3}}^{\mathrm{o}\mathrm{2}}$) by ±5 % of their Opt_sel value in Table B1. Background gray shadings show O₂ concentrations, with horizontal lines highlighting O₂=1 mmol m⁻³ (dotted lines) and O₂=10 mmol m⁻³ (dashed lines).

5.4 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₂ 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₂ 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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O, and on the fixed N removal and ${\mathrm{NO}}_{\mathrm{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₂ supply from below the OMZ, leading to increasing O₂ 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₂ 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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and the N deficit in the OMZ core. Notably, as the OMZ reaches the brink of anoxia, i.e., as the minimum O₂ concentration falls within the N₂O production window, the upper and lower N₂O maxima merge into a single N₂O spike with particularly high N₂O concentrations, reflecting the largest imbalance between production and consumption.

Opposite changes are observed for a reduction in both w_up and K_v, which result in an expansion of the OMZ layer; increased ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and N₂O reduction; a larger OMZ ${\mathrm{NO}}_{\mathrm{2}}^{-}$ peak; and a broader separation of the upper and lower N₂O maxima. The interplay between the position of the oxygen-deficient layer, sinking particle fluxes, and transport processes further modulates the response of tracer profiles. For example, as anoxic waters expand upwards following a reduction in K_v, they intercept a higher concentration of sinking organic matter, which in turn fuels higher remineralization rates. Together with reduction in diffusive fluxes, this likely favors the strengthening of the upper N₂O maximum at low K_v observed in Fig. 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.

Because the supply of POC to the OMZ controls the overall magnitude of remineralization reactions, including O₂ consumption and denitrification, the model is particularly sensitive to the sinking POC flux at the upper model boundary (${\mathrm{\Phi }}_{\mathrm{poc}}^{\mathrm{top}}$, Table B3; Fig. 12, top panel). Increasing ${\mathrm{\Phi }}_{\mathrm{poc}}^{\mathrm{top}}$ causes a greater remineralization rate, which reduces available O₂, and drives a progressive thickening of the OMZ, with a series of cascading impacts on tracers similar to the ones discussed above. In contrast, decreasing ${\mathrm{\Phi }}_{\mathrm{poc}}^{\mathrm{top}}$ reduces the R_rem rates and increases O₂ to the point that anoxic conditions and their signature disappear.

Similar changes can also be driven by variations in the bottom-boundary O₂ concentration, which directly controls upward O₂ supply by upwelling (Fig. 12, bottom panel). Increasing bottom O₂ progressively decreases the thickness of the OMZ, shifting it upwards and eventually eroding the anoxic layer. Conversely, decreasing bottom O₂ leads to a downward expansion of the OMZ and an intensification of anoxic conditions and the resulting anaerobic reactions.

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

6 Discussion and conclusions

We developed a model of the N cycle in low O₂ 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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O, and their transformation 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 model representation of the N cycle in and around O₂-deficient environments.

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₂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₂ 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₂ 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₂ shows systematic variations, with weaker O₂ inhibition for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction and stronger for N₂O reduction (Fig. 6). Accordingly, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ tends to occur at higher O₂ concentrations than ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction to N₂O, which in turn occurs at higher O₂ concentrations than N₂O reduction to N₂. This result is consistent with tracer incubation experiments (Dalsgaard et al., 2014). However, we note that the specific value of these O₂ 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₂ supports an O₂-dependent window for N₂O production, which allows accumulation of N₂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 (Figs. 11–12). 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 ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂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 ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to N₂O may dominate production of N₂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 proportional 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).

Based on its modular design, the model can be naturally expanded to represent new processes that, while thought to be relevant in OMZ, are still uncertain. These include (1) additional known N cycle pathways and their sensitivity to environmental variability, such as DNRA (Lam et al., 2009), hybrid N₂O production from AOA (Stieglmeier et al., 2014), and direct ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to N₂O (Ji et al., 2018a; Frey et al., 2020); (2) alternative oxidation pathways, for example, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation with iodate or ${\mathrm{NO}}_{\mathrm{2}}^{-}$ disproportionation reactions (Babbin et al., 2020; Buchwald et al., 2015a; Sun et al., 2021a); (3) coupling of N tracers with the cycles of other elements, e.g., carbon, sulfur, and iron, such as chemolithotrophic denitrification coupled to hydrogen sulfide (H₂S) oxidation or anaerobic ${\mathrm{NO}}_{\mathrm{2}}^{-}$-based methane (CH₄) oxidation (Azhar et al., 2014; Scholz et al., 2016; Thamdrup et al., 2019; Callbeck et al., 2021); (4) explicit representation of chemolithotrophy and its effects on organic matter fixation (Swan et al., 2011); (5) explicit coupling to the inorganic carbon cycle by inclusion of CO₂ and alkalinity changes associated with N cycle reactions (Cinay et al., 2022); (6) the cycling of nitric oxide (NO) (Ward and Zafiriou, 1988; Lutterbeck et al., 2018); and (7) a more detailed representation of the microbial ecology underlying the N cycle (Louca et al., 2016; Zakem et al., 2018; Penn et al., 2019).

Appendix A: NitrOMZ equations

A1 Heterotrophic rate equations

$$\begin{array}{}\text{(A1)}& {\displaystyle }{R}_{\mathrm{rem}}=k_{\mathrm{rem}}\cdot {\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{rem}}^{{\mathrm{o}}_{\mathrm{2}}}+{\mathrm{O}}_{\mathrm{2}}}}\cdot \mathrm{POC}\text{(A2)}& {\displaystyle }{R}_{\mathrm{den}\mathrm{1}}=k_{\mathrm{den}\mathrm{1}}\cdot {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}^{-}}{K_{\mathrm{den}\mathrm{1}}^{{\mathrm{no}}_{\mathrm{3}}}+{\mathrm{NO}}_{\mathrm{3}}^{-}}}\cdot e^{\frac{-{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{den}\mathrm{1}}^{{\mathrm{o}}_{\mathrm{2}}}}}\cdot \mathrm{POC}\text{(A3)}& {\displaystyle }{R}_{\mathrm{den}\mathrm{2}}=k_{\mathrm{den}\mathrm{2}}\cdot {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{2}}^{-}}{K_{\mathrm{den}\mathrm{2}}^{{\mathrm{no}}_{\mathrm{2}}}+{\mathrm{NO}}_{\mathrm{2}}^{-}}}\cdot e^{\frac{-{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{den}\mathrm{2}}^{{\mathrm{o}}_{\mathrm{2}}}}}\cdot \mathrm{POC}\text{(A4)}& {\displaystyle }{R}_{\mathrm{den}\mathrm{3}}=k_{\mathrm{den}\mathrm{3}}\cdot {\displaystyle \frac{{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}{K_{\mathrm{den}\mathrm{3}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}+{\mathrm{N}}_{\mathrm{2}}\mathrm{O}}}\cdot e^{\frac{-{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{den}\mathrm{3}}^{{\mathrm{o}}_{\mathrm{2}}}}}\cdot \mathrm{POC}\end{array}$$

A2 Chemolithotrophic rate equations

$$\begin{array}{}\text{(A5)}& {\displaystyle }{R}_{\mathrm{ao}}=k_{\mathrm{ao}}\cdot {\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{ao}}^{{\mathrm{o}}_{\mathrm{2}}}+{\mathrm{O}}_{\mathrm{2}}}}\cdot {\displaystyle \frac{{\mathrm{NH}}_{\mathrm{4}}^{+}}{K_{\mathrm{ao}}^{{\mathrm{nh}}_{\mathrm{4}}}+{\mathrm{NH}}_{\mathrm{4}}^{+}}}\text{(A6)}& {\displaystyle }{R}_{\mathrm{no}}=k_{\mathrm{no}}\cdot {\displaystyle \frac{{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{no}}^{{\mathrm{o}}_{\mathrm{2}}}+{\mathrm{O}}_{\mathrm{2}}}}\cdot {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{2}}^{-}}{K_{\mathrm{no}}^{{\mathrm{no}}_{\mathrm{2}}}+{\mathrm{NO}}_{\mathrm{2}}^{-}}}\text{(A7)}& {\displaystyle }{R}_{\mathrm{ax}}=k_{\mathrm{ax}}\cdot {\displaystyle \frac{{\mathrm{NH}}_{\mathrm{4}}^{+}}{K_{\mathrm{ax}}^{{\mathrm{nh}}_{\mathrm{4}}}+{\mathrm{NH}}_{\mathrm{4}}^{+}}}\cdot {\displaystyle \frac{{\mathrm{NO}}_{\mathrm{2}}^{-}}{K_{\mathrm{ax}}^{{\mathrm{no}}_{\mathrm{2}}}+{\mathrm{NO}}_{\mathrm{2}}^{-}}}\cdot e^{\frac{-{\mathrm{O}}_{\mathrm{2}}}{K_{\mathrm{ax}}^{{\mathrm{o}}_{\mathrm{2}}}}}\end{array}$$

A3 Aerobic N₂O production

Production of N₂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₂ concentrations. The partitioning between N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production from R_ao is calculated using the function proposed by Nevison et al. (2003), which was derived by fitting measured N₂O and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ yields ($Y_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}$ and $Y_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}$, 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, 2018a; Santoro et al., 2021):

$$\begin{array}{}\text{(A8)}& {\displaystyle \frac{Y_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}}{Y_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}}}=\left({\displaystyle \frac{{\mathrm{Ji}}_a}{\left[{\mathrm{O}}_{\mathrm{2}}\right]}}+{\mathrm{Ji}}_b\right)\cdot \mathrm{0.01}.\end{array}$$

Nitrification-derived ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂O production rates ($R_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}$ and $R_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}$, respectively; pathways 2a and 2b in Fig. 1) are therefore represented as

$$\begin{array}{}\text{(A9)}& {\displaystyle }{R}_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}=R_{\mathrm{ao}}\cdot Y_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}},\text{(A10)}& {\displaystyle }{R}_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}=R_{\mathrm{ao}}\cdot Y_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}.\end{array}$$

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): C₁₀₆H₁₇₅O₄₂N₁₆P. This chemical composition can be arbitrarily adjusted in NitrOMZ. For example, studies in the eastern tropical South Pacific suggest a C : N ratio closer to 83:1 (Teng et al., 2014). Furthermore, organic matter degradation reactions may also differentially remineralize C, N, and P. For instance, denitrification may preferentially involve degradation of amino acids and thus impact the N : P ratio of remineralization differently from aerobic respiration (Van Mooy et al., 2002).

Based on the stoichiometry of Anderson and Sarmiento (1994), the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ : C and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ : C ratios during aerobic respiration are

$$\begin{array}{}\text{(A11)}& {\displaystyle }{Q}_{\mathrm{rem}}^{\mathrm{N}:\mathrm{C}}={\displaystyle \frac{\mathrm{16}}{\mathrm{106}}},\text{(A12)}& {\displaystyle }{Q}_{\mathrm{rem}}^{\mathrm{P}:\mathrm{C}}={\displaystyle \frac{\mathrm{1}}{\mathrm{106}}}.\end{array}$$

As a result of the POC composition, a total of 472 electrons are required to oxidize POC to CO₂. With four electrons required to reduce O₂ to H₂O, the oxygen : carbon remineralization ratio for aerobic remineralization to ${\mathrm{NH}}_{\mathrm{4}}^{+}$ is represented as

$$\begin{array}{}\text{(A13)}& Q_{\mathrm{rem}}^{\mathrm{O}:\mathrm{C}}={\displaystyle \frac{\mathrm{472}}{\mathrm{106}\cdot \mathrm{4}}}.\end{array}$$

This yields a respiration quotient of $r_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{C}}$ of 1.11, which is within the range of direct chemical measurements of $r_{-{\mathrm{O}}_{\mathrm{2}}:\mathrm{C}}$ from Moreno et al. (2020, 2022). For nitrification, the oxygen : nitrogen ratios for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (R_ao and R_no, respectively) are based on the stoichiometry of the relevant redox reactions:

$$\begin{array}{}\text{(A14)}& {\displaystyle }{Q}_{\mathrm{ao}}^{\mathrm{O}:\mathrm{N}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}},\text{(A15)}& {\displaystyle }{Q}_{\mathrm{no}}^{\mathrm{O}:\mathrm{N}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

For denitrification, two electrons are required for each respective reduction step (${\mathrm{NO}}_{\mathrm{3}}^{-}$ to ${\mathrm{NO}}_{\mathrm{2}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to $\frac{\mathrm{1}}{\mathrm{2}}$N₂O, and N₂O to N₂); thus the corresponding ratios are

$$\begin{array}{}\text{(A16)}& Q_{\mathrm{den}}^{\mathrm{N}:\mathrm{C}}={\displaystyle \frac{\mathrm{472}}{\mathrm{106}\cdot \mathrm{2}}}.\end{array}$$

Finally, for anammox, ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ are combined in 1:1 ratios to produce N₂. The above ratios are then applied to the tracer equations in Appendix A5.

A5 Tracer source-minus-sink equations

$$\begin{array}{}\text{(A17)}& {\displaystyle \frac{\mathrm{d}\left[\mathrm{POC}\right]}{\mathrm{d}t}}=-(R_{\mathrm{rem}}+R_{\mathrm{den}\mathrm{1}}+R_{\mathrm{den}\mathrm{2}}+R_{\mathrm{den}\mathrm{3}})\text{(A18)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\mathrm{d}t}}& =(Q_{\mathrm{rem}}^{\mathrm{O}:\mathrm{C}}\cdot R_{\mathrm{rem}})-(Q_{\mathrm{ao}}^{\mathrm{O}:\mathrm{N}}\cdot R_{\mathrm{ao}})\\ & -(Q_{\mathrm{no}}^{\mathrm{O}:\mathrm{N}}\cdot R_{\mathrm{no}})\end{array}\text{(A19)}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]}{\mathrm{d}t}}=R_{\mathrm{no}}-(Q_{\mathrm{den}}^{\mathrm{N}:\mathrm{C}}\cdot R_{\mathrm{den}\mathrm{1}})\text{(A20)}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]}{\mathrm{d}t}}=Q_{\mathrm{rem}}^{\mathrm{P}:\mathrm{C}}\cdot (R_{\mathrm{rem}}+R_{\mathrm{den}\mathrm{1}}+R_{\mathrm{den}\mathrm{2}}+R_{\mathrm{den}\mathrm{3}})\text{(A21)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{NH}}_{\mathrm{4}}^{+}\right]}{\mathrm{d}t}}& =Q_{\mathrm{rem}}^{\mathrm{N}:\mathrm{C}}\cdot (R_{\mathrm{rem}}+R_{\mathrm{den}\mathrm{1}}+R_{\mathrm{den}\mathrm{2}}+R_{\mathrm{den}\mathrm{3}})\\ & -(R_{\mathrm{ao}}+R_{\mathrm{ax}})\end{array}\text{(A22)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]}{\mathrm{d}t}}& =R_{\mathrm{ao}}^{{\mathrm{no}}_{\mathrm{2}}}+Q_{\mathrm{den}}^{\mathrm{N}:\mathrm{C}}\cdot (R_{\mathrm{den}\mathrm{1}}-R_{\mathrm{den}\mathrm{2}})\\ & -(R_{\mathrm{no}}+R_{\mathrm{ax}})\end{array}\text{(A23)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\mathrm{d}\left[{\mathrm{N}}_{\mathrm{2}}\mathrm{O}\right]}{\mathrm{d}t}}& =\mathrm{0.5}\cdot (R_{\mathrm{ao}}^{{\mathrm{n}}_{\mathrm{2}}\mathrm{o}}+Q_{\mathrm{den}}^{\mathrm{N}:\mathrm{C}}\cdot R_{\mathrm{den}\mathrm{2}})\\ & -(Q_{\mathrm{N}:\mathrm{C}}^{\mathrm{den}}\cdot R_{\mathrm{den}\mathrm{3}})\end{array}\text{(A24)}& {\displaystyle \frac{\mathrm{d}\left[{\mathrm{N}}_{\mathrm{2}}\right]}{\mathrm{d}t}}=(Q_{\mathrm{den}}^{\mathrm{N}:\mathrm{C}}\cdot R_{\mathrm{den}\mathrm{3}})+R_{\mathrm{ax}}\end{array}$$

Appendix B: NitrOMZ parameters and configurations

Nevison et al. (2003)Nevison et al. (2003)

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

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

Table B3 ETSP boundary conditions.

Table B4 Optimized ETSP parameter sets.

Appendix C: supplemental figures

Figure C1 (a) Opt_sel vertical diffusion (K_v) and (b) POC profiles.

Figure C2 Parameter distributions from the 382 CMA-ES-optimized ETSP solutions. Red markers denote Opt_sel values.

Figure C3 The same as in Fig. 5 but for O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, and N*.