Interactive comment on “ Modeling oceanic nitrite concentrations and isotopes using a 3 D inverse N cycle model ”

Abstract. Nitrite (NO2−) is a key intermediate in the marine nitrogen (N) cycle and a substrate in nitrification, which produces nitrate (NO3−), as well as water column N loss processes, denitrification and anammox. In models of the marine N cycle, NO2− is often not considered as a separate state variable, since NO3− occurs in much higher concentrations in the ocean. In oxygen deficient zones (ODZs), however, NO2− represents a substantial fraction of the bioavailable N, and modeling its production and consumption is important to understanding the N cycle processes occurring there, especially those where bioavailable N is lost from or retained within the water column. Here we present the expansion of a global 3D inverse N cycle model to include NO2− as a reactive intermediate as well as the processes that produce and consume NO2− in marine ODZs. NO2− accumulation in ODZs is accurately represented by the model involving NO3− reduction, NO2− reduction, NO2− oxidation, and anammox. We model both 14N and 15N and use a compilation of oceanographic measurements of NO3− and NO2− 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 NO2− oxidation, and we find that the larger (more negative) inverse isotope effects for NO2− oxidation along with relatively high rates of NO2− oxidation give a better simulation of NO3− and NO2− concentrations and isotopes in marine ODZs.



Introduction
Nitrogen (N) is an important nutrient to consider when assessing the 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 of nitrogen 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, 25 2008).
There are several chemical species in which N can be found in the ocean.The largest pool of bioavailable N is nitrate (NO3 -), a dissolved inorganic species, which can be taken up by microbes for use in assimilatory or dissimilatory processes.Another Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.effects used for given processes.Though there are estimates of isotope effects for processes based on both environmental measurements and in-lab studies, there is not always agreement between them.For example, laboratory cultures of NO2 - oxidizers indicate an N isotope effect of 15 ε = -13‰ (Casciotti, 2009), while measured concentrations and isotopes of NO3 - and NO2 -in ODZs indicate that an isotope effect closer to -30‰ is more consistent with the observations (Buchwald et al., 2015;Casciotti et al., 2013).5 Here we present an expansion of an existing global ocean 3D inverse isotope-resolving N cycling model (DeVries et al., 2013) to include NO2 -and its isotopes as tracers.The addition of NO2 -allows us to include many additional internal N cycling processes, as well as a more nuanced and realistic version of the processes occurring in ODZs.This should result in a reasonably accurate global representation of NO3 -and NO2 -concentrations and isotopes.We used a database of NO3 -and NO2 - 10   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-NO3 -reduction, NO2 -reduction, and NO2 -oxidation-to discern what isotope effect estimates result in the best fit to the observations.2 Methods 15

Inverse nitrogen cycle model overview
The model used here is a steady-state inverse model that solves for the concentration and δ 15 N of NO3 -, NO2 -, 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.However, on inter-annual timescales the N cycle is thought to be approximately in balance (Gruber, 2004;Bianchi et al., 2012).The residence time of 20 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 is 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 25 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.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.

30
The model's uncertain biological parameters were determined through an optimization process that minimizes the difference between the modeled and observed NO3 -and NO2 -concentration and isotope data.Computational time limits the number of Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.parameters that we were able to optimize.We therefore focus 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 exhibiting 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 nitrogen in the ocean in 5 general and in the ODZs in particular.
The sensitivity analysis revealed that the modeled distribution of 15 N was very sensitive to isotope effects, parameters that control the relative rates of 15 N and 14 N in chemical and biological processes.There are literature estimates each of the isotope effects of interest in this work, although there is often a discrepancy between isotope effects estimated in laboratory studies 10 and those expressed in oceanographic measurements (Casciotti et al., 2013;Buchwald et al., 2015;Bourbonnais et al., 2015;Martin and Casciotti, 2017;Fuchsman 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.15 In addition to the optimized parameters and isotope effects, there were some non-sensitive 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, 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 20
The model uses a uniform 2˚x2˚ 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-minute gridded bathymetry (ETOPO2v2) that was then interpolated to the model grid.Our linear N cycle model relies on the transport of dissolved N species (NO3 -, NO2 -, 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 N species (NO3 -, NO2 -, and DON) between boxes.The 25 original version of the tracer data-assimilation procedure used to generate the transport operator for dissolved species (Tf) 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.There are two organic N (ON) pools: DON and PON.
There are also two dissolved inorganic N (DIN) pools: NO3 -and NO2 -.We excluded ammonium (NH4 + ) from the model 30 because it typically only accumulates in low concentrations throughout the ocean, and scarcity of data (especially 15  Because we are interested in both 14 N and 15 N concentrations 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 α that is specific to each process and substrate.By solving for steady-state solutions to both 14 N and 15 N concentrations, we were able to model global distributions 5 of [NO3 -], [NO2 -], and their corresponding δ 15 N values.

N cycle parameterization
We will first address 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.10 The governing equations for the 14 N-containing ON and DIN state variables can be written as follows: 123 − J -. 56678,:;+ − J -.

6LM
The model is assumed to be in steady state, thus the # #$ term is 0. The J terms represent the source and sink processes for each state variable, expressed in units of mmol/m 3 /yr and will be described in more detail below.Briefly, J -. 123 is the spatiallyvariable deposition of NO3 -from the atmosphere to the sea surface.J -. 56678,:;+ and J -. 56678,:;C represent the assimilation of NO3 - 20 and NO2 -, respectively, by phytoplankton in the upper two box levels.This assimilated NO3 -produces DON and PON, with proportions set by a spatially-variable s term.Assimilation in the DON and PON equations is represented by J -. 56678,J;< and is dependent on WOA [NO3 -] as described below.N2 fixation is split between DON and PON with the same s term.NO3 - reduction (J -. :<= ), NO2 -reduction (J -. :D= ), NO2 -oxidation (J -. :>= ), and anammox (J -. <B> ) act on the NO3 -and NO2 -pools.J -.

621
represents the removal of NO3 -via benthic denitrification.J -. 6LM represents the dissolution of PON into DON.J -. N287O is related 25 to ammonia oxidation (J -. <B; ) and J -. <B> as described below.
Through the use of these J terms, the governing equations are all linear with respect to the state variables.In order to introduce dependence of rates on the concentrations of multiple state variables, we run the ON equations first and the DIN equations annual data interpolated to the model grid.

Atmospheric deposition
Deposition of N to the ocean from the atmosphere is one of two sources of bioavailable N supplied to the ocean in this model.5 N deposition is assumed to only occur in the top box of the model, and we assume that most of the N deposited is as NO3 -, and that the other species would be rapidly oxidized to NO3 -in the oxic surface waters.
VWX as (1 - -U VWX ).The units of S 123 are given in mg N/m 2 /yr, which we convert to mmol NO3 -/m 3 /yr via dimensional analysis and by dividing by the depth of the surface box.This source term of N to the model is independent of the modeled N terms.depth of the top model grid box to give units of mmol Fe/m 3 /yr.Fe and PO4 3-are assumed to limit N2 fixation at low concentrations via Michaelis-Menten kinetics.KFe and KP are their respective half-saturation constants.Additionally, there is a term that allows us to set the isotopic ratio of newly fixed N, r -.

N2 fixation
'7I , which is the fractional abundance of 14 N in newly fixed N and is calculated as in Equation 6 from δ 15 Nfix (-1‰; Table 1).All of the N2 fixation parameters are fixed rather than optimized 10 (Table 1).Due to the use of non-optimized parameters and an input NO3 -field rather than modeled NO3 -, N2 fixation serves as an independent check that our modeled N cycle produces reasonable N concentrations and overall N loss rates.However, N2 fixation is not explicitly modeled here and is instead taken as a fixed, though spatially variable, input field.
In the model, N2 fixation and NO3 -assimilation (Section 2.3.3) are assumed to be the two processes that create exportable OM. 15 A fraction, σ, of this OM is portioned into DON rather than PON (Equations 3-4).In order to create spatial variability in this constant, we assumed (1-σ), the fraction of assimilated N partitioned to PON, is equal to the particle export (pe) ratio.This pe ratio is the ratio of particle export to primary production, and is equivalent to the fraction of OM that is exported from the euphotic zone as particulate matter rather than recycled or solubilized into DON.The pe ratio is calculated for each model grid square from the mixed layer temperature (Tml) and net primary production (NPP) as described by Dunne et al. (2005): 20 8. pe = ϕT 8M + 0.582 log(NPP) + 0.419 ɸ has a constant value of -0.0101 ˚C-1 .NPP estimates in units of mmol carbon/m 2 /yr were taken from a satellite-derived productivity model (Westberry et al., 2008), annually averaged, and interpolated onto the model grid.Tml is calculated from the 2013 WOA 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 Tml.As temperature increases, the pe ratio decreases and less PON is exported, 25 resulting in more DON recycling in the surface with several possible explanatory mechanisms discussed in greater detail by Dunne et al. (2005).As NPP increases, the pe ratio increases and more PON is exported; NPP explains 74% of the observed variance in particle export (Dunne et al., 2005).

Assimilation of nitrate and nitrite 30
Assimilation accounts for the uptake of DIN and its incorporation into OM.Since assimilation affects both the ON and DIN pools, we must account for it in both sets of model runs.Since the ON model is run first and the assimilation rates are dependent Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.on DIN concentrations, assumptions must be made about the DIN field in order to account for assimilation prior to the DIN model runs.We will first address assimilation in the ON equations.
The assimilation rates for DON and PON must be calculated using observed surface [NO3 -], rather than modeled [NO3 -].For this assumption to be valid, our modeled surface [NO3 -] must be close to the observed values, which we will test in Section 5 3.1.We also assume from the perspective of ON that only NO3 -is being assimilated, since NO2 -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 pe ratio as previously described and shown in Equations 3 and 4. 9.  -.
, varies spatially and is determined using observations of surface [NO3 -] and satellite derived net primary production (NPP) estimates (Westberry et al., 2008).The rate constant is converted to N units using the ratio of carbon (C) to N in organic matter (rC:N) which we assume to be 106:16 (Redfield et al., 1963).The value of the rate constant is only non-zero in the top two boxes of the model, where we assume primary production to be occurring.The same rate 15 constant is used in both the ON and DIN assimilation equations.
The setup for assimilation in the DIN equations is similar, but can use modeled [NO3 -] and [NO2 -] rather than the WOA values.
is calculated as described above and is assumed to be the same for both NO3 -and NO2 -.We justify using only [NO3 -] to 20 parameterize  ‚ƒƒ"… -.
because NO3 -generally makes up the bulk of DIN available for assimilation at the surface, but this assumption will be discussed in more detail below.

Solubilization
Solubilization is the transformation of PON to DON, and is dependent only on [PON] and a solubilization rate constant ( 14 ksol), which is optimized (Table 2).
[  -. ] The solubilization of PON, together with the particle transport operator (Tp), produces a particle flux attenuation curve similar 30 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

Remineralization
Remineralization, or ammonification, is the release of DON into the DIN pool.This is determined using the concentration of 5 DON and a remineralization rate constant ( 14 kremin), which is optimized (Table 2).

š𝐷𝑂 𝑁 -. oe
The removal of this remineralized DON, since it does not accumulate as NH4 + , is either through ammonia oxidation (AMO) or anammox (AMX), depending on [O2] as described below and in Section 2.3.4.We use the same remineralization rate constant regardless of the utilized electron acceptor (e.g.O2, NO3 -).Since particle flux attenuation is observed to be somewhat 10 weaker in oxygen deficient zones compared with oxygenated water (Van Mooy et al., 2002), this may slightly overestimate the rates of heterotrophic remineralization occurring in ODZs.

Ammonia oxidation
Ammonia oxidation (AMO) uses ammonia (NH3) as a substrate.Since we do not include NH3 or NH4 + in the model system, 15 we treat remineralized DON as the substrate for AMO.In order to maintain consistency between the ON and DIN model runs, remineralized DON is routed either to AMO or AMX (lost from the system) based on the O2 dependencies of AMO and AMX.
Rather than using a strict O2 cutoff for AMO, it is limited by O2 using Michaelis-Menten kinetics.The half-saturation constant for O2, Km AMO (Table 1), sets the O2 concentration at which AMO reaches half of its maximal value.
Recent studies have shown that AMO and NO2 -oxidation (NXR), both O2-requiring processes, have very low O2 half saturation constants and can occur down to nM levels of [O2] (Peng et al., 2015;Bristow et al., 2016b).In contrast, O2-inhibited processes such as AMX are only allowed to occur at O2 concentrations below a given threshold.The handling of O2 thresholds for anaerobic processes is discussed in more detail below (Section 2.3.4).Briefly, the O2 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 25 year its WOA 2013 falls below the [O2] threshold for anammox (O2 AMX , Table 1).If, for example, [O2] in a given grid box is always above the threshold for AMX, ηAMX = 0 and all of the remineralized DON (represented by J14 rem ) will be oxidized via AMO.If [O2] is less than O2 AMX , ηAMX will be non-zero 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 O2 threshold for anammox.30 Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.

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

N sink processes
Nitrate and nitrite reduction NO3 -reduction (NAR) and NO2 -reduction (NIR) are two processes within the stepwise reductive pathway of canonical denitrification.The end result of denitrification is the conversion of DIN as N2 gas, rendering it bioavailable to only a restricted 10 set of marine organisms.Although there are intermediate gaseous products between NO2 -and N2, we treat NO2 -reduction 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.Since both NAR and NIR are heterotrophic processes, they consume organic matter in addition to their main N substrates/electron receptors.In order to maintain levels of 15 NAR and NIR that are bounded both by the available NO3 -or NO2 -and the available organic matter in a linear model, it was necessary to run ON and DIN equations separately, since it is not possible to include dependencies on two state variables (e.g.DON and NO3 -) in the linear system.Both NAR and NIR are dependent on the remineralization rate ( -. •W…"™ ) that is calculated in the ON model run.In model boxes where NAR and NIR are occurring, some of the remineralization is carried out with electron acceptors other than O2.For simplicity, we assume that  -.
•W…"™ means that kNAR and kNIR are not first order rate constants and have different units than kPON, kDON, and kNXR (Table 2).25 The inhibition of NAR and NIR by O2, like AMX, is parameterized by a parameter η, which inhibits these processes when [O2] is above their maximum threshold.Originally, we treated this term as a binary operator where it would be 0 if the empirically-corrected 2013 WOA annually averaged [O2] is above the threshold for the process and 1 if [O2] is below the threshold.On further refinement, we wanted to account for the possibility of seasonal shifts in [O2] in ODZs.Thus, for each 30 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 (7 µM and 5 µM, respectively; Table 1).Since we do not explicitly model O2, [O2] was predetermined using the 2013 WOA monthly O2 climatology (Garcia et al., 2014) interpolated to the model grid.We also applied an empirical correction that improves the fit of WOA [O2] data to observed suboxic measurements (Bianchi et al., 2012).

Anammox 5
Anammox (AMX) catalyzes the production of N2 from NH4 + and NO2 -.Since we do not use NH4 + as a variable in our N cycling equations, we substituted remineralized DON ( -. •W…"™ ) as a proxy for NH4 + availability.As described above in Section 2.3.3,remineralized DON is routed through either AMO or AMX depending on [O2] and the O2 dependencies of AMO and AMX.

Sedimentary denitrification
Sedimentary (or benthic) denitrification ( -. ƒWV ) is an important loss term for N in marine sediments, 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.) OE

25
NPP is derived from the productivity modeling of Westberry et al. (2008) as described in Section 2.3.2.The pe ratio is calculated as previously described in Section 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 2nd 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).30 The transfer function for sedimentary denitrification was originally described using a non-linear dependence of the rate on In organic-rich shelf sediments an additional term is introduced that reduces the sedimentary denitrification rate by a globally averaged 27% if the depth of the bottom model box is less than 1000m due to the potential for efflux of NH4 + into the water column (Bohlen et al., 2012).This decreases overall sedimentary denitrification by approximately 6 Tg N/yr.This 10 transfer function also assumes that all of the NH4 + efflux is immediately oxidized to NO3 -and does not alter its isotopic composition in bottom water.

N isotope implementation
In our model, we are interested in using the isotopic composition of NO3 -and NO2 -to constrain the rates of N cycling and loss from the global ocean.As DON and PON are ultimate substrates for NO2 -and NO3 -production, it is essential to track the 15  N 15 in the ON pools as well.The matrix setup is similar to that for the 14 N species, but the rates were changed as follows: •‹´Wƒƒ is the rate of each relevant 14 N process as described above, and  -U X•‹´Wƒƒ is the rate of each 15 N process. X•‹´Wƒƒ 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 20 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 (Section 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  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 30 isotopic parameterizations were also required for atmospheric deposition, N2 fixation and anammox, as described below.

Atmospheric deposition
For atmospheric deposition of N, we prescribe a constant δ 15 N value of -4‰ (Table 1), which can be related to the fractional abundance of 14 N, previously described in Section 2.3.2 as  -.
VWX , and the fractional abundance of 15 N ( -U VWX ) in deposited N.
We multiply  -U VWX by S dep , the estimated rate of total N deposition, to obtain  -U VWX . 5

Nitrogen fixation
Similar to atmospheric deposition, newly fixed N has a constant prescribed δ 15 N value (-1‰; Table 1).In Section 2.3.2 we described  -.
¹"º , the fractional abundance of 14 N in newly fixed N.Here we multiply the fractional abundance of 15 N,  -U ¹"º , by the other terms in the N2 fixation equation (Equation 6) to obtain the rate of 15 N fixation. 10

Anammox
Anammox is the most complicated to parameterize isotopically because it has three different isotope effects associated with it.
There is an isotope effect on both substrates converted to N2 (NO2 -and NH4 + ), as well as for the associated NO2 -oxidation to NO3 -.We assume that the fractionation factor for ammonium oxidation via AMX (αAMX,NH4) is 1, setting it to match the fractionation factor for AMO (αAMO; Table 1), both with no expressed fractionation since NH4 + does not accumulate in the 15 model.Since 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. 15NO2 -is removed with the isotope effects of NO2 -reduction (αAMX,NIR) and NO2 -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 20 were of the general form Ax = b.All 200,160 model boxes are accounted for in the matrices.Matrix A (400,320 x 400,320) contained the rate constants and other parameters that are multiplied by the vector of state variables, x (400,320 x 1).Vector x contained the state variables (i.e.[NO3 -] and [NO2 -] or [DON] and [PON]) to be solved for by the linear solver.Vector b (400,320 x 1) contained the rates that were independent of the state variables, such as N2 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 25 In MATLAB, we used METIS ordering, which is part of SuiteSparse (http://faculty.cse.tamu.edu/davis/suitesparse.html) 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 5
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.Others, such as the rate constants, were objects of our investigation and were optimized against available observations.For our optimization, we compared model output using different parameter values to a database of NO3 -and NO2 -concentrations and isotopes.The database was originally compiled by Rafter et al. (in prep.) and has been expanded to include some additional unpublished data (Table S1).All of the database observations were 10 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.

15
The optimization procedure used the MATLAB function fminunc to obtain values for the non-fixed parameters that minimized a cost function.In each iteration of the optimization, the model system was solved by running the 14 N-ON model, 15 N-ON model, 14 N-DIN model, and 15  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: 25 The w terms are weighting terms intended to scale the contributions of the four observed parameters in order to equalize their 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.We chose to iterate over these parameters rather than optimize them since there is a large range of estimates from the literature as to what these parameters might be (Table 3).
We assigned different possible values for each of these parameters (Table 3), resulting in 12 possible combinations.The 5 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 NO2 -distribution and its isotopic composition are the most unique features of this model in comparison to 10 many other existing global models of the marine N cycle.As such, NO2 -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 NO2 -accumulating at 200m in the major ODZs of the Eastern Tropical North Pacific (ETNP), Eastern Tropical South Pacific (ETSP), and the Arabian Sea (AS) (Figure 2), which is consistent with observations and expected based on the low O2 conditions found there.However, accumulation of NO2 -in the model ETSP was lower than expected.The model also 15 accumulated NO2 -in the Bay of Bengal, which is a low-O2 region off the east coast of India that does not generally accumulate NO2 -or support water column denitrification, but is thought to be near the "tipping point" for allowing N loss to occur (Bristow et al., 2016a).The Bay of Bengal will be discussed further in Section 4.2.
The model optimization described above yielded a set of isotope effects that best fit the global dataset of [NO3 -], [NO2 -], 20 d 15 NNO3 and d 15 NNO2.The best fit was achieved for isotope effects of 13‰ for NO3 -reduction (εNAR), 0‰ for NO2 -reduction (εNIR), and -13‰ for NO2 -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 (Figure S3), indicating that we did not overfit the training data.25 In the test set, there were some low [O2] points where our model [NO3 -] exceeded observations (Figure 3a, filled black circles); these are largely within the Eastern Tropical South Pacific (ETSP), where insufficient NO3 -reduction occurred in the model.
In contrast, the Arabian Sea tended to show slightly lower modeled [NO3 -] than expected.In addition, the [NO2 -] accumulation (Figure 3b) and d 15 NNO3 signals (Figure 3c) in the ETSP were generally too low compared with observations.These signals 30 are likely tied to insufficient NO3 -reduction in the ETSP, although excessive consumption of NO2 -may also play a role.Another consideration is that there may be a mismatch in resolution between the model and the space and time scales needed to resolve the high NO2 -accumulations observed sporadically (Anderson et al., 1982;Codispoti et al., 1985;1986).Overall, the fits of δ 15 NNO3 were fairly good (RMSE = 2.4‰), though there were some points above δ 15 NNO3 = 10‰ where the modeled δ 15 NNO3 exceeded the observed δ 15 NNO3 and others where modeled δ 15 NNO3 was lower than observations (Figure 3c).Many of the points with overestimated δ 15 NNO3-were located within the Arabian Sea ODZ, where there may be too much NAR 5 occurring, leading to excess enrichment in 15 N-NO3 -.The fits of δ 15 NNO2 were also fairly good (RMSE = 8.6‰), though the modeled δ 15 NNO2 was generally not low enough (Figure 3d), indicating an underestimated sink of 'heavy' NO2 -.

Oxygen deficient zone model-data comparison using station profiles
To further investigate the distribution of model N species within the three main ODZs, we selected representative grid boxes within each ODZ that contained observations to directly compare with model results in station profiles.Overall, the modeled 10 NO3 -and NO2 -concentration and isotope profiles in the AS and ETNP were consistent with the observations, with [NO3 -] slightly under estimated in the Arabian Sea ODZ and overestimated in the ETSP (Figure 4).As [O2] goes to zero, the O2intolerant processes NAR, NIR, and AMX are released from inhibition.These processes result in a decrease in [NO3 -] (via NAR) which corresponds to an increase in δ 15 NNO3, since NAR has a normal isotope effect.NO2 -also starts to accumulate in the SNM as a result of NAR.δ 15 NNO2 is lower than δ 15 NNO3 since light NO2 -is preferentially created via NAR, and this 15 discrepancy 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 [NO3 -] depletion and [NO2 -] accumulation in the model were lower than observed.
In order the gauge the model results for N loss, we also calculated N*, a measure of the availability of DIN relative to PO4 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.30 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  (Martin and Casciotti, 2017).Although, in the Arabian Sea, modeled δ 15 NNO3 were too high in the ODZ, likely in part due to overpredicted rates of NAR, which also resulted in lower modeled [NO3 -] (Figure 4). 5

Model-data comparison in GEOTRACES sections
We also investigated the agreement between global best fit modeled [NO3 -] and d 15 NNO3 with data from two GEOTRACES cruise sections: GP16 in the South Pacific, and GA03 in the North Atlantic.For GP16, we see that [NO3 -] matches fairly well in the surface waters, but diverged below 500m as well as at the eastern edge of the ETSP ODZ (Figure 5).Although the patterns are generally correct, insufficient NO3 -is being accumulated in the deep waters of the model Pacific, which could be 10 due to an underestimate of preformed NO3 -(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 [NO3 -] are 5-10 µM lower than observations (Figure S4), which could be enough to explain the lower than expected  In the GP16 section, we also see that the insufficient depletion of NO3 -in the ETSP ODZ (Figure 5b) extends beyond the single grid box highlighted earlier (Figure 4).Along with this, there is also minimal increase of δ 15 NNO3-in the ETSP ODZ and the upper thermocline in the eastern part of the section, consistent with an underestimate of NO3 -reduction.Surface δ 15 NNO3 values were also not as high in the model as in the observations (Figure 5), which could result from insufficient NO3 -assimilation or too low supplied δ 15 NNO3 (Peters et al., 2018a).However, we do see a similar depth range for high surface δ 15 NNO3 and a local 20 δ 15 NNO3 minimum between the surface and ODZ propagating westward in both the model and observations, indicating that the physical processes affecting δ 15 NNO3 are represented by the model.Additionally, the model shows slightly elevated δ 15 NNO3 in the thermocline depths (200-500m) west of the ODZ, which is consistent with the observations (Figure 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 In the north Atlantic along GEOTRACES section GA03, we see good agreement between the observed and modeled [NO3 -] (Figure 6).There is generally low surface [NO3 -] with a distinct area of high [NO3 -] propagating from near the African coast.
Deep water (> 2000 m) [NO3 -] is lower than we see in the Pacific section, and the model matches well with the Atlantic observations.Again, there is not quite enough NO3 -present in Southern Ocean-sourced intermediate waters (500-1500 m; 5 Marconi et al., 2015).Modeled δ 15 NNO3 values at first glance appear higher than observed values at the surface (Figure 6).However, many of the surface [NO3 -] were below the operating limit for d 15 NNO3 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 NNO3 values in upper thermocline waters in both the model and observations, likely corresponding to low δ 15 N contributions from N2 fixation that is remineralized at depth and accumulated in North Atlantic Central Water (Marconi et al., 2015; Knapp 10 et al., 2008).The model input includes high rates of N2 fixation in the North Atlantic that are consistent with this observation (Martin et al., in prep).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-50% of N input along the cruise track.

Assumption checks
As previously mentioned (Section 2.3), ON and DIN were modeled separately in order to introduce dependence on both ON and substrate availability for the heterotrophic processes NAR and NIR.This separate modeling requires several assumptions to be made regarding the processes that impact both ON and DIN, namely assimilation and remineralization.

20
The first assumption is that the rates of N assimilation are equal between the ON and DIN model runs.The ON model run used WOA surface [NO3 -] to estimate the contribution of DIN assimilation to the production of ON, whereas the DIN model run used modeled [NO3 -] and [NO2 -] to estimate DIN removal via assimilation.Though these two methods utilized the same rate constants for assimilation, differences in [DIN] 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 ON produced in the organic 25 N model (Figure S5).This could be due in part to assimilation of NO2 -in the top two boxes, since NO2 -assimilation is unaccounted for in the ON model.This is largely an issue in the oligotrophic gyres, where surface [NO3 -] is very low and NO2 - accumulates to low but non-zero concentrations (Figure 2).Assimilation of NO2 -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.In other regions, modeled surface [NO3 -] may be higher than the WOA surface [NO3 -] data that are supplied to the 30 ON 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 ON production rates do tend to have modeled [NO3 -] that was higher than observed  S4).We conclude that though the assimilation rates are not identical in the ON 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 5 split model.

Model dependency on input O2
The modeled concentration and isotope profiles for the ETSP, unlike in the AS and ETNP, reflected a significant underestimation of water column denitrification in the best-fit model.In ETSP [NO3 -] measurements, there is a clear deficit in [NO3 -], coincident with the SNM and N* minimum (Figure 4).In our modeled profiles, this NO3 -deficit is missing, and 10 although a SNM is present, its magnitude is lower than observed (Figure 4).The model also does not capture the negative N* excursion (Figure 4), which we think reflects a model underestimation of NO3 -and NO2 -reduction in the ETSP.The cause of this missing denitrification is likely to be poor representation of the ETSP O2 conditions in the model grid space.Since our model grid is fairly coarse (2˚x2˚), only a few boxes within the ETSP had averaged [O2] 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 15 as in the AS and ETNP (Figure S6); therefore, this region in particular was less compatible with the model grid.In order to test whether the parameterization of O2 dependence was the cause of the low N loss, we ran the model using the globally optimized parameters (Table 3) but with higher O2 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 20 stimulated the development of a NO3 -deficit, larger secondary NO2 -maximum, and N* minimum within the ODZ (Figure 7).
As previously mentioned in Section 3.1, modeled [NO2 -] in the Bay of Bengal is higher than observations.The accumulation of NO2 -here in the model is likely due to O2 concentrations falling below the set threshold for NAR but above the threshold for NIR.Although AMX and NXR occur there, the modeled rates are rather slow, which leads to more NO2 -being produced 25 than consumed.This is in contrast to observations that NO2 -production via NAR is tightly matched with NO2 -consumption via NXR, which limits NO2 -accumulation and N loss in the Bay of Bengal (Bristow et al., 2016a).The fact that the model over predicts NAR in the Arabian Sea may also be connected with over-prediction of NAR in the Bay of Bengal, leading to accumulation of NO2 -despite realistic NO2 -consumption rates.Further work on oxygen sensitivities of N cycle processes will be addressed in a companion study (Martin et  ], [NO2 -], δ 15 NNO3, and δ 15 NNO2.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 5 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 [NO3 -] and [NO2 -] distributions that were required to achieve this roughly balanced budget are well within 10 the range of observed values.Some interesting take-home messages from this work are 1) a relatively low isotope effect for NO3 -reduction (εNAR = 13‰) gives a good fit to δ 15 NNO3 data, similar to that concluded in some recent studies (Marconi et al., 2017;Casciotti et al., 2013), 2) low O2 half-saturation constants for NO2 -oxidation allowing NXR to occur in parallel with NAR, NIR, and AMX were needed to achieve the correct distributions of NO3 -, NO2 -, and their isotopes, in the oceans water column ODZs.15 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 O2 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 20 represent some important temporal dynamics.While we attempted to account for seasonal variation in the strength of the ODZs, we did not simulate seasonal variations in NPP and the strength of the biological pump.
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.As discussed in Section 2.6, the rate parameters were optimized using a variety of different 25 combinations of isotope effects for NAR, NIR, and NXR to explore model uncertainty.As presented in Section 3.1, the larger magnitude isotope effect for NXR best fit the ETNP and Arabian Sea ODZs, where most of the ODZ volume resided.The larger magnitude isotope effect also resulted in optimized rate constants for NAR, NIR, and NXR that were lower than the global best fit rate constants.The lower rate constants and larger isotope effects resulted in better fits to observations of δ 15 N and DIN concentrations.This reinforces the importance of obtaining realistic isotope effect estimates for each process that are 30 relevant on an environmental scale.   kPON (yr -1 ) 3.9 3.9 0 3.9 14 kDON (yr N data) would make model validation difficult.A schematic diagram of the N cycle in the model is shown in Figure 1.Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.
second.When [DON] is found in the [DIN] governing equations, that [DON] value has already been determined for each grid Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.box from the ON model.When [NO3 -] is found in the DON governing equations, it is drawn from World Ocean Atlas (WOA13)

20 N2
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 N2 fixation in the model of DeVries et al. (2013), with partial inhibition by NO3 -(Holl and Montoya, 2005) and dependence on iron (Fe) and phosphate (PO4

5 Fe
, 2005) divided by the depth of the top model box.NO3,obs is the 2013 World Ocean Atlas (WOA) annually averaged surface NO3 -interpolated to the model grid (Garcia et al., 2014).λ is an inhibition constant for N2 fixation in the presence of NO3 -.Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.The temperature (T) terms scale the rate of N2 fixation based on the observed temperature (Tobs), maximum observed sea surface temperature (Tmax), and the minimum preferred growth temperature for Trichodesmium (T0; Capone et al., 2005).The temperature data were taken from 2013 WOA annually averaged temperature interpolated to the model grid (Locarnini et al., 2013). is the modeled deposition of soluble Fe interpolated to the model grid (mmol Fe/m 2 /yr; Chien et al., 2016) divided by the flux attenuation is somewhat longer in ODZs compared to oxygenated portions of the water column varies regionally (Berelson Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.et al., 2002; Buesseler et al., 2008; Buesseler and Boyd, 2009), our model uses a spatially-invariant 14 ksol.This is a refinement that could be introduced in future model versions.

(
[O2] -[NO3 -]).In order for sedimentary denitrification to be properly implemented in our linear model, we broke the original Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.non-linear relationship into three roughly linear segments.We obtained three linear relationships between ([O2] -[NO3 -]) and sedimentary denitrification rate, each applicable across a given range of ([O2] -[NO3 -]) values (Figure S1).Due to the nature of our linear model, we needed to express our cutoff points between the segments in terms of O2 rather than ([O2] -[NO3 -]).Therefore, a linear relationship between O2 and ([O2] -[NO3 -]) was determined using the 2013 WOA annually averaged data (Garcia et al., 2014; Figure S2).The cutoff points were determined to be 75 µM O2 and 175 µM O2.The linear relationships 5 were then rearranged in order to estimate sedimentary denitrification rate as a function of RRPOC, [O2], and [NO3 -].These equations were then further broken down into a component that is dependent on [NO3 -] and a component that is dependent on [O2].
N/ 14 N from 25 measured d 15 N, and assumed that [ 14 N] ~ [ 14 N] + [ 15 N], the measured concentration of the modeled N species.
produce and consume NO3 -that are also dependent on [NO3 -].The top right corner of matrix A would contain rate constants for processes that produce and consume NO3 -but are dependent on [NO2 -].The bottom left corner of matrix A would contain rate constants for processes that produce and consume NO2 -but are dependent on [NO3 -].The bottom right corner of matrix A would contain rate constants that produce and consume NO2 -and are also dependent on [NO2 -].The top half of vector x would be [NO3 -] for each model box, and the bottom half of vector x would be [NO2 -] for each model box.The top half of vector b 30 would be independent processes that produce or consume [NO3 -], and the bottom half of vector b would be independent processes that produce or consume [NO2 -].Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.
N-DIN model.The modeled output [NO3 -], [NO2 -], δ 15 NNO3, and δ 15 NNO2 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 NO3 -values were useful in constraining the parameters that control PON solubilization and DON 20 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 30 data set (e.g. the number of [NO3 -] data points = nNO3).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 [NO3 -] data points within the training set).Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.
3- 20 compared to Redfield ratio stoichiometry (N* = [NO3 -] + [NO2 -] -16 * [PO4 3-]; Deutsch et al., 2001).Negative N* values are associated with N loss due to AMX or NIR, while positive N* values are associated with input of new N through N2 fixation (Gruber and Sarmiento, 1997).Although we did not model PO4 3-, we used the modeled [NO3 -] and [NO2 -] together with WOA PO4 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 25 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.

d 15 NNO3
across the basin.Peters et al. (2018a) and Rafter et al. (2013) also postulated that these elevated δ 15 NNO3 values were 25 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 other N isotope studies (Rafter et al., 2013; Peters et al., 2018).The model also shows slightly elevated δ 15 NNO3 in the intermediate depths (500-1500m), which is consistent with observations, again reflecting remineralization of PON with d 15 N greater than mean ocean d 15 NNO3.Overall, the patterns of δ 15 NNO3 are correct but the magnitudes of isotopic variation are muted, largely due to the lack 30 of N loss in the ODZ and modeled surface δ 15 NNO3 values that are lower than observations.Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.
[NO3 -] (Figure Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.S5).Likewise, points with relatively lower DIN assimilation had modeled [NO3 -] less than observed [NO3 -].However, the majority of DIN assimilation estimates were within 10 µM/yr of the ON production estimates, with an overall DIN assimilation rate of 9235 µM/yr.We also found that the WOA surface NO3 -values are fairly well represented by our modeled surface NO3 - (Figure al., in prep).30 Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License. 5 Conclusions A global inverse ocean model modified to include [NO2 -], d 15 NNO3, and δ 15 NNO2 as state variables, as well as adding the processes required to describe the cycling of NO2 -in the global ocean resulted in a globally representative distribution of [NO3 - Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-397Manuscript under review for journal Biogeosciences Discussion started: 17 September 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 1 : 5 (
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 (NO3 -) and nitrite (NO2 -).N source processes are nitrogen (N2) fixation and atmospheric deposition.N sink processes are sedimentary denitrification, NO2 -reduction 5

Figure 4 :
Figure 4: Depth profiles comparing model results with binned and averaged database observations from a model water column.Results are shown for the three main oxygen deficient zones (ODZs): the Arabian Sea (a-e), the Eastern Tropical North Pacific (ETNP; f-j), and the Eastern Tropical South Pacific (ETSP; k-o).Average modeled nitrate concentration ([NO3]), nitrite concentration ([NO2 -]), and N* are shown in black.Gray error lines around the black line show the 2σ spread from the average from 5

Figure 5 :
Figure 5: Section profiles of NO3 -concentrations and isotopes over the GP16 cruise track (panel (a) inset) in the South Pacific.Profiles are presented from east (right) to west (left).Comparison of (a) observed [NO3 -] to (b) modeled [NO3 -] is presented over the full depth range (0-6000m).Comparison of (c) observed δ 15 NNO3-to (d) modeled δ 15 NNO3-is presented over a shortened depth range (0-1000m) to better assess surface and ODZ values.GEOTRACES data are from Peters et al. (2018a) and available from BCO-5

Figure 6 :
Figure 6: Section profiles of NO3 -concentrations and isotopes over the GA03 cruise track (panel (a) inset) in the North Atlantic.Profiles are presented from east (right) to west (left) from 0-6000km section distance, and then from south to north.Comparison of (a) observed [NO3 -] to (b) modeled [NO3 -] is presented over the full depth range (0-6000m).Comparison of (c) observed δ 15 NNO3-to (d) modeled δ 15 NNO3-is presented over a shortened depth range (0-1000m) to better assess surface and the low δ 15 NNO3-contribution 5

Figure 7 : 5 [
Figure 7: Plot of DIN concentrations and N* from the ETSP ODZ comparing modified O2 thresholds for N loss.In the original optimized version of the model, there is insufficient N loss and NO2 -accumulation in the ETSP.To demonstrate that this issue may be caused by issues with the gridded averages of O2 and model grid size in the ETSP, we raised the O2 thresholds for N loss-related processes (NAR, NIR, and AMX) to 15 µM.This effectively lowers the observed [O2] in order to stimulate N loss.The resulting (a) 5