Modeling oceanic nitrate and nitrite concentrations and isotopes using a 3-D inverse N cycle model

Nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) is a key intermediate in the marine nitrogen (N) cycle and a substrate in nitrification, which produces nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), as well as water column N loss processes denitrification and anammox. In models of the marine N cycle, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is often not considered as a separate state variable, since ${\mathrm{NO}}_{\mathrm{3}}^{-}$ occurs in much higher concentrations in the ocean. In oxygen deficient zones (ODZs), however, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ represents a substantial fraction of the bioavailable N, and modeling its production and consumption is important to understand the N cycle processes occurring there, especially those where bioavailable N is lost from or retained within the water column. Improving N cycle models by including ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is important in order to better quantify N cycling rates in ODZs, particularly N loss rates. Here we present the expansion of a global 3-D inverse N cycle model to include ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as a reactive intermediate as well as the processes that produce and consume ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in marine ODZs. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulation in ODZs is accurately represented by the model involving ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation, and anammox. We model both ¹⁴N and ¹⁵N and use a compilation of oceanographic measurements of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations and isotopes to place a better constraint on the N cycle processes occurring. The model is optimized using a range of isotope effects for denitrification and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation, and we find that the larger (more negative) inverse isotope effects for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation, along with relatively high rates of ${\mathrm{NO}}_{\mathrm{2}}^{-}$, oxidation give a better simulation of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations and isotopes in marine ODZs.

1 Introduction

Nitrogen (N) is an important nutrient to consider when assessing biogeochemical cycling in the ocean. The N cycle is intrinsically tied to the carbon (C) cycle, whereby N can be the limiting nutrient for primary production and carbon dioxide uptake (Moore et al., 2004; Codispoti, 1989). Understanding the distribution and speciation of bioavailable N in the ocean allows us to make inferences about the effects on other nutrient cycles and potential roles that N may play in a regime of climate change (Gruber, 2008).

There are several chemical species in which N can be found in the ocean. The largest pool of bioavailable N is nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$), a dissolved inorganic species, which can be taken up by microbes for use in assimilatory or dissimilatory processes. Another dissolved inorganic species, nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$), accumulates in much lower concentrations but is a key intermediate in many N cycling processes. Models of the marine N cycle often include ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ together as a single dissolved inorganic N (DIN) pool, or exclude ${\mathrm{NO}}_{\mathrm{2}}^{-}$ entirely (DeVries et al., 2013; Deutsch et al., 2007; Brandes and Devol, 2002). However, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ does accumulate significantly in oxygen deficient zones (ODZs) in features known as secondary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maxima, and it is an intermediate or substrate in many important N cycle processes occurring there.

ODZs are hotspots for marine N loss (Codispoti et al., 2001; Deutsch et al., 2007), which is driven by processes that result in conversion of bioavailable DIN to dinitrogen gas (N₂). The two main water column N loss processes, denitrification and anammox, use ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as a substrate. Denitrification involves the stepwise reduction of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and then to gaseous nitric oxide (NO), nitrous oxide (N₂O), and N₂. Anammox consists of the anaerobic oxidation of ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) to N₂ using ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as the electron acceptor. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is also oxidized to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ during anammox, representing an alternative fate for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in ODZs. Indeed, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation appears to be prevalent in ODZs, with more ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation occurring than can be explained by anammox alone (Gaye et al., 2013; Peters et al., 2016, 2018b; Babbin et al., 2017; Buchwald et al., 2015; Casciotti et al., 2013; Martin and Casciotti, 2017). ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation results in the regeneration of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ that would otherwise be converted to N₂ and lost from the system. The close coupling between ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation back to ${\mathrm{NO}}_{\mathrm{3}}^{-}$, represents a control valve on the marine N budget (Penn et al., 2016; Bristow et al., 2016). Where ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation can outcompete ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction via denitrification and anammox, bioavailable N is retained. Water column N losses may occur primarily where ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation rates are limited by oxygen availability. Thus, understanding the ${\mathrm{NO}}_{\mathrm{2}}^{-}$ dynamics in ODZ waters is critical to assess the N loss occurring there.

The observed ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations alone do not allow us to fully characterize the N cycling processes occurring in a given region. Stable isotope measurements of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ provide additional insight and constraints on marine N cycle processes. There are two stable isotopes of N: ¹⁴N and ¹⁵N. The isotopic ratios for a given N species, usually expressed in delta notation as δ¹⁵N (‰) = $\left(\right({}^{\mathrm{15}}\mathrm{N}/{}^{\mathrm{14}}\mathrm{N}{)}_{\mathrm{sample}}/({}^{\mathrm{15}}\mathrm{N}/{}^{\mathrm{14}}\mathrm{N}{)}_{\mathrm{standard}}-\mathrm{1})\times \mathrm{1000}$, are an integrated measure of the processes that have produced and consumed that N species. Each process imparts a unique isotope effect (ε (‰) = ${(}^{\mathrm{14}}k{/}^{\mathrm{15}}k-\mathrm{1})\times \mathrm{1000}$, where ¹⁴k and ¹⁵k are the first-order rate constants for the ¹⁴N and ¹⁵N containing molecules, respectively) that impacts the isotopic composition of the substrate and the product (Mariotti et al., 1981). In particular, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ cycling processes have distinct isotope effects, where ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction occurs with normal isotopic fractionation (Bryan et al., 1983; Martin and Casciotti, 2016; Brunner et al., 2013) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation occurs with an unusual inverse kinetic isotope effect (Casciotti, 2009; Buchwald and Casciotti, 2010; Brunner et al., 2013). Thus, the isotopes of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ are sensitive to the relative importance of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction in ${\mathrm{NO}}_{\mathrm{2}}^{-}$ consumption (Casciotti, 2009; Casciotti et al., 2013).

Models of the marine N cycle have employed isotopes and isotope effects in conjunction with N concentrations to elucidate N cycle processes (Brandes and Devol, 2002; Sigman et al., 2009; Somes et al., 2010; DeVries et al., 2013; Casciotti et al., 2013; Buchwald et al., 2015; Peters et al., 2016). A model can either assume a set of processes and infer the underlying isotope effects, or assume isotope effects and infer a set of processes. The latter isotope models are highly dependent on the chosen isotope effects used for given processes. Although there are estimates of isotope effects for processes based on both environmental measurements and laboratory studies, there is not always agreement between them. For example, laboratory cultures of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidizers indicate an N isotope effect of ${}^{\mathrm{15}}\mathit{\epsilon }=-\mathrm{10}$ to −20 ‰ (Casciotti, 2009; Buchwald and Casciotti, 2010), while measured concentrations and isotopes of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in ODZs indicate that isotope effects closer to −30 ‰ are needed to explain the observations (Buchwald et al., 2015; Casciotti et al., 2013; Peters et al., 2016).

Here we present an expansion of an existing global ocean 3-D inverse isotope-resolving N cycling model (DeVries et al., 2013) to investigate the isotopic constraints on N cycling in ODZs and the impact of these regions on global ocean N isotope patterns. An important step was to include ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and its isotopes as tracers. The addition of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ allows us to include additional internal N cycling processes, as well as a more nuanced and realistic version of the processes occurring in ODZs. We used a database of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ observations in order to assess the performance of the model as well as optimize the model N cycle parameters for which we do not have good prior estimates. In the model we employ a variety of isotope effect estimates for three important ODZ processes – ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction, and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation – to discern what isotope effects result in the best fit to the observations.

2 Methods

2.1 Inverse N cycle model overview

The model used here is a steady-state inverse model that solves for the concentration and δ¹⁵N of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, particulate organic N (PON), and dissolved organic N (DON) using a set of linear equations. Because the model assumes that the system is in steady state, it is not able to capture time-dependent properties of the system such as seasonality and anthropogenic change. However, on interannual timescales the N cycle is thought to be approximately in balance (Gruber, 2004; Bianchi et al., 2012). The residence time of N in the ocean, which is thought to be on the order of 2000–3000 years (Gruber, 2008), is sufficiently long to preclude any detectable changes in the global N inventory to date on timescales commensurate with the global overturning circulation. An important advantage of the steady-state assumption for our linear model is that it is possible to find solutions by direct matrix inversion without the need for a spin-up period as required by forward models. The solution to the system provides ¹⁴N and ¹⁵N concentrations of the N species of interest at every grid point in the model system. Working with a linear system imposes some restrictions on how complicated the rate equations can be, but there are improvements in model performance and ease of use, allowing us to test hypotheses about the processes that govern the marine N cycle and budget, particularly those occurring in and around oceanic ODZs. We aimed to produce a realistic N cycle model that represented ODZ processes accurately while limiting the number of free parameters. The description below outlines the dependencies and simplifications employed in this version of the model.

The model's uncertain biological parameters were determined through an optimization process that minimizes the difference between the modeled and observed ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration and isotope data. Computational time limits the number of parameters that we were able to optimize. We therefore focused our investigation on parameters that are poorly constrained by literature values and to which the model solution is most sensitive. In order to determine the parameters for optimization, a sensitivity analysis was performed on each parameter, varying them individually by ±10 % and computing the change in the modeled ¹⁴N and ¹⁵N. Those parameters that resulted in modeled ¹⁴N and ¹⁵N variability of >5 % were chosen for optimization in the model. The sensitivity analysis and the optimal values of the parameters contribute to an improved understanding of the cycling of N in the ocean in general and in the ODZs in particular. The optimization process is discussed in further detail in Sect. 2.6.

Table 1 Non-optimized model parameters.

^* Value used in ETSP O₂ sensitivity test (Sect. 4.2).

The sensitivity analysis revealed that the modeled distribution of ¹⁵N was very sensitive to chosen isotope effects, those parameters that control the relative rates of ¹⁵N and ¹⁴N in chemical and biological processes. There are literature estimates for each of the isotope effects of interest in this work, although there is often a discrepancy between isotope effects estimated in laboratory studies and those expressed in oceanographic measurements (Kritee et al., 2012; Casciotti et al., 2013; Bourbonnais et al., 2015; Martin and Casciotti, 2017; Fuchsman et al., 2017; Marconi et al., 2017; Peters et al., 2018b). Rather than optimizing the isotope effect values, we have chosen to use multiple cases with different combinations of previously estimated isotope effects in order to assess which values best fit the observations.

In addition to the optimized parameters and isotope effects, there were some nonsensitive parameters that were fixed prior to the optimization and whose values were chosen using literature estimates (Table 1). Some N cycle processes are also dependent on prescribed input fields that are not explicitly modeled, such as temperature, phosphate, oxygen, and net primary production. These external input fields will be discussed in detail in the relevant sections for each N cycle process.

2.2 Model grid and transport

The model uses a uniform 2^∘ × 2^∘ grid with 24 depth levels. The thickness of each model layer increases with depth, from 36 m at the top of the water column to 633 m near the bottom. Bottom topography was determined using 2 min gridded bathymetry (ETOPO2v2) that was then interpolated to the model grid. Our linear N cycle model relies on the transport of dissolved N species (${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and DON) in the ocean. For this we use the annual averaged circulation as captured by a tracer transport operator that governs the rate of transport of dissolved species (${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and DON) between boxes. The original version of the tracer data-assimilation procedure used to generate the transport operator for dissolved species (T_f) is described by DeVries and Primeau (2011), and the higher resolution version used here is described by DeVries et al. (2013).

2.3 N cycle

In the N cycling portion of the model, we track four different N species (Fig. 1). There are two organic N (ON) pools: dissolved (DON) and particulate (PON). There are also two dissolved inorganic N (DIN) pools: ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$. We did not explicitly model ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) because it typically occurs in low concentrations throughout the ocean, and scarcity of data (especially δ¹⁵N data) would make model validation difficult. Although ${\mathrm{NH}}_{\mathrm{4}}^{+}$ has been observed to accumulate to micromolar concentrations in some ODZs (Bristow et al., 2016; Hu et al., 2016), this occurs largely in shallow, coastal shelf regions that are not resolved by the model.

Because we used the concentrations of both ¹⁴N and ¹⁵N of each N species to constrain the rate parameters, two sets of governing equations were employed: one that depends on ¹⁴N and another that depends on ¹⁵N. Generally, the rate for ¹⁵N processes was dependent on the rate of ¹⁴N processes and an isotopic fractionation factor ($\mathit{\alpha }{=}^{\mathrm{14}}k{/}^{\mathrm{15}}k$) that is specific to each process and substrate. By solving for steady-state solutions to both ¹⁴N and ¹⁵N concentrations, we were able to model global distributions of [${\mathrm{NO}}_{\mathrm{3}}^{-}$], [${\mathrm{NO}}_{\mathrm{2}}^{-}$], and their corresponding δ¹⁵N values.

2.3.1 N cycle parameterization

We first describe the ¹⁴N equations and the general format of the N cycle in the model. Each equation is then broken down into its component parts for further explanation of the biological processes and their parameterization. The ¹⁵N equations and isotope implementation will be discussed in a later section.

Figure 1 Diagram showing the nitrogen (N) cycle processes represented in the model. Two organic N pools are modeled: particulate organic N (PON) and dissolved organic N (DON). Two inorganic N pools are modeled: nitrate (${\mathrm{NO}}_{\mathrm{3}}^{-}$) and nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$). N source processes are nitrogen (N₂) fixation and atmospheric deposition. N sink processes are sedimentary denitrification, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (NIR), and anammox (AMX). Internal cycling processes that transform N from one species to another are solubilization, remineralization, assimilation, ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (NAR), ammonia oxidation (AMO), and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (NXR). Neither ammonia (NH₃) nor ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) are tracked in this model, since they are assumed to not accumulate. N₂ is also not explicitly accounted for in the model.

The governing equations for the ¹⁴N-containing DIN (${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$) and organic N (DON and PON) state variables can be written as follows:

$$\begin{array}{ll}{\displaystyle }\left({\displaystyle \frac{\partial }{\partial t}}+T_{\mathrm{f}}\right)\left[{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}^{-}\right]& {\displaystyle }=J_{\mathrm{14}}^{\mathrm{dep}}-J_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{3}}}-J_{\mathrm{14}}^{\mathrm{NAR}}\\ \text{(1)}& {\displaystyle }& {\displaystyle }+J_{\mathrm{14}}^{\mathrm{NXR}}+\mathrm{0.3}{J}_{\mathrm{14}}^{\mathrm{AMX}}-J_{\mathrm{14}}^{\mathrm{sed}},{\displaystyle }\left({\displaystyle \frac{\partial }{\partial t}}+T_{\mathrm{f}}\right)\left[{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{2}}^{-}\right]& {\displaystyle }=J_{\mathrm{14}}^{\mathrm{AMO}}-J_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{2}}}+J_{\mathrm{14}}^{\mathrm{NAR}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }-J_{\mathrm{14}}^{\mathrm{NXR}}-J_{\mathrm{14}}^{\mathrm{NIR}}-\mathrm{1.3}{J}_{\mathrm{14}}^{\mathrm{AMX}},{\displaystyle }\left({\displaystyle \frac{\partial }{\partial t}}+T_{\mathrm{f}}\right)\left[\mathrm{DO}{}^{\mathrm{14}}\mathrm{N}\right]& {\displaystyle }=\mathit{\sigma }(J_{\mathrm{14}}^{\mathrm{fix}}+J_{\mathrm{14}}^{\mathrm{assim},\mathrm{WOA}})\\ \text{(3)}& {\displaystyle }& {\displaystyle }+J_{\mathrm{14}}^{\mathrm{sol}}-J_{\mathrm{14}}^{\mathrm{remin}},\end{array}$$

$$\begin{array}{ll}{\displaystyle }\left({\displaystyle \frac{\partial }{\partial t}}+T_{\mathrm{p}}\right)\left[\mathrm{PO}{}^{\mathrm{14}}\mathrm{N}\right]& {\displaystyle }=(\mathrm{1}-\mathit{\sigma })(J_{\mathrm{14}}^{\mathrm{fix}}+J_{\mathrm{14}}^{\mathrm{assim},\mathrm{WOA}})\\ \text{(4)}& {\displaystyle }& {\displaystyle }-J_{\mathrm{14}}^{\mathrm{sol}}.\end{array}$$

The model is designed to represent a steady state, thus the $\frac{\partial }{\partial t}$ term is 0. The J terms represent the source and sink processes for each state variable, expressed in units of mmol m⁻³ yr⁻¹ and will be described in more detail below. Briefly, $J_{\mathrm{14}}^{\mathrm{dep}}$ is the spatially variable deposition of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ from the atmosphere to the sea surface. In the DIN model equations, $J_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{3}}}$ and $J_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{2}}}$ represent the assimilation of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$, respectively, by phytoplankton in the upper two box levels. This assimilated ${\mathrm{NO}}_{\mathrm{3}}^{-}$ produces DON and PON, with proportions set by a spatially variable term, σ. Assimilation in the DON and PON equations is represented by $J_{\mathrm{14}}^{\mathrm{assim},\mathrm{WOA}}$ and is dependent on 2013 World Ocean Atlas (WOA) [${\mathrm{NO}}_{\mathrm{3}}^{-}$] interpolated to the model grid. N₂ fixation ($J_{\mathrm{14}}^{\mathrm{fix}}$) is split between DON and PON with the same σ term. ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction ($J_{\mathrm{14}}^{\mathrm{NAR}}$), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction ($J_{\mathrm{14}}^{\mathrm{NIR}}$), ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation ($J_{\mathrm{14}}^{\mathrm{NXR}}$), and anammox ($J_{\mathrm{14}}^{\mathrm{AMX}}$) act on the ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ pools. $J_{\mathrm{14}}^{\mathrm{sed}}$ represents the removal of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ via benthic denitrification. $J_{\mathrm{14}}^{\mathrm{sol}}$ represents the dissolution of PON into DON. $J_{\mathrm{14}}^{\mathrm{remin}}$ represents the degradation of DON, which feeds into ammonia oxidation ($J_{\mathrm{14}}^{\mathrm{AMO}}$) and $J_{\mathrm{14}}^{\mathrm{AMX}}$ as described below.

Through the use of these J terms, the governing equations are all linear with respect to the state variables. However, in order to introduce dependence of rates on the concentrations of multiple state variables, for example allowing heterotrophic ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction to be dependent on organic N as well as ${\mathrm{NO}}_{\mathrm{3}}^{-}$, we run the organic N equations and the DIN equations seperately. When [DON] is found in the [DIN] governing equations, that [DON] value has already been determined for each grid box from the organic N model. When [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is found in the DON governing equations, it is drawn from 2013 World Ocean Atlas annual data interpolated to the model grid.

2.3.2 N source processes

Atmospheric deposition and N₂ fixation are the two largest sources of N to the ocean (Gruber and Galloway, 2008) and the only sources of new bioavailable N in the model. We do not consider the third largest source of N, riverine fluxes, in the model due to a lack of coastal resolution and the expectation that much of the river-derived N is denitrified in the shelf sediments (Nixon et al., 1996; Seitzinger and Giblin, 1996). Representing these processes may be possible in a future version of the model, but is beyond the scope of the current model, given its coarse resolution near the coasts.

Atmospheric deposition

N deposition is assumed to only occur in the top box of the model. We assume that most of the N deposited is as ${\mathrm{NO}}_{\mathrm{3}}^{-}$, and that the other species would be rapidly oxidized to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ in the oxic surface waters.

$$\begin{array}{}\text{(5)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{dep}}=r_{\mathrm{14}}^{\mathrm{dep}}{S}^{\mathrm{dep}}\end{array}$$

To calculate $J_{\mathrm{14}}^{\mathrm{dep}}$, the atmospheric deposition rate of ¹⁴N, we use modeled total inorganic N deposition for 1993, S^dep (Galloway et al., 2004; Dentener et al., 2006; data available online at https://daac.ornl.gov/CLIMATE/guides/global_N_deposition_maps.html, last access: November 2017), which was interpolated to our model grid. This term, S^dep, is then multiplied by a prescribed fractional abundance of ¹⁴N in the deposited N ($r_{\mathrm{14}}^{\mathrm{dep}}$), which is calculated from the isotopic composition of deposited N (δ¹⁵N_dep, −4 ‰; Eq. 6), to yield the deposition of ¹⁴N to the sea surface in each box ($J_{\mathrm{14}}^{\mathrm{dep}}$). To calculate $r_{\mathrm{14}}^{\mathrm{dep}}$ from δ¹⁵N_dep, we first calculate $r_{\mathrm{15}}^{\mathrm{dep}}$ using $r_{\mathrm{15}}^{\mathrm{air}}$, a standard with a value of 0.003676 (Eq. 6; Mariotti, 1983).

$$\begin{array}{}\text{(6)}& {\displaystyle }{r}_{\mathrm{15}}^{\mathrm{dep}}=\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{dep}}}{\mathrm{1000}}}+\mathrm{1}\right)\times r_{\mathrm{15}}^{\mathrm{air}}\end{array}$$

Then, using the approximation that ¹⁵N∕¹⁴N = ${}^{\mathrm{15}}\mathrm{N}/{(}^{\mathrm{15}}\mathrm{N}{+}^{\mathrm{14}}\mathrm{N})$, we calculate $r_{\mathrm{14}}^{\mathrm{dep}}$ as ($\mathrm{1}-r_{\mathrm{15}}^{\mathrm{dep}}$). The units of S^dep are given in mg N m⁻² yr⁻¹, which we convert to mmol ${\mathrm{NO}}_{\mathrm{3}}^{-}$ m⁻³ yr⁻¹ by dividing by the depth of the surface box. This source term of N to the model is spatially variable but independent of the modeled N terms.

N₂ fixation

N₂ fixation is the other source of new N to the model, and is assumed to only occur in the top box of the model. It is parameterized similarly to N₂ fixation in the model of DeVries et al. (2013), with partial inhibition by ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Holl and Montoya, 2005) and dependence on iron (Fe) and phosphate (${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$) availability (Monteiro et al., 2011).

$$\begin{array}{}\text{(7)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{fix}}=r_{\mathrm{14}}^{\mathrm{fix}}{F}_{\mathrm{0}}{e}^{-{\mathrm{NO}}_{\mathrm{3},\mathrm{obs}}/\mathit{\lambda }}{e}^{\frac{T_{\mathrm{obs}}-T_{\max }}{T_{\mathrm{0}}}}{\displaystyle \frac{\mathrm{Fe}}{\mathrm{Fe}+K_{\mathrm{Fe}}}}{\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{{\mathrm{PO}}_{\mathrm{4}}+K_{\mathrm{P}}}}\end{array}$$

F₀ is the maximum rate of N₂ fixation (1.5 mmol m⁻³ yr⁻¹; Table 1) and is calculated from the estimated areal rate of N₂ fixation in the western tropical Atlantic (Capone et al., 2005) divided by the depth of the top model box. NO_3,obs is the 2013 World Ocean Atlas annually averaged surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ interpolated to the model grid (Garcia et al., 2013b). The parameter λ is an inhibition constant for N₂ fixation in the presence of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Table 1).

The temperature (T) terms scale the rate of N₂ fixation based on the observed temperature (T_obs), maximum observed sea surface temperature (T_max), and the minimum preferred growth temperature for Trichodesmium (T₀; Capone et al., 2005). The temperature data were taken from 2013 World Ocean Atlas annually averaged temperature interpolated to the model grid (Locarnini et al., 2013). We recognize that this will likely provide a conservative estimate of N₂ fixation, given the growing recognition of N₂ fixation outside of the tropical and subtropical ocean by organisms other than Trichodesmium (Shiozaki et al., 2017; Harding et al., 2018; Landolfi et al., 2018).

Fe is the modeled deposition of soluble Fe interpolated to the model grid (mmol Fe m⁻² yr⁻¹; Chien et al., 2016) divided by the depth of the top model grid box to give units of mmol Fe m⁻³ yr⁻¹. Fe and ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ are assumed to limit N₂ fixation at low concentrations via Michaelis–Menten kinetics. K_Fe and K_P are their respective half-saturation constants. Additionally, there is a term that allows us to set the isotopic ratio of newly fixed N, $r_{\mathrm{14}}^{\mathrm{fix}}$, which is the fractional abundance of ¹⁴N in newly fixed N and is calculated as in Eq. (6) from δ¹⁵N_fix (−1 ‰; Table 1). All of the N₂ fixation parameters are fixed rather than optimized (Table 1). Due to the use of non-optimized parameters and an input ${\mathrm{NO}}_{\mathrm{3}}^{-}$ field rather than modeled ${\mathrm{NO}}_{\mathrm{3}}^{-}$, N₂ fixation serves as an independent check that our modeled N cycle produces reasonable N concentrations and overall N loss rates. However, N₂ fixation is not explicitly modeled here and is instead taken as a fixed, though spatially variable, input field (Fig. S1 in the Supplement). The global rate of N₂ fixation produced by this parameterization is 131 Tg N yr⁻¹, which is in line with several current estimates (Table S1 in the Supplement).

In the model, N₂ fixation and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ assimilation (Sect. 2.3.3) are assumed to be the two processes that create exportable organic N. A fraction, σ, of this organic N is partitioned into DON rather than PON (Eqs. 3–4). In order to create spatial variability in this constant, we assumed that 1−σ, the fraction of assimilated N partitioned to PON, is equal to the particle export (P_e) ratio. This P_e ratio is the ratio of particle export to primary production, and is equivalent to the fraction of organic N that is exported from the euphotic zone as particulate matter rather than recycled or solubilized into DON. The P_e ratio is calculated for each model grid square from the mixed layer temperature (T_ml) and net primary production (NPP) as described by Dunne et al. (2005):

$$\begin{array}{}\text{(8)}& {\displaystyle }{P}_e=\mathit{\varphi }{T}_{\mathrm{ml}}+\mathrm{0.582}\log \left(\mathrm{NPP}\right)+\mathrm{0.419}.\end{array}$$

The constant Φ has a value of −0.0101 ^∘C⁻¹ as determined by Dunne et al. (2005). Net primary production estimates (in units of mmol carbon m⁻² yr⁻¹) were taken from a satellite-derived productivity model (Westberry et al., 2008), annually averaged, and interpolated onto the model grid. T_ml is calculated from the 2013 World Ocean Atlas annual average (Locarnini et al., 2013), which has been interpolated to the model grid. The temperature of the top two model boxes were averaged to give T_ml. As temperature increases, the P_e ratio decreases and less PON is exported, resulting in more DON recycling in the surface with several possible explanatory mechanisms discussed in greater detail by Dunne et al. (2005). As net primary production increases, the P_e ratio increases and relatively more PON is exported; net primary production explains 74 % of the observed variance in particle export (Dunne et al., 2005).

2.3.3 Internal N cycling processes

Assimilation of nitrate and nitrite

Assimilation accounts for the uptake of DIN and its incorporation into organic matter in the shallowest two layers of the global model. Since assimilation affects both the organic and inorganic N pools, we must account for it in both sets of model runs. We will first address assimilation in the organic N model (Eqs. 9 and 10).

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{assim},\mathrm{WOA}}{=}^{\mathrm{14}}{k}_{\mathrm{assim}}[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{obs}}\text{(10)}& {\displaystyle }& {}^{\mathrm{14}}{k}_{\mathrm{assim}}={\displaystyle \frac{\mathrm{NPP}}{r_{\mathrm{C}:\mathrm{N}}[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{obs}}}}\end{array}$$

Since the organic N model is run first and the assimilation rates are dependent on DIN concentrations, assumptions must be made about the DIN field in order to account for assimilation prior to the DIN model runs. Here we used observed [${\mathrm{NO}}_{\mathrm{3}}^{-}$] from the 2013 World Ocean Atlas annual product interpolated to the model grid [${\mathrm{NO}}_{\mathrm{3}}^{-}$]_obs (Garcia et al., 2013b) to calculate the assimilation rates for DON and PON production ($J_{\mathrm{14}}^{\mathrm{assim},\mathrm{WOA}}$). For this assumption to be valid, our modeled surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] must be close to the observed values, which we will test in Sect. 3.1. The rate constant for assimilation, ¹⁴k_assim, varies spatially and is determined using observations of surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and satellite-derived net primary production (NPP; Westberry et al., 2008). The rate constant is converted to N units using the ratio of carbon (C) to N in organic matter (r_C:N), which we assume to be 106:16 (Redfield et al., 1963). The value of the rate constant is only nonzero in the top two boxes of the model, where we assume primary production to be occurring. The same rate constant is used in both the organic N and DIN assimilation equations. We also assume from the perspective of organic N that only ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is being assimilated, since ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is present at relatively low concentrations in the surface ocean, and it may be characterized as recycled production. Assimilated N is partitioned between PON and DON using the P_e ratio as previously described and shown in Eqs. (3) and (4).

The setup for assimilation in the DIN model (Eqs. 11 and 12) is similar, but can use modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] rather than the World Ocean Atlas values. In order to appropriately reflect surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations, both ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ are assimilated. ¹⁴k_assim is calculated as described above and is assumed to be the same for both ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$. We justify using only [${\mathrm{NO}}_{\mathrm{3}}^{-}$] to parameterize ¹⁴k_assim because ${\mathrm{NO}}_{\mathrm{3}}^{-}$ generally makes up the bulk of DIN available for assimilation at the surface, but this assumption will be discussed in more detail below.

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{3}}}{=}^{\mathrm{14}}{k}_{\mathrm{assim}}{[}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}^{-}]\text{(12)}& {\displaystyle }& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{assim},{\mathrm{NO}}_{\mathrm{2}}}{=}^{\mathrm{14}}{k}_{\mathrm{assim}}{[}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{2}}^{-}]\end{array}$$

Solubilization

Solubilization is the transformation of PON to DON, and is dependent only on [PON] and a solubilization rate constant (¹⁴k_sol), which is optimized (Table 2).

$$\begin{array}{}\text{(13)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{sol}}{=}^{\mathrm{14}}{k}_{\mathrm{sol}}\left[\mathrm{PO}{}^{\mathrm{14}}\mathrm{N}\right]\end{array}$$

The solubilization of PON, together with the particle transport operator (T_p), produces a particle flux attenuation curve similar to a Martin curve with exponent $b=-\mathrm{0.858}$ (Table 1; Martin et al., 1987). While in the real world, the length scale for particle flux attenuation is somewhat longer in ODZs compared to oxygenated portions of the water column, and also varies regionally (Berelson, 2002; Buesseler et al., 2008; Buesseler and Boyd, 2009), our model uses a spatially invariant ¹⁴k_sol. A spatially variable ¹⁴k_sol that accounts for lower apparent values in ODZs is a refinement that could be introduced in future model versions.

Table 2 Optimized model parameters.

Remineralization

Remineralization, or ammonification, is the release of DON into the DIN pool. This is parameterized using the concentration of DON and a remineralization rate constant (¹⁴k_remin), which is optimized (Table 2).

$$\begin{array}{}\text{(14)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{remin}}{=}^{\mathrm{14}}{k}_{\mathrm{remin}}\left[{\mathrm{DO}}^{\mathrm{14}}\mathrm{N}\right]\end{array}$$

The removal of this remineralized DON, since it does not accumulate as ${\mathrm{NH}}_{\mathrm{4}}^{+}$, is either through ammonia oxidation (AMO) or anammox (AMX), depending on [O₂] as described below and in Sect. 2.3.4. We use the same remineralization rate constant regardless of the utilized electron acceptor (e.g., O₂, ${\mathrm{NO}}_{\mathrm{3}}^{-}$). Since particle flux attenuation is observed to be somewhat weaker in ODZs compared with oxygenated water (Van Mooy et al., 2002), this may slightly overestimate the rates of heterotrophic remineralization occurring in ODZs.

Ammonia oxidation (AMO)

AMO uses ammonia (NH₃) as a substrate. Since we do not include NH₃ or ${\mathrm{NH}}_{\mathrm{4}}^{+}$ in the model system, we treat remineralized DON as the substrate for AMO. In order to maintain consistency between the organic N and DIN model runs, remineralized DON is routed either to AMO or AMX (lost from the system) based on the O₂ dependencies of AMO and AMX. Rather than using a strict O₂ cutoff for AMO, it is limited by O₂ using Michaelis–Menten kinetics. The half-saturation constant for O₂, $K_m^{\mathrm{AMO}}$ (Table 1), sets the O₂ concentration at which AMO reaches half of its maximal value.

$$\begin{array}{}\text{(15)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{AMO}}=(\mathrm{1}-{\mathit{\eta }}_{\mathrm{AMX}})J_{\mathrm{14}}^{\mathrm{remin}}+{\mathit{\eta }}_{\mathrm{AMX}}{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_m^{\mathrm{AMO}}}}{J}_{\mathrm{14}}^{\mathrm{remin}}\end{array}$$

Recent studies have shown that AMO and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (NXR), both O₂-requiring processes, have very low O₂ half saturation constants and can occur down to nanomolar levels of [O₂] (Peng et al., 2016; Bristow et al., 2016). In contrast, O₂-inhibited processes such as AMX are only allowed to occur at O₂ concentrations below a given threshold. The handling of O₂ thresholds for anaerobic processes is discussed in more detail below (Sect. 2.3.4), though we describe it briefly here due to the interplay between AMO and AMX in the model. Briefly, the O₂ dependence of AMX is represented by the parameter η_AMX, which has a value between 0 and 1 for a given grid box depending on the average number of months in a year its 2013 World Ocean Atlas [O₂] falls below the [O₂] threshold for anammox (${\mathrm{O}}_{\mathrm{2}}^{\mathrm{AMX}}$, Table 1). If, for example, the [O₂] in a given grid box is always above the threshold for AMX, η_AMX=0 and all of the remineralized DON (represented by $J_{\mathrm{14}}^{\mathrm{rem}}$) will be oxidized via AMO. If [O₂] is less than ${\mathrm{O}}_{\mathrm{2}}^{\mathrm{AMX}}$, η_AMX will be nonzero and a smaller fraction of the remineralized DON will be oxidized via AMO. The fraction ultimately oxidized by AMO is thus determined by the Michaelis–Menten parameterization of AMO, as well as the O₂ threshold for anammox.

Nitrite oxidation

The rates of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (NXR) are dependent on the availability of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as well as O₂. Similar to AMO, we parameterize O₂ dependence using Michaelis–Menten kinetics and a fixed half-saturation constant for O₂ ($K_m^{\mathrm{NXR}}$, Table 1). $K_m^{\mathrm{NXR}}$ was taken to be 0.8 µM O₂, based on kinetics experiments performed with natural populations of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidizing bacteria (Bristow et al., 2016). Finally, we employ an optimized rate constant (¹⁴k_NXR, Table 2) to fit the available data.

$$\begin{array}{}\text{(16)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{NXR}}{=}^{\mathrm{14}}{k}_{\mathrm{NXR}}\left[{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{2}}^{-}\right]{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_m^{\mathrm{NXR}}}}\end{array}$$

2.3.4 N sink processes

Nitrate and nitrite reduction

${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (NAR) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (NIR) are two processes within the stepwise reductive pathway of canonical denitrification. The end result of denitrification is the conversion of DIN to N₂ gas, rendering it bioavailable to only a restricted set of marine organisms. Although there are intermediate gaseous products between ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and N₂, we treat NIR as the rate-limiting step in the denitrification pathway, where DIN is removed from the system.

For both NAR and NIR, we introduce a dependency on two state variables, their respective N substrates, and organic matter availability. Where NAR and NIR occur heterotrophically, they consume organic matter in addition to their main N substrates or electron receptors. When NAR occurs chemoautotrophically, it would be dependent primarily on the presence of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and an electron donor, such as hydrogen sulfide (Lavik et al., 2009). Since we do not model the production of reduced sulfur species in our model, our estimates of denitrification would not explicitly include the effects of this process. However, chemolithotrophic denitrification could be tacitly accounted for in the optimization process, since the rate constants that control the rates of NAR and NIR are optimized in order to best fit the observations, and the isotope effect for chemolithotrophic denitrification is thought to be similar to that of heterotrophic denitrification (Frey et al., 2014).

In order to maintain levels of heterotrophic NAR and NIR that are dependent on both the available ${\mathrm{NO}}_{\mathrm{3}}^{-}$ or ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and the available organic matter in a linear model, it was necessary to run organic N and DIN equations separately, since it is not possible to include dependencies on two state variables (e.g., DON and ${\mathrm{NO}}_{\mathrm{3}}^{-}$) in the linear system. Both NAR and NIR are dependent on the remineralization rate ($J_{\mathrm{14}}^{\mathrm{remin}}$) that is calculated in the organic N model run. In model boxes where NAR and NIR are occurring, some of the remineralization is carried out with electron acceptors other than O₂. As mentioned above, we assume that $J_{\mathrm{14}}^{\mathrm{remin}}$ does not depend on the choice of electron acceptor.

$$\begin{array}{}\text{(17)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{NAR}}& {\displaystyle }={\mathit{\eta }}_{\mathrm{NAR}}^{\mathrm{14}}{k}_{\mathrm{NAR}}\left[{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{3}}^{-}\right]J_{\mathrm{14}}^{\mathrm{remin}}\text{(18)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{NIR}}& {\displaystyle }={\mathit{\eta }}_{\mathrm{NIR}}^{\mathrm{14}}{k}_{\mathrm{NIR}}\left[{}^{\mathrm{14}}{\mathrm{NO}}_{\mathrm{2}}^{-}\right]J_{\mathrm{14}}^{\mathrm{remin}}\end{array}$$

The rate coefficients for NAR (¹⁴k_NAR) and NIR (¹⁴k_NIR) are optimized rather than fixed (Table 2). Further, the dependence of $J_{\mathrm{14}}^{\mathrm{NAR}}$ and $J_{\mathrm{14}}^{\mathrm{NIR}}$ on $J_{\mathrm{14}}^{\mathrm{remin}}$ means that k_NAR and k_NIR are not first-order rate constants and have different units than k_PON, k_DON, and k_NXR (Table 2).

The inhibition of NAR and NIR by O₂, like AMX, is parameterized by a parameter η, which inhibits these processes when [O₂] is above their maximum threshold. Originally, we treated this term as a binary operator that would be set to 0 if the empirically corrected 2013 World Ocean Atlas annually averaged [O₂] was above the threshold for the process and 1 if [O₂] was below the threshold. On further refinement, we wanted to account for the possibility of seasonal shifts in [O₂] in ODZs. Thus, for each month, we assigned a value of 0 or 1 to each model grid box. These values were then averaged over the 12 months of the year to give a sliding value of η between 0 and 1 for each grid box. The O₂ thresholds used to calculate η_NAR and η_NIR were fixed (7 and 5 µM, respectively; Table 1). Since we do not explicitly model O₂, [O₂] was predetermined using the 2013 World Ocean Atlas monthly O₂ climatology (Garcia et al., 2013a) interpolated to the model grid. We also applied an empirical correction that improves the fit of World Ocean Atlas [O₂] data to observed suboxic measurements (Bianchi et al., 2012).

Anammox

Anammox (AMX) catalyzes the production of N₂ from ${\mathrm{NH}}_{\mathrm{4}}^{+}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$. Since we do not use ${\mathrm{NH}}_{\mathrm{4}}^{+}$ as a variable in our N cycling equations, we substituted remineralized DON ($J_{\mathrm{14}}^{\mathrm{remin}}$) as a proxy for ${\mathrm{NH}}_{\mathrm{4}}^{+}$ availability. As described above in Sect. 2.3.3, remineralized DON is routed through either AMO or AMX depending on [O₂] and the O₂ dependencies of AMO and AMX.

$$\begin{array}{}\text{(19)}& {\displaystyle }{J}_{\mathrm{14}}^{\mathrm{AMX}}={\mathit{\eta }}_{\mathrm{AMX}}\left(\mathrm{1}-{\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}\right]}{\left[{\mathrm{O}}_{\mathrm{2}}\right]+K_m^{\mathrm{AMO}}}}\right)J_{\mathrm{14}}^{\mathrm{remin}}\end{array}$$

The O₂ threshold used to calculate η_AMX from monthly O₂ climatology is fixed (10 µM; Table 1). In order to maintain mass balance on remineralized DON, we do not include dependence on [${\mathrm{NO}}_{\mathrm{2}}^{-}$] in Eq. (19), although $J_{\mathrm{14}}^{\mathrm{AMX}}$ removes ${\mathrm{NO}}_{\mathrm{2}}^{-}$ (Eq. 2). This parameterization inherently assumes that AMX is limited primarily by [${\mathrm{NH}}_{\mathrm{4}}^{+}$] supply and not [NO₂], which may not always be correct (Bristow et al., 2016). Anammox also produces 0.3 moles of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ via associated NXR for every 1 mole of N₂ gas produced (Strous et al., 1999). For this reason, anammox appears in the state equation for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Eq. 1).

Sedimentary denitrification

Sedimentary (or benthic) denitrification ($J_{\mathrm{14}}^{\mathrm{sed}}$) is an important loss term for N in the marine environment, and in order to encapsulate it within the model grid we assume that it is occurring within the bottom depth box for any particular model water column. The parameterization for sedimentary denitrification is based on a transfer function described by Bohlen et al. (2012). The original transfer function was dependent on bottom water [O₂], bottom water [${\mathrm{NO}}_{\mathrm{3}}^{-}$], and the rain rate of particulate organic carbon (RRPOC). Here, RRPOC was calculated via a Martin curve (Martin et al., 1987) using the P_e ratio, net primary production (NPP), depth (z), euphotic zone depth (z_eu), and a Martin curve exponent (b):

$$\begin{array}{}\text{(20)}& {\displaystyle }\mathrm{RRPOC}=\mathrm{NPP}\cdot P_{\mathrm{e}}\cdot {\left({\displaystyle \frac{z}{z_{\mathrm{eu}}}}\right)}^b.\end{array}$$

Net primary production is derived from the productivity modeling of Westberry et al. (2008) as described in Sect. 2.3.2. The P_e ratio is calculated as previously described in Sect. 2.3.2. The depth for any given model box is assumed to be the depth at the bottom of the box. The euphotic zone depth is the bottom depth of the second box (73 m), since all production is assumed to be occurring in the top two boxes. As described above, the Martin curve exponent, b, is a fixed value in our model ($b=-\mathrm{0.858}$; Table 1), though this may result in underestimation of the particulate matter reaching the seafloor below ODZs (Van Mooy et al., 2002).

The transfer function for sedimentary denitrification was originally described using a nonlinear dependence of the rate on $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$. In order for sedimentary denitrification to be properly implemented in our linear model, we broke the original nonlinear relationship into three roughly linear segments to create a piecewise relationship between $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ and sedimentary denitrification rate (Fig. S2). We obtained three linear relationships between $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ and sedimentary denitrification rate, each applicable across a given range of $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ values. Due to the nature of our linear model, we needed to express the interval cutoff points that define the transition between the piecewise relationship segments in terms of O₂ rather than $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$. Therefore, a linear relationship between O₂ and $\left[{\mathrm{O}}_{\mathrm{2}}\right]-\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]$ was determined using the 2013 World Ocean Atlas annually averaged data (Garcia et al., 2013a, b; Fig. S3). The cutoff points were determined to be 75 and 175 µM O₂. The linear relationships were then rearranged in order to estimate sedimentary denitrification rate as a function of RRPOC, [O₂], and [${\mathrm{NO}}_{\mathrm{3}}^{-}$]. These equations were then further broken down into a component that is dependent on [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and a component that is dependent on [O₂] (see Supplement).

An additional term is introduced that reduces the sedimentary denitrification rate by 27 % if the depth of the bottom model box is less than 1000 m. This term represents the potential for efflux of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ into the water column from shallow, organic rich shelf sediments (Bohlen et al., 2012). This decreases overall sedimentary denitrification by approximately 6 Tg N yr⁻¹. This transfer function also assumes that all of the ${\mathrm{NH}}_{\mathrm{4}}^{+}$ efflux is immediately oxidized to ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and does not alter its isotopic composition in bottom water. This is a conservative estimate of the effects of benthic N loss on water column ${\mathrm{NO}}_{\mathrm{3}}^{-}$ isotopes, as several studies suggest that benthic N processes may contribute to water column nitrate ¹⁵N-enrichment (Lehmann et al., 2007; Granger et al., 2011; Somes and Oschlies, 2015; Brown et al., 2015). However, our current model parameterization does not require enhanced fractionation during benthic N loss to fit deep ocean ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$. Additionally, our spatial resolution does not well represent regions where this effect might be significant on bottom water ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, such as the shallow shelves.

2.4 N isotope implementation

In our model, we are interested in using the isotopic composition of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ to constrain the rates of N cycling and loss from the global ocean. As DON and PON are ultimate substrates for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NO}}_{\mathrm{3}}^{-}$ production, it is essential to track the ¹⁵N in the organic N pools as well. The matrix setup for ¹⁵N is similar to that for the ¹⁴N species, but the rates were changed as follows:

$$\begin{array}{}\text{(21)}& {\displaystyle }{J}_{\mathrm{15}}^{\mathrm{process}}=\mathrm{1}/{\mathit{\alpha }}_{\mathrm{process}}{\displaystyle \frac{\left[{}^{\mathrm{15}}{\mathrm{N}}_{\mathrm{substrate}}\right]}{\left[{}^{\mathrm{14}}{\mathrm{N}}_{\mathrm{substrate}}\right]}}{J}_{\mathrm{14}}^{\mathrm{process}}.\end{array}$$

$J_{\mathrm{14}}^{\mathrm{process}}$ is the rate of each relevant ¹⁴N process as described above, and $J_{\mathrm{15}}^{\mathrm{process}}$ is the rate of each ¹⁵N process. The α_process is the fractionation factor for a given process, which is given by the ratio between the rate constants for ¹⁴N and ¹⁵N ($\mathit{\alpha }{=}^{\mathrm{14}}k{/}^{\mathrm{15}}k$). A fractionation factor greater than 1 indicates a normal isotope effect and a fractionation factor less than 1 indicates an inverse isotope effect. Several of these fractionation factors are well known, but others are more poorly constrained, especially when values are calculated from in situ concentration and isotope ratio measurements (Hu et al., 2016; Casciotti et al., 2013; Ryabenko et al., 2012). For this reason, we ran several model cases with different fractionation factors for NAR, NIR, and NXR during the optimization process (Sect. 2.6, Table 3). The other fractionation factors were fixed (Table 1). In order to produce the ¹⁵N concentrations of N species from our observations to constrain the model, we calculated ¹⁵N∕¹⁴N from measured δ¹⁵N and multiplied by the measured concentration of each modeled N species, assuming that $\left[{}^{\mathrm{14}}\mathrm{N}\right]\sim \left[{}^{\mathrm{14}}\mathrm{N}\right]+\left[{}^{\mathrm{15}}\mathrm{N}\right]$.

Table 3 Isotope effect cases.

This simple ¹⁵N implementation was used with fixed fractionation factors for remineralization (α_remin=1), solubilization (α_sol=1), assimilation (α_assim=1.004), sedimentary denitrification (α_sed=1), and AMO (α_AMO=1) (Table 1). Isotope effects for NAR (ε_NAR), NIR (ε_NIR), and NXR (ε_NXR) were varied in different combinations during model optimization (Table 3). Distinct isotopic parameterizations were also required for atmospheric deposition, N₂ fixation, and anammox, as described below.

Atmospheric deposition

For atmospheric deposition of N, we prescribe fixed δ¹⁵N value of −4 ‰ (Table 1), which can be related to the fractional abundance of ¹⁴N ($r_{\mathrm{14}}^{\mathrm{dep}}$), previously described in Sect. 2.3.2, as well as the fractional abundance of ¹⁵N ($r_{\mathrm{15}}^{\mathrm{dep}}$) in deposited N. We multiply $r_{\mathrm{15}}^{\mathrm{dep}}$ by S^dep, the estimated rate of total N deposition to obtain $J_{\mathrm{15}}^{\mathrm{dep}}$.

Nitrogen fixation

Similar to atmospheric deposition, newly fixed N has a δ¹⁵N value (−1 ‰; Table 1). In Sect. 2.3.2 we described $r_{\mathrm{14}}^{\mathrm{fix}}$, the fractional abundance of ¹⁴N in newly fixed N. Here we multiply the fractional abundance of ¹⁵N, $r_{\mathrm{15}}^{\mathrm{fix}}$, by the other terms in the N₂ fixation equation (Eq. 6) to obtain the rate of ¹⁵N fixation.

Anammox

Anammox is the most complicated process to parameterize isotopically because it has three different N isotope effects associated with it. There is an isotope effect on both substrates that are converted to N₂ (${\mathrm{NO}}_{\mathrm{2}}^{-}$ and ${\mathrm{NH}}_{\mathrm{4}}^{+}$), as well as for the associated ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation to ${\mathrm{NO}}_{\mathrm{3}}^{-}$. We assume that the fractionation factor for ammonium oxidation via AMX (${\mathit{\alpha }}_{\mathrm{AMX},{\mathrm{NH}}_{\mathrm{4}}}$) is 1, setting it to match the fractionation factor for AMO (α_AMO; Table 1), both with no expressed fractionation since ${\mathrm{NH}}_{\mathrm{4}}^{+}$ does not accumulate in the model. Since all remineralized DON must be routed either through AMO or AMX, this simplifies the mass balance and ensures that all remineralized ¹⁴N and ¹⁵N is accounted for. ${}^{\mathrm{15}}{\mathrm{NO}}_{\mathrm{2}}^{-}$ is removed with the isotope effects of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (α_AMX,NIR) and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (α_AMX,NXR) in the expected 1:0.3 proportion (Brunner et al., 2013).

2.5 Model inversion

Once our N cycle equations were set up as described above, we input them into MATLAB in block matrix form. The equations were of the general form Ax=b. All model ocean boxes (200 160 in total) are accounted for in the matrices. Matrix A (400 320 × 400 320) contained the rate constants and other parameters that are multiplied by the vector of state variables, x (400 320 × 1). Vector x contained the state variables (i.e., [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] or [DON] and [PON]) to be solved for by the linear solver. Vector b (400 320 × 1) contained the rates that were independent of the state variables, such as N₂ fixation and N deposition. Let us consider, as an example, the DIN model setup. The top left corner of matrix A would contain rate constants for processes that produce and consume ${\mathrm{NO}}_{\mathrm{3}}^{-}$ that are also dependent on [${\mathrm{NO}}_{\mathrm{3}}^{-}$]. The top right corner of matrix A would contain rate constants for processes that produce and consume ${\mathrm{NO}}_{\mathrm{3}}^{-}$ but are dependent on [${\mathrm{NO}}_{\mathrm{2}}^{-}$]. The bottom left corner of matrix A would contain rate constants for processes that produce and consume ${\mathrm{NO}}_{\mathrm{2}}^{-}$ but are dependent on [${\mathrm{NO}}_{\mathrm{3}}^{-}$]. The bottom right corner of matrix A would contain rate constants that produce and consume ${\mathrm{NO}}_{\mathrm{2}}^{-}$ and are also dependent on [${\mathrm{NO}}_{\mathrm{2}}^{-}$]. The top half of vector x would be [${\mathrm{NO}}_{\mathrm{3}}^{-}$] for each model box, and the bottom half of vector x would be [${\mathrm{NO}}_{\mathrm{2}}^{-}$] for each model box. The top half of vector b would be independent processes that produce or consume [${\mathrm{NO}}_{\mathrm{3}}^{-}$], and the bottom half of vector b would be independent processes that produce or consume [${\mathrm{NO}}_{\mathrm{2}}^{-}$].

In MATLAB, we used METIS ordering, which is part of the SuiteSparse (http://faculty.cse.tamu.edu/davis/suitesparse.html, last access: December 2017) to order our large, sparse matrix A. We then used the built-in function “umfpack” with METIS to factorize matrix A. The built-in matrix solver “mldivide” was then used with the factorized components of matrix A and matrix b to solve for x.

2.6 Parameter optimization

There are many parameters in the model that control the rates of the different N cycle processes (Tables 1–3). Some of these parameters are well constrained by literature values (Table 1). Others, such as the rate constants, were objects of our investigation and were optimized against available observations (Table 2). For our optimization, we compared model output using different parameter values to a database of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentrations and isotopes. The database was originally compiled by Rafter et al. (2019) and has been expanded to include some additional unpublished data (Table S2). All of the database observations were binned and interpolated to the model grid. If multiple observations occurred within the same model grid box, the values were averaged and a standard deviation was calculated. The database was divided randomly into a training set, used for optimization, and a test set, used to assess model performance. The same number of grid points with observations was used in the training and test sets.

The optimization procedure used the MATLAB function “fminunc” to obtain values for the nonfixed parameters that minimized a cost function (Eq. 22). In each iteration of the optimization, the model system was solved by running the ¹⁴N-organic N model, ¹⁵N-organic N model, ¹⁴N-DIN model, and ¹⁵N-DIN model. The modeled output [${\mathrm{NO}}_{\mathrm{3}}^{-}$], [${\mathrm{NO}}_{\mathrm{2}}^{-}$], ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ were compared to values from the database training set. Though DON and PON observations were not used to optimize the model, the open ocean and deep water ${\mathrm{NO}}_{\mathrm{3}}^{-}$ values were useful in constraining the parameters that control PON solubilization and DON remineralization. The entire model was run using a set of initial parameter values (Table 2) and the optimization scheme continued to alter those starting parameters until a minimum in the cost function was attained. We optimized the logarithm of the parameter values rather than the original parameters themselves so the unconstrained optimization returned positive values. The transformed starting parameters and subsequent modified parameter sets were then fed back into the model equations as e^x, where x denotes the log-transformed parameter. The cost function in the optimization procedure is as follows:

$$\begin{array}{ll}{\displaystyle }\mathrm{Cost}& {\displaystyle }={\displaystyle \frac{w_{{\mathrm{NO}}_{\mathrm{3}}}}{n_{{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{sd}}_{{\mathrm{NO}}_{\mathrm{3}}}}}\sum \left(\right[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{model}}-[{\mathrm{NO}}_{\mathrm{3}}^{-}{]}_{\mathrm{training}}{)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{w_{{\mathrm{NO}}_{\mathrm{2}}}}{n_{{\mathrm{NO}}_{\mathrm{2}}}{\mathrm{sd}}_{{\mathrm{NO}}_{\mathrm{2}}}}}\sum \left(\right[{\mathrm{NO}}_{\mathrm{2}}^{-}{]}_{\mathrm{model}}-[{\mathrm{NO}}_{\mathrm{2}}^{-}{]}_{\mathrm{training}}{)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }+{\displaystyle \frac{w_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{3}}}}{n_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{3}}}{\mathrm{sd}}_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{3}}}}}\sum ({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}},\mathrm{model}}-{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}},\mathrm{training}}{)}^{\mathrm{2}}\\ \text{(22)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{w_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{2}}}}{n_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{2}}}{\mathrm{sd}}_{\mathit{\delta }{\mathrm{NO}}_{\mathrm{2}}}}}\sum ({\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{model}}-{\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}},\mathrm{training}}{)}^{\mathrm{2}}.\end{array}$$

The w terms are weighting terms introduced to scale the contributions of the four observed parameters to equalize their contributions to the cost function. The n terms and standard deviation (sd) terms were used to normalize the contributions of each measurement type to the cost function. Each n term is equal to the number of each type of measurement in the training dataset (e.g., the number of [${\mathrm{NO}}_{\mathrm{3}}^{-}$] data points $=n_{{\mathrm{NO}}_{\mathrm{3}}}$). The sd term is equal to the standard deviation of all the measurements of a given type (e.g., the standard deviation of all the [${\mathrm{NO}}_{\mathrm{3}}^{-}$] data points within the training set).

In order to account for error in our model parameter estimates, we also iterated over several possible values for three of the most important isotope effects for processes in ODZs: ε_NAR, ε_NIR, and ε_NXR (Table 3). We chose to iterate over these parameters rather than optimize them since there is a large range of estimates for each of these parameters. We assigned different possible values for each of these parameters (Table 3), resulting in 12 possible combinations. The optimization protocol was performed for each of those combinations and unique optimized parameter sets were obtained. The parameter results were then averaged (final values, Table 2) and their spread is categorized as the error (error, Table 2).

Figure 2 Map showing the model-estimated accumulation of nitrite (${\mathrm{NO}}_{\mathrm{2}}^{-}$) at 200 m depth.

3 Results

3.1 Global model–data comparison

The simulations of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration and isotopic composition are the most unique features of this model in comparison to existing global models of the marine N cycle. As such, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulation in ODZs is a feature that should be well represented by the model in order to use it to test hypotheses about processes that control N cycling and loss in ODZs. Overall, we see ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulating at 200 m in the major ODZs of the Eastern Tropical North Pacific (ETNP), Eastern Tropical South Pacific (ETSP), and the Arabian Sea (AS) (Fig. 2), which is consistent with observations and expected based on the low O₂ conditions found there. However, accumulation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the model ETSP was lower than expected. The model also accumulated ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the Bay of Bengal, which is a low-O₂ region off the east coast of India that does not generally accumulate ${\mathrm{NO}}_{\mathrm{2}}^{-}$ or support water column denitrification, but is thought to be near the “tipping point” for allowing N loss to occur (Bristow et al., 2017). Possible reasons for the underestimation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the ETSP and overestimation in the Bay of Bengal will be discussed further in Sect. 4.2.

The model optimization described above yielded a set of isotope effects that best fit the global dataset of [${\mathrm{NO}}_{\mathrm{3}}^{-}$], [${\mathrm{NO}}_{\mathrm{2}}^{-}$], ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$. The best fit was achieved for isotope effects of 13 ‰ for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (ε_NAR), 0 ‰ for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction (ε_NIR), and −13 ‰ for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation (ε_NXR). Figure 3 shows the test set comparison for the global best-fit set of isotope effects overlaid with a 1:1 line, which the data would follow if there was perfect agreement between model results and observations. There is general agreement between model and observations, with most of the data clustering near the 1:1 lines. Agreement between the observations and the training data are similar (Fig. S4), indicating that we did not overfit the training data.

Figure 3 Modeled (a) [${\mathrm{NO}}_{\mathrm{3}}^{-}$], (b) [${\mathrm{NO}}_{\mathrm{2}}^{-}$], (c) ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, and (d) ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ are compared against the corresponding values from the database test set. Shown on each panel is a 1:1 line starting at the origin. Data in black have corresponding [O₂]<10 µM, and data in gray have [O₂]≥10 µM.

In the test set, there were some low [O₂] points where our model [${\mathrm{NO}}_{\mathrm{3}}^{-}$] exceeded observations (Fig. 3a, filled black circles); these are largely within the ETSP. In contrast, the AS tended to show slightly lower modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] than expected. The [${\mathrm{NO}}_{\mathrm{2}}^{-}$] accumulation (Fig. 3b) and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ signals (Fig. 3c) in the ETSP were also generally too low compared with observations. These signals are likely tied to insufficient ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction occurring in the model ETSP. Another consideration is that there may be a mismatch in resolution between the model and the time and space scales needed to resolve the high ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulations observed sporadically (Anderson et al., 1982; Codispoti et al., 1985, 1986).

Overall, the representation of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was fairly good (RMSE = 2.4 ‰), though there were a subset of points above ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}=\mathrm{10}$ ‰ where the modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ exceeded the observed ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, and others where modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ was lower than observations (Fig. 3c). Many of the points with overestimated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ were located within the AS ODZ, where there may be too much ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction occurring, leading to artificially elevated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values. As indicated above, the underestimated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ points largely fell within the ETSP where we believe the model is underestimating ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction. The representation of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ was also fairly good (RMSE = 8.6 ‰), though the modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ was generally not low enough (Fig. 3d), indicating an underestimated sink of “heavy” ${\mathrm{NO}}_{\mathrm{2}}^{-}$.

3.2 ODZ model–data comparison using station profiles

To further investigate the distribution of model N species within the three main ODZs, we selected representative offshore grid boxes within each ODZ that contained observations to directly compare with model results in station profiles. Overall, the modeled ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ concentration and isotope profiles in the AS and ETNP were consistent with the observations, with [${\mathrm{NO}}_{\mathrm{3}}^{-}$] slightly underestimated in the AS ODZ and overestimated in the ETSP (Fig. 4). As [O₂] goes to zero, the O₂-intolerant processes NAR, NIR, and AMX are released from inhibition. These processes result in a decrease in [${\mathrm{NO}}_{\mathrm{3}}^{-}$] (via NAR) which corresponds to an increase in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, since NAR has a normal isotope effect. ${\mathrm{NO}}_{\mathrm{2}}^{-}$ also starts to accumulate in the secondary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maximum as a result of NAR. The ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ is lower than ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ since light ${\mathrm{NO}}_{\mathrm{2}}^{-}$ is preferentially created via NAR, and this fractionation is further reinforced by the inverse isotope effect of NXR (Casciotti, 2009). These patterns are readily observed in the AS and ETNP, but were less apparent in the ETSP, where [${\mathrm{NO}}_{\mathrm{3}}^{-}$] depletion and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] accumulation in the model were lower than observed. This could be due in part to the time-independent nature of this steady-state inverse model, which does not capture the effects of upwelling events in the ETSP on N supply and cycling (Canfield, 2006; Chavez and Messié, 2009).

Figure 4 Depth profiles comparing model results with binned and averaged database observations from a model water column. Results are shown for offshore regions of the three main oxygen deficient zones (ODZs): the Arabian Sea (AS; a–e), the Eastern Tropical North Pacific (ETNP; f–j), and the Eastern Tropical South Pacific (ETSP; k–o). Average modeled nitrate concentration ([${\mathrm{NO}}_{\mathrm{3}}^{-}$]), nitrite concentration ([${\mathrm{NO}}_{\mathrm{2}}^{-}$]), and N^* are shown in black. Gray error lines around the black line show the 2σ spread from the average from the 12 different optimized model results using the different combinations of isotope effects for nitrate reduction (ε_NAR), nitrite reduction (ε_NIR), and nitrite oxidation (ε_NXR). Observed data are shown in yellow in all panels. Modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ are shown for three different combinations of isotope effect. The blue lines represent ε_NAR=13, ${\mathit{\epsilon }}_{\mathrm{NXR}}=-\mathrm{13}$, and ε_NIR=0, which are the best fit isotope effects globally and in the ETSP. The red lines represent ε_NAR=13, ${\mathit{\epsilon }}_{\mathrm{NXR}}=-\mathrm{32}$, and ε_NIR=0, which are the best fit isotope effects in the AS.

Figure 5 Section profiles of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations and isotopes over the GEOTRACES GP16 cruise track (panel a inset) in the South Pacific. Section distance runs from west to east in each panel. Comparison of (a) observed [${\mathrm{NO}}_{\mathrm{3}}^{-}$] to (b) modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is presented over the full depth range (0–6000 m). Comparison of (c) observed ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ to (d) modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ is presented over a shortened depth range (0–1000 m) to better assess surface and ODZ values. GEOTRACES data are from Peters et al. (2018a) and available from Biological and Chemical Oceanography Data Management Office (BCO-DMO).

In order to gauge the model results for N loss, we also calculated N^*, a measure of the availability of DIN relative to ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ compared to Redfield ratio stoichiometry (${\mathrm{N}}^{*}=\left[{\mathrm{NO}}_{\mathrm{3}}^{-}\right]+\left[{\mathrm{NO}}_{\mathrm{2}}^{-}\right]-\mathrm{16}\cdot \left[{\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}\right]$; Deutsch et al., 2001). Negative N^* values are associated with N loss due to AMX or NIR or release of ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ from anoxic sediments (Noffke et al., 2012), while positive N^* values are associated with input of new N through N₂ fixation (Gruber and Sarmiento, 1997). Although we did not model ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$, we used the modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] together with World Ocean Atlas ${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$ data interpolated to the model grid to calculate N^* resulting from the model. Both the AS and ETNP showed a decrease in model N^* in the ODZ, as expected for water column N loss. Below the ODZ, N^* increased again and returned to expected deep water values. Modeled N^* in the ETSP, however, did not follow the observed trend, consistent with an underestimate of N loss in the model ETSP.

Though the global best fit isotope effects for NAR, NIR, and NXR produced good agreement to the data in general, the isotope effects that best fit individual ODZ regions differed when the cost function was restricted to observations from a given ODZ. For the ETSP, the best fit isotope effects were the same as the previously stated global best fit. For the AS, the best fit isotope effects were ε_NAR=13 ‰, ε_NIR=0 ‰, and ${\mathit{\epsilon }}_{\mathrm{NXR}}=-\mathrm{32}$ ‰. For the ETNP, the best fit isotope effects were ε_NAR=13 ‰, ε_NIR=15 ‰, and ${\mathit{\epsilon }}_{\mathrm{NXR}}=-\mathrm{32}$ ‰, though the performance is only marginally better than with ε_NIR=0 ‰. The lower (more inverse) value for ε_NXR resulted in higher ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ and lower ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$, which better fit the ODZ ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ data compared to the global best fit ${\mathit{\epsilon }}_{\mathrm{NXR}}=-\mathrm{13}$ ‰. These results are consistent with earlier isotope modeling studies in the ETSP (Casciotti et al., 2013; Peters et al., 2016, 2018b) and in the AS (Martin and Casciotti, 2017). Although, in the AS, modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values were too high, likely in part due to overpredicted rates of NAR, which also resulted in lower modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] (Fig. 4).

3.3 Model–data comparison in GEOTRACES sections

We also investigated the agreement between global best fit model concentration and isotope distributions with data from two GEOTRACES cruise sections: GP16 in the South Pacific, and GA03 in the North Atlantic. For GP16, we see that [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is low in surface waters and increases to a mid-depth maximum between 1000 and 2000 m. The highest [${\mathrm{NO}}_{\mathrm{3}}^{-}$] are found at mid-depth in the eastern boundary of the section. The model reproduces the general patterns, matching observations fairly well in the surface waters, but diverges below 500 m (Fig. 5). Although the patterns are generally correct, insufficient ${\mathrm{NO}}_{\mathrm{3}}^{-}$ is accumulated in the deep waters of the model Pacific. This could be due to an underestimate of preformed ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (over estimate of assimilation in the Southern Ocean), or inadequate supply of organic matter to be remineralized at depth. In the Southern Ocean, model surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] are 5–10 µM lower than observations (Fig. S5), which could be enough to explain the lower-than-expected [${\mathrm{NO}}_{\mathrm{3}}^{-}$] in the deep Pacific, which is largely sourced from the Southern Ocean (Rafter et al., 2013; Sigman et al., 2009; Peters et al., 2018a, b).

In the GP16 section, we also see that there are elevated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values in the model surface waters and in the ETSP ODZ (Fig. 5d), as expected from observations (Fig. 5c). However, we can also see that the insufficient depletion of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and increase in ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ in the ETSP ODZ (Fig. 5b and d) extends beyond the single grid box highlighted earlier (Fig. 4). The less-than-expected increase of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ in the ETSP ODZ and the upper thermocline in the eastern part of the section is consistent with an underestimate of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction. In GP16 we were also able to compare modeled and observed [${\mathrm{NO}}_{\mathrm{2}}^{-}$] and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ (Fig. S6). Patterns of modeled [${\mathrm{NO}}_{\mathrm{2}}^{-}$] and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ showed accumulation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the ODZ, with an appropriate ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ value (Fig. S6). Although, generally lower modeled concentrations of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the ODZ also support an underestimate of NAR (Fig. S6).

Figure 6 Section profiles of ${\mathrm{NO}}_{\mathrm{3}}^{-}$ concentrations and isotopes over the GA03 cruise track (panel a inset) in the North Atlantic. In each panel, section distance runs from west to east for the first 6000 km, and then runs from south to north. Comparison of (a) observed [${\mathrm{NO}}_{\mathrm{3}}^{-}$] to (b) modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is presented over the full depth range (0–6000 m). Comparison of (c) observed ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ to (d) modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ is presented over a shortened depth range (0–1000 m) to better assess surface and the low ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ contribution from N₂ fixation. GEOTRACES data are from Marconi et al. (2015) and available from BCO-DMO.

Surface ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values were also not as high in the model as in the observations (Fig. 5), which could result from insufficient ${\mathrm{NO}}_{\mathrm{3}}^{-}$ assimilation or too low supplied ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ (Peters et al., 2018a). However, we do see a similar depth range for high surface ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ and a local ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ minimum between the surface and ODZ propagating westward in both the model and observations, indicating that the physical and biogeochemical processes affecting ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ are represented by the model. Additionally, the model shows slightly elevated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ in the thermocline depths (200–500 m) west of the ODZ, which is consistent with the observations (Fig. 5c), though not of the correct magnitude. This is partly related to the muted ODZ signal as mentioned above and its lessened impact on thermocline ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ across the basin. Peters et al. (2018a) and Rafter et al. (2013) also postulated that these elevated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values were in part driven by remineralization of organic matter with high δ¹⁵N. The δ¹⁵N of sinking PON in the model (6 ‰–10 ‰) was similar to those observed in the South Pacific (Raimbault et al., 2008), as well as those predicted from aforementioned N isotope studies (Rafter et al., 2013; Peters et al., 2018a). The model also shows slightly elevated ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ in the intermediate depths (500–1500 m), which is consistent with observations, again reflecting remineralization of PON with δ¹⁵N greater than mean ocean ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$. Overall, the patterns of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ for the model GP16 are correct but the magnitudes of isotopic variation are muted, largely due to the lack of N loss in the ODZ and modeled surface ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values that are lower than observations. The simplification of ${\mathrm{NH}}_{\mathrm{4}}^{+}$ dynamics in the model could also contribute to underestimation of ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values if there was a large flux of ¹⁵N-enriched ${\mathrm{NH}}_{\mathrm{4}}^{+}$ from sediments (Granger et al., 2011), or if ¹⁵N-depleted ${\mathrm{NH}}_{\mathrm{4}}^{+}$ was preferentially transferred to the N₂ pool via anammox. While the isotope effect on ${\mathrm{NH}}_{\mathrm{4}}^{+}$ during anammox (Brunner et al., 2013) is indeed higher than that applied here, we chose to balance this with a low isotope effect during aerobic ${\mathrm{NH}}_{\mathrm{4}}^{+}$ oxidation (Table 1).

In the North Atlantic along GEOTRACES section GA03, we see good agreement between the observed and modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] (Fig. 6). There is generally low surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] with a distinct area of high [${\mathrm{NO}}_{\mathrm{3}}^{-}$] propagating from near the African coast. Deep water (>2000 m) [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is lower than we see in the Pacific section, and the model matches well with the Atlantic observations. Again, there is not quite enough ${\mathrm{NO}}_{\mathrm{3}}^{-}$ present in Southern Ocean-sourced intermediate waters (500–1500 m; Marconi et al., 2015). Modeled ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values at first glance appear higher than observed values at the surface (Fig. 6). However, many of the surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] were below the operating limit for ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ analysis and were not determined. Focusing on areas where both measurements and model results are present yields excellent agreement. For example, we do see low ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ values in upper thermocline waters in both the model and observations, likely corresponding to low δ¹⁵N contributions from N₂ fixation that is remineralized at depth and accumulated in North Atlantic Central Water (Marconi et al., 2015; Knapp et al., 2008). The model input includes significant rates of N₂ fixation in the North Atlantic that are consistent with this observation (Fig. S1). However, rates of N deposition in the North Atlantic are also fairly high and can contribute to the low δ¹⁵N signal (Knapp et al., 2008). In our model, atmospheric N deposition contributed between 0 % and 50 % of N input along the cruise track.

4 Discussion

4.1 Assumption checks

As previously mentioned (Sect. 2.3), organic N and DIN were modeled separately in order to introduce dependence on both organic N and substrate availability for the heterotrophic processes NAR and NIR. These separate model runs required several assumptions to be made regarding the processes that impact both organic N and DIN, namely assimilation and remineralization.

The first assumption was that the rates of N assimilation are equal between the organic N and DIN model runs. The organic N model run uses World Ocean Atlas surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] to estimate the contribution of DIN assimilation to the production of organic N, whereas the DIN model uses modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] to estimate DIN removal via assimilation. Though these two methods used the same rate constants for assimilation, differences in DIN concentrations could cause some discrepancies between the overall rates. Analysis of the results revealed that slightly more overall DIN assimilation occurred in the DIN model run than organic N produced in the organic N model (Fig. S7). This could be due in part to assimilation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the top two boxes of the DIN model, since ${\mathrm{NO}}_{\mathrm{2}}^{-}$ assimilation is unaccounted for in the organic N model. This is largely an issue in the oligotrophic gyres, where surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is very low and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulates to low but nonzero concentrations (Fig. 2). Assimilation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accounts for a significant fraction of DIN assimilation in these regions, but the overall assimilation rates there are low and the resulting influence on the whole system is also low. Another source of discrepancy would be where modeled surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] is higher than the World Ocean Atlas surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] that is supplied to the organic N model, which would result in higher assimilation rates in the DIN model run. Indeed, points at which the DIN assimilation rates are higher than the organic N production rates do tend to have modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] that was higher than observed [${\mathrm{NO}}_{\mathrm{3}}^{-}$] (Fig. S7). Likewise, points with relatively lower DIN assimilation had modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$] less than observed [${\mathrm{NO}}_{\mathrm{3}}^{-}$]. However, the majority of DIN assimilation estimates were within 10 µM yr⁻¹ of the organic N production estimates, with an average offset of approximately 3.5 % compared to DIN assimilation. The total global assimilation rates were within 0.4 %, with some spatially variable differences due to offset between surface [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and modeled [${\mathrm{NO}}_{\mathrm{3}}^{-}$]. However, we find that the World Ocean Atlas surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ values are fairly well represented by our modeled surface ${\mathrm{NO}}_{\mathrm{3}}^{-}$ (Fig. S5). We conclude that though the assimilation rates are not identical in the organic N and DIN model runs, the discrepancy in modeled DIN assimilation is less than 0.1 %, and there is unlikely to be significant creation or loss of N as a result of the split model.

4.2 Model dependency on input O₂

The modeled concentration and isotope profiles for the ETSP, unlike in the AS and ETNP, reflected an underestimation of water column denitrification in the best-fit model. In ETSP measurements, there is a clear deficit in [${\mathrm{NO}}_{\mathrm{3}}^{-}$], coincident with the secondary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maximum and N^* minimum (Fig. 4). In our modeled profiles, this ${\mathrm{NO}}_{\mathrm{3}}^{-}$ deficit is missing, and although a secondary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maximum is present, its magnitude is lower than observed (Fig. 4). The model also does not capture the negative N^* excursion (Fig. 4), which we think reflects a model underestimation of NAR and NIR in the ETSP. The cause of this missing denitrification is likely to be poor representation of the ETSP O₂ conditions in the model grid space. Since our model grid is fairly coarse ($\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$), only a few boxes within the ETSP had averaged [O₂] below the thresholds that would allow processes such as NAR and NIR to occur. The anoxic region of the ETSP is adjacent to the coast and not as spatially extensive as in the AS and ETNP (Fig. S8); therefore, this region in particular was less compatible with the model grid. In order to test whether the parameterization of O₂ dependence was the cause of the low N loss, we ran the model using the globally optimized parameters (Table 3) but with higher O₂ thresholds (15 µM) for NAR, NIR, and AMX (Table 1). This extended the region over which ODZ processes could occur and resulted in an increase in water column N loss from 6 to 32 Tg N yr⁻¹ in the ETSP, which is more consistent with previous estimates (DeVries et al., 2012; Deutsch et al., 2001). This change also stimulated the development of a ${\mathrm{NO}}_{\mathrm{3}}^{-}$ deficit, larger secondary ${\mathrm{NO}}_{\mathrm{2}}^{-}$ maximum, and N^* minimum within the ODZ (Fig. 7).

Figure 7 Plots of DIN concentrations and N^* from the ETSP ODZ comparing modified O₂ thresholds for N loss. In the original optimized version of the model, there is insufficient N loss and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulation in the ETSP. To demonstrate that this issue may be caused by the parameterization of O₂ and dependence and/or model grid size in the ETSP, we raised the O₂ thresholds for N loss-related processes (NAR, NIR, and AMX) to 15 µM. This effectively lowers the observed [O₂] in order to stimulate N loss. The resulting (a) [${\mathrm{NO}}_{\mathrm{3}}^{-}$], (b) [${\mathrm{NO}}_{\mathrm{2}}^{-}$], and (c) N^* are shown with the observed values from the database (yellow), original optimized model values (black), and lowered O₂ threshold model values (red).

As previously mentioned (Sect. 3.1), modeled [${\mathrm{NO}}_{\mathrm{2}}^{-}$] in the Bay of Bengal is higher than observations. The accumulation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ here in the model is likely due to O₂ concentrations falling below the set threshold for NAR but above the threshold for NIR, so ${\mathrm{NO}}_{\mathrm{2}}^{-}$ can accumulate via NAR but cannot be consumed via NIR. Although AMX and NXR occur there, the modeled rates of their ${\mathrm{NO}}_{\mathrm{2}}^{-}$ consumption are rather low, which supports a higher accumulation of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the steady-state model. This is in contrast to observations that ${\mathrm{NO}}_{\mathrm{2}}^{-}$ production is tightly matched with ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation in the Bay of Bengal, which limits ${\mathrm{NO}}_{\mathrm{2}}^{-}$ accumulation and N loss there (Bristow et al., 2017). The fact that the model over predicts NAR in the AS may also be connected with overprediction of NAR in the Bay of Bengal. It is possible that the oxygen thresholds for ODZ processes are not the same in all ODZs, and further work on oxygen sensitivities of N cycle processes will be addressed in a companion study (Martin et al., 2019).

5 Conclusions

A global inverse ocean model was modified to include ¹⁴N and ¹⁵N in both ${\mathrm{NO}}_{\mathrm{3}}^{-}$ and ${\mathrm{NO}}_{\mathrm{2}}^{-}$ as state variables. Adding the processes required to describe the cycling of ${\mathrm{NO}}_{\mathrm{2}}^{-}$ in the global ocean, including oxic and anoxic processes, resulted in a globally representative distribution of [${\mathrm{NO}}_{\mathrm{3}}^{-}$], [${\mathrm{NO}}_{\mathrm{2}}^{-}$], ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$, and ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$. In particular, the patterns of variation in both oxic and anoxic waters are generally consistent with observations, though some magnitudes of variation were somewhat muted by the model. This could be due to an underestimation of a process rate, due to parameterization or model resolution, or an underestimation of the isotope effect involved.

Importantly, we were able to generate a roughly balanced steady-state ocean N budget without the need for an artificial restoring force. The [${\mathrm{NO}}_{\mathrm{3}}^{-}$] and [${\mathrm{NO}}_{\mathrm{2}}^{-}$] distributions that were required to achieve this roughly balanced budget are well within the range of observed values. Some interesting take-home messages from this work are the following: (1) a relatively low isotope effect for ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction (ε_NAR=13 ‰) gives a good fit to ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{3}}}$ data, similar to that concluded in some recent studies (Marconi et al., 2017; Bourbonnais et al., 2015; Casciotti et al., 2013); (2) low O₂ half-saturation constants for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation allowing ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation to occur in parallel with ${\mathrm{NO}}_{\mathrm{3}}^{-}$ reduction, ${\mathrm{NO}}_{\mathrm{2}}^{-}$ reduction, and anammox were needed to achieve reasonable distributions of ${\mathrm{NO}}_{\mathrm{3}}^{-}$, ${\mathrm{NO}}_{\mathrm{2}}^{-}$, and their isotopes in the ocean water column ODZs.

Though we have been able to adequately represent and assess N cycling in ODZs, there are many areas in which this model could be improved in order to expand its usefulness. Improving resolution of the model, particularly in coastal regions where there are steep gradients in nutrient and O₂ concentrations, would improve the accuracy of the model in regions such as the ETSP. Further, in regions that have high seasonal or interannual variability, an annually averaged steady-state model may not represent some important temporal dynamics. While we attempted to account for seasonal variation in the strength of the ODZs through use of monthly O₂ climatologies, we did not simulate seasonal variations in net primary production and the strength of the biological pump. Variations in these parameters are likely to drive variations in N loss (Kalvelage et al., 2013; Ward, 2013; Babbin et al., 2014).

In addition to the dependency on external static nutrient and parameter fields, this N cycle model is highly dependent on isotope effects for N cycle processes. Previous work has shown that the laboratory-derived isotope effects for some N cycle processes are not the same as their expressed isotope effects in environmental samples or under conditions relevant to environmental samples (Casciotti et al., 2013; Bourbonnais et al., 2015; Buchwald et al., 2015; Kritee et al., 2012; Marconi et al., 2017). Further probing the isotope effects using an inverse model such as this could provide insight into the expressed isotope effects that should be used in other modeling efforts involving field data. As presented in Sect. 3.1, the larger magnitude isotope effect for ${\mathrm{NO}}_{\mathrm{2}}^{-}$ oxidation best fit the ETNP and AS ODZs, where most of the ODZ volume resides. However, most model ${\mathit{\delta }}^{\mathrm{15}}{\mathrm{N}}_{{\mathrm{NO}}_{\mathrm{2}}}$ values still do not reach the lowest values observed on the edges of marine ODZs, indicating that further work is needed to understand the expression of these isotope effects.

The larger isotope effects resulted in better fits to observations of δ¹⁵N and DIN concentrations with lower rates of N cycling. This reinforces the importance of obtaining realistic isotope effect estimates for each process that are relevant on an environmental scale. Additionally, this highlights the need for critical consideration of isotope effects used in N cycle models that use isotope balance to predict N cycling rates. Though isotopes provide us with a useful tool to assess the relative contributions of different processes, these estimates are highly subject to the isotope effects employed. Also, as illustrated by the regional optimizations, the isotope effect for a given process may vary, or be expressed differently, in different regions.

This model provides an excellent framework for further testing hypotheses about controls on the marine N inventory and cycling of N on a global scale. The distribution and sensitivities of N cycle rates resulting from this model will be explored in a companion manuscript (Martin et al., 2019). Incorporation of variable environmental input data, such as temperature, productivity, and [O₂], could also help us predict how the N cycle might be affected by past and future environmental changes.