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

Nitrite (NO−2 ) is a key intermediate in the marine nitrogen (N) cycle and a substrate in nitrification, which produces nitrate (NO−3 ), as well as water column N loss processes denitrification and anammox. In models of the marine N cycle, NO−2 is often not considered as a separate state variable, since NO−3 occurs in much higher concentrations in the ocean. In oxygen deficient zones (ODZs), however, NO−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 NO−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 NO−2 as a reactive intermediate as well as the processes that produce and consume NO−2 in marine ODZs. NO − 2 accumulation in ODZs is accurately represented by the model involving NO−3 reduction, NO − 2 reduction, NO − 2 oxidation, and anammox. We model both 14N and 15N and use a compilation of oceanographic measurements of NO−3 and NO − 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 NO−2 oxidation, and we find that the larger (more negative) inverse isotope effects for NO−2 oxidation, along with relatively high rates of NO−2 , oxidation give a better simulation of NO − 3 and NO−2 concentrations and isotopes in marine ODZs.


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 (NO − 3 ), a dissolved inorganic species, which can be taken up by microbes for use in assimilatory or dissimilatory processes.Another dissolved inorganic species, nitrite (NO − 2 ), accumulates in much lower concentrations but is a key intermediate in many N cycling processes.Models of the marine N cycle often include NO − 3 and NO − 2 together as a single dissolved inorganic N (DIN) pool, or exclude NO − 2 entirely (DeVries et al., 2013;Deutsch et al., 2007;Brandes and Devol, 2002).However, NO − 2 does accumulate significantly in oxygen deficient zones (ODZs) in features known as secondary NO − 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 2 ).The two main water column N loss processes, denitrification and anammox, use NO − 2 as a substrate.Denitrification involves the stepwise reduction of NO − 3 to NO − 2 and then to gaseous nitric oxide (NO), nitrous oxide (N 2 O), and N 2 .Anammox consists of the anaerobic oxidation of ammo-Published by Copernicus Publications on behalf of the European Geosciences Union.
T. S. Martin et al.: Modeling oceanic nitrate and nitrite isotopes nium (NH + 4 ) to N 2 using NO − 2 as the electron acceptor.NO − 2 is also oxidized to NO − 3 during anammox, representing an alternative fate for NO − 2 in ODZs.Indeed, NO − 2 oxidation appears to be prevalent in ODZs, with more NO − 2 oxidation occurring than can be explained by anammox alone (Gaye et al., 2013;Peters et al., 2016Peters et al., , 2018b;;Babbin et al., 2017;Buchwald et al., 2015;Casciotti et al., 2013;Martin and Casciotti, 2017).NO − 2 oxidation results in the regeneration of NO − 3 that would otherwise be converted to N 2 and lost from the system.The close coupling between NO − 3 reduction to NO − 2 and NO − 2 oxidation back to NO − 3 , represents a control valve on the marine N budget (Penn et al., 2016;Bristow et al., 2016).Where NO − 2 oxidation can outcompete NO − 2 reduction via denitrification and anammox, bioavailable N is retained.Water column N losses may occur primarily where NO − 2 oxidation rates are limited by oxygen availability.Thus, understanding the NO − 2 dynamics in ODZ waters is critical to assess the N loss occurring there.
The observed NO − 3 and NO − 2 concentrations alone do not allow us to fully characterize the N cycling processes occurring in a given region.Stable isotope measurements of NO − 3 and NO − 2 provide additional insight and constraints on marine N cycle processes.There are two stable isotopes of N: 14 N and 15 N.The isotopic ratios for a given N species, usually expressed in delta notation as δ 15 N (‰) = (( 15 N/ 14 N) sample /( 15 N/ 14 N) standard −1)×1000, are an integrated measure of the processes that have produced and consumed that N species.Each process imparts a unique isotope effect (ε (‰) = ( 14 k/ 15 k −1)×1000, where 14 k and 15 k are the first-order rate constants for the 14 N and 15 N containing molecules, respectively) that impacts the isotopic composition of the substrate and the product (Mariotti et al., 1981).In particular, NO − 2 cycling processes have distinct isotope effects, where NO − 2 reduction occurs with normal isotopic fractionation (Bryan et al., 1983;Martin and Casciotti, 2016;Brunner et al., 2013) and NO − 2 oxidation occurs with an unusual inverse kinetic isotope effect (Casciotti, 2009;Buchwald and Casciotti, 2010;Brunner et al., 2013).Thus, the isotopes of NO − 2 are sensitive to the relative importance of NO − 2 oxidation and NO − 2 reduction in NO − 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 NO − 2 oxidizers indicate an N isotope effect of 15 ε = −10 to −20 ‰ (Casciotti, 2009;Buchwald and Casciotti, 2010), while measured concentrations and isotopes of NO − 3 and NO − 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 NO − 2 and its isotopes as tracers.The addition of NO − 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 NO − 3 and NO − 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 -NO − 3 reduction, NO − 2 reduction, and NO − 2 oxidation -to discern what isotope effects result in the best fit to the observations.

Inverse N cycle model overview
The model used here is a steady-state inverse model that solves for the concentration and δ 15 N of NO − 3 , NO − 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 timedependent 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 14 N and 15 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 be-low 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 NO − 3 and NO − 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 14 N and 15 N.Those parameters that resulted in modeled 14 N and 15 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.
The sensitivity analysis revealed that the modeled distribution of 15 N was very sensitive to chosen isotope effects, those parameters that control the relative rates of 15 N and 14 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.

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 (NO − 3 , NO − 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 (NO − 3 , NO − 2 , and DON) between boxes.The original version of the tracer dataassimilation 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).

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: NO − 3 and NO − 2 .We did not explicitly model ammonium (NH + 4 ) because it typically occurs in low concentrations throughout the ocean, and scarcity of data (especially δ 15 N data) would make model validation difficult.Although NH + 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 14 N and 15 N of each N species to constrain the rate parameters, two sets of governing equations were employed: one that depends on 14 N and another that depends on 15 N. Generally, the rate for 15 N processes was dependent on the rate of 14 N processes and an isotopic fractionation factor (α= 14 k/ 15 k) that is specific to each process and substrate.By solving for steadystate solutions to both 14 N and 15 N concentrations, we were able to model global distributions of [NO − 3 ], [NO − 2 ], and their corresponding δ 15 N values.

N cycle parameterization
We first describe the 14 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 15 N equations and isotope implementation will be discussed in a later section.
The governing equations for the 14 N-containing DIN (NO − 3 and NO − 2 ) and organic N (DON and PON) state variables can be written as follows:

∂ ∂t
The model is designed to represent a steady state, thus the ∂ ∂t term is 0. The J terms represent the source and sink processes for each state variable, expressed in units of mmol m −3 yr −1 and will be described in more detail below.Briefly, J dep 14 is the spatially variable deposition of NO − 3 from the atmosphere to the sea surface.In the DIN model equations, J ) and J AMX 14 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 NO − 3 reduction to be dependent on organic N as well as NO − 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 [NO − 3 ] is found in the DON governing equations, it is drawn from 2013 World Ocean Atlas annual data interpolated to the model grid.

N source processes
Atmospheric deposition and N 2 fixation are the two largest sources of N to the ocean (Gruber and Galloway, 2008) and 3 ) and nitrite (NO − 2 ).N source processes are nitrogen (N 2 ) fixation and atmospheric deposition.N sink processes are sedimentary denitrification, NO − 2 reduction (NIR), and anammox (AMX).Internal cycling processes that transform N from one species to another are solubilization, remineralization, assimilation, NO − 3 reduction (NAR), ammonia oxidation (AMO), and NO − 2 oxidation (NXR).Neither ammonia (NH 3 ) nor ammonium (NH + 4 ) are tracked in this model, since they are assumed to not accumulate.N 2 is also not explicitly accounted for in the model.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 NO − 3 , and that the other species would be rapidly oxidized to NO − 3 in the oxic surface waters.(5) To calculate J dep 14 , the atmospheric deposition rate of 14 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 14 N in the deposited N (r dep 14 ), which is calculated from the isotopic composition of deposited N (δ 15 N dep , −4 ‰; Eq. 6), to yield the deposition of 14 N to the sea surface in each box (J dep 14 ).To calculate r dep 14 from δ 15 N dep , we first calculate r dep 15 using r air 15 , a standard with a value of 0.003676 (Eq.6; Mariotti, 1983).3 (Holl and Montoya, 2005) and dependence on iron (Fe) and phosphate (PO 3− 4 ) availability (Monteiro et al., 2011).
F 0 is the maximum rate of N 2 fixation (1.5 mmol m −3 yr −1 ; Table 1) and is calculated from the estimated areal rate of N 2 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 NO − 3 interpolated to the model grid (Garcia et al., 2013b).The parameter λ is an inhibition constant for N 2 fixation in the presence of NO − 3 (Table 1).The temperature (T ) terms scale the rate of N 2 fixation based on the observed temperature (T obs ), maximum observed sea surface temperature (T max ), and the minimum preferred growth temperature for Trichodesmium (T 0 ; 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 2 fixation, given the growing recognition of N 2 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 −2 yr −1 ; Chien et al., 2016) divided by the depth of the top model grid box to give units of mmol Fe m −3 yr −1 .Fe and PO 3− 4 are assumed to limit N 2 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 fix 14 , which is the fractional abundance of 14 N in newly fixed N and is calculated as in Eq. ( 6) from T. S. Martin et al.: Modeling oceanic nitrate and nitrite isotopes δ 15 N fix (−1 ‰; Table 1).All of the N 2 fixation parameters are fixed rather than optimized (Table 1).Due to the use of non-optimized parameters and an input NO − 3 field rather than modeled NO − 3 , N 2 fixation serves as an independent check that our modeled N cycle produces reasonable N concentrations and overall N loss rates.However, N 2 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 2 fixation produced by this parameterization is 131 Tg N yr −1 , which is in line with several current estimates (Table S1 in the Supplement).
In the model, N 2 fixation and NO − 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 .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): The constant has a value of −0.0101 • C −1 as determined by Dunne et al. (2005).Net primary production estimates (in units of mmol carbon m −2 yr −1 ) were taken from a satellitederived 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).

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).
J assim,WOA 14 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 [NO − 3 ] from the 2013 World Ocean Atlas annual product interpolated to the model grid [NO − 3 ] obs (Garcia et al., 2013b) to calculate the assimilation rates for DON and PON production (J assim,WOA

14
).For this assumption to be valid, our modeled surface [NO − 3 ] must be close to the observed values, which we will test in Sect.3.1.The rate constant for assimilation, 14 k assim , varies spatially and is determined using observations of surface [NO − 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 NO − 3 is being assimilated, since NO − 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 [NO − 3 ] and [NO − 2 ] rather than the World Ocean Atlas values.In order to appropriately reflect surface NO − 3 and NO − 2 concentrations, both NO − 3 and NO − 2 are assimilated. 14k assim is calculated as described above and is assumed to be the same for both NO − 3 and NO − 2 .We justify using only [NO − 3 ] to parameterize 14 k assim because NO − 3 generally makes up the bulk of DIN available for assimilation at the surface, but this assumption will be discussed in more detail below.J assim,NO 3 14

Solubilization
Solubilization is the transformation of PON to DON, and is dependent only on [PON] and a solubilization rate constant ( 14 k sol ), which is optimized (Table 2).
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 = −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 14 k sol .A spatially variable 14 k sol that accounts for lower apparent values in ODZs is a refinement that could be introduced in future model versions.

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 ( 14 k remin ), which is optimized (Table 2).
J remin 14 = 14 k remin DO 14 N (14) The removal of this remineralized DON, since it does not accumulate as NH + 4 , is either through ammonia oxidation (AMO) or anammox (AMX), depending on [O 2 ] 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 2 , NO − 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 3 ) as a substrate.Since we do not include NH 3 or NH + 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 2 dependencies of AMO and AMX.Rather than using a strict O 2 cutoff for AMO, it is limited by O 2 using Michaelis-Menten kinetics.The halfsaturation constant for O 2 , K AMO m (Table 1), sets the O 2 concentration at which AMO reaches half of its maximal value.
Recent studies have shown that AMO and NO − 2 oxidation (NXR), both O 2 -requiring processes, have very low O 2 half saturation constants and can occur down to nanomolar levels of [O 2 ] (Peng et al., 2016;Bristow et al., 2016).In contrast, O 2 -inhibited processes such as AMX are only allowed to occur at O 2 concentrations below a given threshold.The handling of O 2 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 2 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  , η 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 2 threshold for anammox.

Nitrite oxidation
The rates of NO − 2 oxidation (NXR) are dependent on the availability of NO − 2 as well as O 2 .Similar to AMO, we parameterize O 2 dependence using Michaelis-Menten kinetics and a fixed half-saturation constant for O 2 (K NXR m , Table 1).K NXR m was taken to be 0.8 µM O 2 , based on kinetics experiments performed with natural populations of NO − 2 oxidizing bacteria (Bristow et al., 2016).Finally, we employ an optimized rate constant ( 14 k NXR , Table 2) to fit the available data.

N sink processes
Nitrate and nitrite reduction NO − 3 reduction (NAR) and NO − 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 2 gas, rendering it bioavailable to only a restricted set of marine organisms.Although there are intermediate gaseous products between NO − 2 and N 2 , 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 T. S. Martin et al.: Modeling oceanic nitrate and nitrite isotopes on the presence of NO − 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 NO − 3 or NO − 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 NO − 3 ) in the linear system.Both NAR and NIR are dependent on the remineralization rate (J remin 14) 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 2 .As mentioned above, we assume that J remin 14 does not depend on the choice of electron acceptor.
The rate coefficients for NAR ( 14 k NAR ) and NIR ( 14 k NIR ) are optimized rather than fixed (Table 2).Further, the dependence of J NAR 14 and J NIR 14 on J remin 14 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 2 , like AMX, is parameterized by a parameter η, which inhibits these processes when [O 2 ] 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 2 ] was above the threshold for the process and 1 if [O 2 ] was below the threshold.On further refinement, we wanted to account for the possibility of seasonal shifts in [O 2 ] 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 2 thresholds used to calculate η NAR and η NIR were fixed (7 and 5 µM, respectively; Table 1).Since we do not explicitly model O 2 , [O 2 ] was predetermined using the 2013 World Ocean Atlas monthly O 2 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 2 ] data to observed suboxic measurements (Bianchi et al., 2012).

Anammox
Anammox (AMX) catalyzes the production of N 2 from NH + 4 and NO − 2 .Since we do not use NH + 4 as a variable in our N cycling equations, we substituted remineralized DON (J remin 14 ) as a proxy for NH + 4 availability.As described above in Sect.2.3.3,remineralized DON is routed through either AMO or AMX depending on [O 2 ] and the O 2 dependencies of AMO and AMX.
The O 2 threshold used to calculate η AMX from monthly O 2 climatology is fixed (10 µM; Table 1).In order to maintain mass balance on remineralized DON, we do not include dependence on [NO − 2 ] in Eq. ( 19), although J AMX 14 removes NO − 2 (Eq.2).This parameterization inherently assumes that AMX is limited primarily by [NH + 4 ] supply and not [NO 2 ], which may not always be correct (Bristow et al., 2016).Anammox also produces 0.3 moles of NO − 3 via associated NXR for every 1 mole of N 2 gas produced (Strous et al., 1999).For this reason, anammox appears in the state equation for NO − 3 (Eq.1).

Sedimentary denitrification
Sedimentary (or benthic) denitrification (J sed 14) 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 2 ], bottom water [NO − 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): 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 = −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 3 ] 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 2 .The linear relationships were then rearranged in order to estimate sedimentary denitrification rate as a function of RRPOC, [O 2 ], and [NO − 3 ].These equations were then further broken down into a component that is dependent on [NO − 3 ] and a component that is dependent on [O 2 ] (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 NH + 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 −1 .This transfer function also assumes that all of the NH + 4 efflux is immediately oxidized to NO − 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 NO − 3 isotopes, as several studies suggest that benthic N processes may contribute to water column nitrate 15 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 δ 15 N NO 3 .Additionally, our spatial resolution does not well represent regions where this effect might be significant on bottom water δ 15 N NO 3 , such as the shallow shelves.

N isotope implementation
In our model, we are interested in using the isotopic composition of NO − 3 and NO − 2 to constrain the rates of N cycling and loss from the global ocean.As DON and PON are ultimate substrates for NO − 2 and NO − 3 production, it is essential to track the 15 N in the organic N pools as well.The matrix setup for 15 N is similar to that for the 14 N species, but the rates were changed as follows: . (21) is the rate of each relevant 14 N process as described above, and J process 15 is the rate of each 15 N process.The α process is the fractionation factor for a given process, which is given by the ratio between the rate constants for 14 N and 15 N (α= 14 k/ 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 15 N concentrations of N species from our observations to constrain the model, we calculated 15 N/ 14 N from measured δ 15 N and multiplied by the measured concentration of each modeled N species, assuming that This simple 15 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 2 fixation, and anammox, as described below.

Atmospheric deposition
For atmospheric deposition of N, we prescribe fixed δ 15 N value of −4 ‰ (Table 1), which can be related to the fractional abundance of 14 N (r

Nitrogen fixation
Similar to atmospheric deposition, newly fixed N has a δ 15 N value (−1 ‰; Table 1).In Sect.2.3.2 we described r fix 14 , the fractional abundance of 14 N in newly fixed N.Here we multiply the fractional abundance of 15 N, r fix 15 , by the other terms in the N 2 fixation equation (Eq.6) to obtain the rate of 15 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 2 (NO − 2 and NH + 4 ), as well as for the associated NO − 2 oxidation to NO − 3 .We assume that the fractionation factor for ammonium oxidation via AMX (α AMX,NH 4 ) is 1, setting it to match the fractionation factor for AMO (α AMO ; Table 1), both with no expressed fractionation since NH + 4 does not accumulate in the model.Since  Buchwald and Casciotti (2010), Casciotti et al. (2013) all remineralized DON must be routed either through AMO or AMX, this simplifies the mass balance and ensures that all remineralized 14 N and 15 N is accounted for. 15NO − 2 is removed with the isotope effects of NO − 2 reduction (α AMX,NIR ) and NO − 2 oxidation (α AMX,NXR ) in the expected 1 : 0.3 proportion (Brunner et al., 2013).

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., [NO − 3 ] and [NO − 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 2 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 NO − 3 that are also dependent on [NO 2 ] for each model box.The top half of vector b would be independent processes that produce or consume [NO − 3 ], and the bottom half of vector b would be independent processes that produce or consume [NO − 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.

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 NO − 3 and NO − 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 14 N-organic N model, 15   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: 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 [NO − 3 ] data points = n NO 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 [NO − 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).

Global model-data comparison
The simulations of NO − 2 concentration and isotopic composition are the most unique features of this model in compari-son to existing global models of the marine N cycle.As such, NO − 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 NO − 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 2 conditions found there.However, accumulation of NO − 2 in the model ETSP was lower than expected.The model also accumulated NO − 2 in the Bay of Bengal, which is a low-O 2 region off the east coast of India that does not generally accumulate NO − 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 NO − 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 [NO − 3 ], [NO − 2 ], δ 15 N NO 3 and δ 15 N NO 2 .The best fit was achieved for isotope effects of 13 ‰ for NO − 3 reduction (ε NAR ), 0 ‰ for NO − 2 reduction (ε NIR ), and −13 ‰ for NO − 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.
In the test set, there were some low [O 2 ] points where our model [NO − 3 ] exceeded observations (Fig. 3a, filled black circles); these are largely within the ETSP.In contrast, the AS tended to show slightly lower modeled [NO − 3 ] than expected.The [NO − 2 ] accumulation (Fig. 3b) and δ 15 N NO 3 signals (Fig. 3c) in the ETSP were also generally too low compared with observations.These signals are likely tied to insufficient NO − 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 NO − 2 accumulations observed sporadically (Anderson et al., 1982;Codispoti et al., 1985Codispoti et al., , 1986)).
Overall, the representation of δ 15 N NO 3 was fairly good (RMSE = 2.4 ‰), though there were a subset of points above δ 15 N NO 3 = 10 ‰ where the modeled δ 15 N NO 3 exceeded the observed δ 15 N NO 3 , and others where modeled δ 15 N NO 3 was lower than observations (Fig. 3c).Many of the points with overestimated δ 15 N NO 3 were located within the AS ODZ, where there may be too much NO − 3 reduction occurring, leading to artificially elevated δ 15 N NO 3 values.As indicated above, the underestimated δ 15 N NO 3 points largely fell within the ETSP where we believe the model is underestimating NO − 3 reduction.The representation of δ 15 N NO 2 was also fairly good (RMSE = 8.6 ‰), though the modeled δ 15 N NO 2 was generally not low enough (Fig. 3d), indicating an underestimated sink of "heavy" NO − 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 NO − 3 and NO − 2 concentration and isotope profiles in the AS and ETNP were consistent with the observations, with [NO − 3 ] slightly underestimated in the AS ODZ and overestimated in the ETSP (Fig. 4).As [O 2 ] goes to zero, the O 2 -intolerant processes NAR, NIR, and AMX are released from inhibition.These processes result in a decrease in [NO − 3 ] (via NAR) which corresponds to an increase in δ 15 N NO 3 , since NAR has a normal isotope effect.NO − 2 also starts to accumulate in the secondary NO − 2 maximum as a result of NAR.The δ 15 N NO 2 is lower than δ 15 N NO 3 since light NO − 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 [NO − 3 ] depletion and [NO − 2 ] accumulation in the model were lower than observed.This could be due in part to the timeindependent 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).
In order to gauge the model results for N loss, we also calculated N * , a measure of the availability of DIN relative to PO 3− 4 compared to Redfield ratio stoichiometry Deutsch et al., 2001).Negative N * values are associated with N loss due to AMX or NIR or release of PO 3− 4 from anoxic sediments (Noffke et al., 2012), while positive N * values are associated with input of new N through N 2 fixation (Gruber and Sarmiento, 1997).Although we did not model PO 3− 4 , we used the modeled [NO − 3 ] and 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 ε NXR = −32 ‰.For the ETNP, the best fit isotope effects were ε NAR = 13 ‰, ε NIR = 15 ‰, and ε NXR = −32 ‰, though the performance is only marginally better than with ε NIR = 0 ‰.The lower (more inverse) value for ε NXR resulted in higher δ 15 N NO 3 and lower δ 15 N NO 2 , which better fit the ODZ δ 15 N NO 2 data compared to the global best fit ε NXR = −13 ‰.These results are consistent with earlier isotope modeling studies in the ETSP (Casciotti et al., 2013;Peters et al., 2016Peters et al., , 2018b) ) and in the AS (Martin and Casciotti, 2017).Although, in the AS, modeled δ 15 N NO 3 values were too high, likely in part due to overpredicted rates of NAR, which also resulted in lower modeled [NO − 3 ] (Fig. 4).

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 [NO − 3 ] is low in surface waters and increases to a mid-depth maximum between 1000 and 2000 m.The highest [NO − 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 NO − 3 is accumulated in the deep waters of the model Pacific.This could be due to an underestimate of preformed NO − 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 [NO − 3 ] are 5-10 µM lower than observations (Fig. S5), which could be enough to explain the lower-than-expected [NO − 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 δ 15 N NO 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 NO − 3 and increase in δ 15 N NO 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 δ 15 N NO 3 in the www.biogeosciences.net/16/347/2019/Biogeosciences, 16, 347-367, 2019 ETSP ODZ and the upper thermocline in the eastern part of the section is consistent with an underestimate of NO − 3 reduction.In GP16 we were also able to compare modeled and observed [NO − 2 ] and δ 15 N NO 2 (Fig. S6).Patterns of modeled [NO − 2 ] and δ 15 N NO 2 showed accumulation of NO − 2 in the ODZ, with an appropriate δ 15 N NO 2 value (Fig. S6).Although, generally lower modeled concentrations of NO − 2 in the ODZ also support an underestimate of NAR (Fig. S6).
Surface δ 15 N NO 3 values were also not as high in the model as in the observations (Fig. 5), which could result from insufficient NO − 3 assimilation or too low supplied δ 15 N NO 3 (Peters et al., 2018a).However, we do see a similar depth range for high surface δ 15 N NO 3 and a local δ 15 N NO 3 minimum between the surface and ODZ propagating westward in both the model and observations, indicating that the physical and biogeochemical processes affecting δ 15 N NO 3 are represented by the model.Additionally, the model shows slightly elevated δ 15 N NO 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 δ 15 N NO 3 across the basin.Peters et al. (2018a) and Rafter et al. (2013) also postulated that these elevated δ 15 N NO 3 values were in part driven by remineralization of organic matter with high δ 15 N.The δ 15 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 δ 15 N NO 3 in the intermediate depths (500-1500 m), which is consistent with observations, again reflecting remineralization of PON with δ 15 N greater than mean ocean δ 15 N NO 3 .Overall, the patterns of δ 15 N NO 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 δ 15 N NO 3 values that are lower than observations.The simplification of NH + 4 dynamics in the model could also contribute to underestimation of δ 15 N NO 3 values if there was a large flux of 15 N-enriched NH + 4 from sediments (Granger et al., 2011), or if 15 N-depleted NH + 4 was preferentially transferred to the N 2 pool via anammox.While the isotope effect on NH + 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 NH + 4 oxidation (Table 1).
In the North Atlantic along GEOTRACES section GA03, we see good agreement between the observed and modeled [NO − 3 ] (Fig. 6).There is generally low surface [NO − 3 ] with a distinct area of high [NO − 3 ] propagating from near the African coast.Deep water (> 2000 m) [NO − 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 NO − 3 present in Southern Ocean-sourced intermediate waters (500-1500 m; Marconi et al., 2015).Modeled δ 15 N NO 3 values at first glance appear higher than observed values at the surface (Fig. 6).However, many of the surface [NO − 3 ] were below the operating limit for δ 15 N NO 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 δ 15 N NO 3 values in upper thermocline waters in both the model and observations, likely corresponding to low δ 15 N contributions from N 2 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 2 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 δ 15 N signal (Knapp et al., 2008).In our model, atmospheric N deposition contributed between 0 % and 50 % of N input along the cruise track.

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 [NO − 3 ] to estimate the contribution of DIN assimilation to the production of organic N, whereas the DIN model uses modeled [NO − 3 ] and [NO − 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 NO − 2 in the top two boxes of the DIN model, since NO − 2 assimilation is unaccounted for in the organic N model.This is largely an issue in the oligotrophic gyres, where surface [NO − 3 ] is very low and NO − 2 accumulates to low but nonzero concentrations (Fig. 2).Assimilation of NO − 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 [NO − 3 ] is higher than the World Ocean Atlas surface [NO − 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 mod-

eled [NO −
3 ] that was higher than observed [NO − 3 ] (Fig. S7).Likewise, points with relatively lower DIN assimilation had modeled [NO − 3 ] less than observed [NO − 3 ].However, the majority of DIN assimilation estimates were within 10 µM yr −1 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 [NO − 3 ] and modeled [NO − 3 ].However, we find that the World Ocean Atlas surface NO − 3 values are fairly well represented by our modeled surface NO − 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.

Model dependency on input O 2
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 [NO − 3 ], coincident with the secondary NO − 2 maximum and N * minimum (Fig. 4).In our modeled profiles, this NO − 3 deficit is missing, and although a secondary NO − 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 2 conditions in the model grid space.Since our model grid is fairly coarse (2 • × 2 • ), only a few boxes within the ETSP had averaged [O 2 ] 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 2 dependence was the cause of the low N loss, we ran the model using the globally optimized parameters (Table 3) but with higher O 2 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 −1 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 NO − 3 deficit, larger secondary NO − 2 maximum, and N * minimum within the ODZ (Fig. 7).
As previously mentioned (Sect.3.1), modeled [NO − 2 ] in the Bay of Bengal is higher than observations.The accumulation of NO − 2 here in the model is likely due to O 2 concentrations falling below the set threshold for NAR but above the threshold for NIR, so NO − 2 can accumulate via NAR but cannot be consumed via NIR.Although AMX and NXR occur there, the modeled rates of their NO − 2 consumption are rather low, which supports a higher accumulation of NO − in the steady-state model.This is in contrast to observations that NO − 2 production is tightly matched with NO − 2 oxidation in the Bay of Bengal, which limits NO − 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).

Conclusions
A global inverse ocean model was modified to include 14 N and 15 N in both NO − 3 and NO − 2 as state variables.Adding the processes required to describe the cycling of NO − 2 in the global ocean, including oxic and anoxic processes, resulted in a globally representative distribution of [NO − 3 ], [NO − 2 ], δ 15 N NO 3 , and δ 15 N NO 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 [NO − 3 ] and [NO − 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 NO − 3 reduction (ε NAR = 13 ‰) gives a good fit to δ 15 N NO 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 2 half-saturation constants for NO − 2 oxidation allowing NO − 2 oxidation to occur in parallel with NO − 3 reduction, NO − 2 reduction, and anammox were needed to achieve reasonable distributions of NO − 3 , NO − 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 2 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 2 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 NO − 2 oxidation best fit the ETNP and AS ODZs, where most of the ODZ volume resides.However, most model δ 15 N NO 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 δ 15 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 2 ], could also help us predict how the N cycle might be affected by past and future environmental changes.

Figure 1 .
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 (NO −3 ) and nitrite (NO − 2 ).N source processes are nitrogen (N 2 ) fixation and atmospheric deposition.N sink processes are sedimentary denitrification, NO − 2 reduction (NIR), and anammox (AMX).Internal cycling processes that transform N from one species to another are solubilization, remineralization, assimilation, NO − 3 reduction (NAR), ammonia oxidation (AMO), and NO − 2 oxidation (NXR).Neither ammonia (NH 3 ) nor ammonium (NH + 4 ) are tracked in this model, since they are assumed to not accumulate.N 2 is also not explicitly accounted for in the model.
14 N = 15 N/( 15 N+ 14 N), we calculate r dep 14 as (1 − r dep 15 ).The units of S dep are given in mg N m −2 yr −1 , which we convert to mmol NO − 3 m −3 yr −1 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 2 fixation N 2 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 2 fixation in the model of De-Vries et al. (2013), with partial inhibition by NO − rate on [O 2 ] − [NO − 3 ].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 [O 2 ] − [NO − 3 ] and sedimentary denitrification rate (Fig. S2).We obtained three linear relationships between [O 2 ]−[NO − 3 ] and sedimentary denitrification rate, each applicable across a given range of [O 2 ] − [NO − 3 ] 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 2 rather than [O 2 ]−[NO − 3 ].Therefore, a linear relationship between O 2 and [O 2 ]−[NO − dep , the estimated rate of total N deposition to obtain J dep 15 .
N-organic N model, 14 N-DIN model, and 15 N-DIN model.The modeled output [NO − 3 ], [NO − 2 ], δ 15 N NO 3 , and δ 15 N NO 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 NO − 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

Figure 2 .
Figure 2. Map showing the model-estimated accumulation of nitrite (NO − 2 ) at 200 m depth.

Figure 4 .
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 ([NO − 3 ]), nitrite concentration ([NO − 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 δ 15 N NO 3 and δ 15 N NO 2 are shown for three different combinations of isotope effect.The blue lines represent ε NAR = 13, ε NXR = −13, and ε NIR = 0, which are the best fit isotope effects globally and in the ETSP.The red lines represent ε NAR = 13, ε NXR = −32, and ε NIR = 0, which are the best fit isotope effects in the AS.

Figure 5 .
Figure 5. Section profiles of NO − 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 [NO − 3 ] to (b) modeled [NO − 3 ] is presented over the full depth range (0-6000 m).Comparison of (c) observed δ 15 N NO 3 to (d) modeled δ 15 N NO 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).

Figure 6 .
Figure 6.Section profiles of NO − 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 [NO − 3 ] to (b) modeled [NO − 3 ] is presented over the full depth range (0-6000 m).Comparison of (c) observed δ 15 N NO 3 to (d) modeled δ 15 N NO 3 is presented over a shortened depth range (0-1000 m) to better assess surface and the low δ 15 N NO 3 contribution from N 2 fixation.GEOTRACES data are from Marconi et al. (2015) and available from BCO-DMO.

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