Isotopic constraints on the pre-industrial oceanic nitrogen budget

The size of the bioavailable (i.e., “fixed”) nitrogen inventory in the ocean influences global marine productivity and the biological carbon pump. Despite its importance, the pre-industrial rates for the major source and sink terms of the oceanic fixed nitrogen budget, N 2 fixation and denitrification, respectively, are not well known. These processes leave distinguishable imprints on the ratio of stable nitrogen isotopes, δ15N, which can therefore help to infer their patterns and rates. Here we use δ15N observations from the water column and a new database of seafloor measurements to constrain rates of N 2 fixation and denitrification predicted by a global three-dimensional Model of Ocean Biogeochemistry and Isotopes (MOBI). Sensitivity experiments were performed to quantify uncertainties associated with the isotope effect of denitrification in the water column and sediments. They show that the level of nitrate utilization in suboxic zones, that is the balance between nitrate consumption by denitrification and nitrate replenishment by circulation and mixing (dilution effect), significantly affects the isotope effect of water column denitrification and thus global mean δNO−3 . Experiments with lower levels of nitrate utilization within the suboxic zone (i.e., higher residual water column nitrate concentrations, ranging from 20 to 32 μM) require higher ratios of benthic to water column denitrification, BD : WCD= 0.75–1.4, to satisfy the global mean NO −3 and δNO−3 constraints in the modern ocean. This suggests that nitrate utilization in suboxic zones plays an important role in global nitrogen isotope cycling. Increasing the net fractionation factorεBD for benthic denitrification ( εBD = 0–4 ‰) requires even higher ratios, BD : WCD = 1.4–3.5. The model experiments that best reproduce observed seafloor δ15N support the middle to high-end estimates for the net fractionation factor of benthic denitrification ( εBD = 2–4 ‰). Assuming a balanced fixed nitrogen budget, we estimate that preindustrial rates of N2 fixation, water column denitrification, and benthic denitrification were between 195–350 (225), 65– 80 (76), and 130–270 (149) Tg N yr −1, respectively, with our best model estimate ( εBD = 2 ‰) in parentheses. Although uncertainties still exist, these results suggest that marine N 2 fixation is occurring at much greater rates than previously estimated and the residence time for oceanic fixed nitrogen is between∼ 1500 and 3000 yr.


Introduction
Biotically available "fixed" nitrogen (fixed-N) is one of the major nutrients limiting phytoplankton growth in the ocean.Its generally low abundance in sunlit surface waters limits the primary production that forms the base of ocean ecosystems and provides energy for more complex, higher-level organisms (e.g., marine animals).Thereby, fixed-N also limits the biological sequestration of atmospheric carbon dioxide (CO 2 ) into biomass, part of which subsequently sinks towards the deep ocean before being respired back into CO 2 .This so-called "biological carbon pump" affects the partitioning of CO 2 among the ocean and atmosphere.It has been suggested that large changes in the oceanic fixed-N inventory can modulate the strength of the biological carbon pump and thereby influence atmospheric CO 2 over glacial/interglacial timescales (McElroy, 1983;Falkowski, 1997).Since processes that control the size of the fixed-N inventory are sensitive to climate (Gruber, 2004;Galbraith and Kienast, 2013), they may have an important feedback on atmospheric CO 2 concentrations in past and future climates.
Published by Copernicus Publications on behalf of the European Geosciences Union.
N 2 fixation (NFix), the conversion of N 2 gas into fixed-N by specialized microorganisms (diazotrophs), provides the ocean with most of its fixed-N.Other contributions to the fixed-N pool are from river input and atmospheric N deposition, which are estimated to be approximately an order of magnitude lower than NFix in pre-industrial times (Galloway et al., 2004;Codispoti, 2007;Duce et al., 2008;Gruber, 2008).However, industrial N emissions and fertilizer production through the Haber-Bosch process that eventually cycles fixed-N into the atmosphere and rivers have been steadily increasing in recent decades and are estimated to become comparable to "natural" NFix in following decades (Galloway et al., 2004).
Denitrification and anammox (anaerobic ammonium oxidation) are the most important fixed-N removal processes in the ocean.They convert fixed-N into N 2 gas under suboxic conditions (O 2 < ∼ 10 µM) in the water column and seafloor sediments."Canonical" denitrification occurs when heterotrophic bacteria replace O 2 consumption with the reduction of nitrate (NO − 3 ) to dinitrogen gas (N 2 ) as the electron acceptor during respiration, once O 2 is no longer available in sufficient quantity.Anammox, on the other hand, is a chemoautotrophic process that converts nitrite (NO − 2 ) and ammonium (NH + 4 ) into N 2 gas (Thamdrup and Dalsgaard, 2002;Kuypers et al., 2003).Since both denitrification and anammox processes occur in oxygen minimum zones (OMZs), the relative importance of each process is uncertain.Recent studies have found that water column denitrification (WCD) dominates N-loss in the Arabian Sea (Ward et al., 2009;Bulow et al., 2010), while anammox is more important in the eastern tropical South Pacific (ETSP) (Lam et al., 2009).However, Lam et al. (2009) estimate that nitrate reduction, the first step of canonical denitrification, provides at least two-thirds of the nitrite that anammox consumes, suggesting that canonical denitrification may be the most important driver of total N-loss in OMZs.Therefore, we refer to denitrification as the major fixed-N loss process in this paper.
While the major source and sink processes of the fixed-N budget, N 2 fixation and denitrification, respectively, are relatively well known, estimating their global rates remains a challenge.Model estimates for N 2 fixation have a large range between ∼ 100 and 300 Tg N yr −1 and predict different spatial patterns (Gruber and Sarmiento, 1997;Brandes and Devol, 2002;Deutsch et al., 2007;Monteiro et al., 2011;Eugster and Gruber, 2012), whereas estimates based on extrapolation of measured in situ rates agree with the low-end model estimates near 100 Tg N yr −1 (Karl et al., 2002a;Mahaffey et al., 2005).However, methods historically used to measure N 2 fixation rates have been found to underestimate N 2 fixation by as much as a factor of 2 (Mohr et al., 2010;Großkopf et al., 2012).Since new important N 2 -fixing species are also still being discovered (Montoya et al., 2004;Foster et al., 2011;Zehr, 2011), more N 2 fixation could be taking place in the ocean than previously thought.
Global estimates of denitrification vary considerably as well.In the water column, studies suggest global rates between 50 and 150 Tg N yr −1 and in the sediments generally between 100 and 300 Tg N yr −1 (Middelburg et al., 1996;Galloway et al., 2004;Gruber, 2004;Codispoti, 2007;Bohlen et al., 2012;DeVries et al., 2013;Eugster and Gruber, 2012).Since the high-end estimates for denitrification are substantially larger than the high-end estimates for N 2 fixation, it has led to the debate whether the ocean could be rapidly losing as much as 400 Tg N yr −1 (Codispoti et al., 2001;Brandes and Devol, 2002;Seitzinger et al., 2006;Codispoti, 2007), about 0.07 % of the total nitrate inventory per year, or whether the nitrogen budget is more balanced (Gruber and Sarmiento, 1997;Gruber, 2004;Altabet, 2007;Bianchi et al., 2012).
The observed global mean δ 15 NO − 3 provides a constraint on the fixed-N budget (Brandes and Devol, 2002) since it is determined by the relative rates of NFix, WCD and BD and their isotope effects.The input of 15 N-depleted nitrogen from NFix must compensate for the preferential removal of 15 N-depleted nitrate by denitrification.Modeling studies using this approach have produced differing results.For example, Brandes and Devol (2002) used a one-box model with fractionation factors of ε BD = 1.5 ‰ and ε WCD = 25 ‰.They estimated a ratio BD : WCD = 3.7 : 1 to match global mean δ 15 NO − 3 , supporting the high-end estimates for BD.Other studies have used more complex box models (Deutsch et al., 2004;Eugster and Gruber, 2012).Assuming ε BD = 0 ‰ and ε WCD = 25 ‰, these studies estimate a lower BD : WCD ratio of 1.8-2.7 compared to Brandes and Devol (2002) by accounting for mixing of locally reduced nitrate concentrations in suboxic zones with surrounding oxic waters.Since the mixed δ 15 NO − 3 value will be weighted towards the higher nitrate concentration of the oxic waters, the local utilization of nitrate reduces the influence of the 15 N-enriched nitrate from WCD on global mean δ 15 NO − 3 , a mechanism referred to as the "dilution effect" (Deutsch et al., 2004).Altabet (2007) used another one-box model with parameterized isotope effects of local nitrate utilization and dilution during WCD.He argues that these effects would further reduce the isotope effect of WCD by ∼ 13 ‰ compared to the inherent microbial process near 25 ‰.Applying this reduced isotope effect for water column denitrification (∼ 12 ‰), and 0 ‰ for BD, respectively, he estimated a BD : WCD ratio of ∼ 1 : 1.These box model studies highlight how sensitive the results can be to different assumptions made for the isotope effects of denitrification.
In this study, we go beyond earlier box model approaches and employ a global coupled three-dimensional circulationbiogeochemistry-isotope model to investigate to what extent rates of N 2 fixation and denitrification can be constrained by δ 15 N observations in the water column and seafloor sediments.In particular, we will investigate the uncertainties associated with (i) the effects of nitrate utilization and dilution on the isotope effect of WCD and (ii) the net fractionation factor associated with BD.In addition to water column δ 15 NO − 3 observations, a new global seafloor sediment δ 15 N database (Tesdal et al., 2013) is used to evaluate the model experiments.The rates of N 2 fixation and denitrification implicit in the model simulations that most closely simulate observed δ 15 NO − 3 and seafloor δ 15 N then provide our best estimate of Nfix and denitrification in the real ocean.

Model description
The global coupled Model of Ocean Biogeochemistry and Isotopes (MOBI) is based on the version of Somes et al. (2010b).A technical description of the model is located in Appendix B and a brief overview is provided below.

Biogeochemical-ecosystem model
MOBI is an improved version of the model used in Somes et al. (2010b) (see Fig. 1).The organic compartments include two classes of phytoplankton, N 2 -fixing diazotrophs (P Diaz ) and other nitrate assimilating phytoplankton (P O ), as well as zooplankton (Z) and sinking detritus (D).The inorganic variables include dissolved oxygen (O 2 ) and two nutrients, nitrate (NO − 3 ) and phosphate (PO 3− 4 ), both of which are consumed by phytoplankton.
This model version deviates from that of Somes et al. (2010b) by including varying elemental stoichiometry.While general phytoplankton nitrogen to phosphorus ratio (N : P) remains at the canonical Redfield ratio (N : P = 16), diazotroph N : P is set to 40, which is in better agreement with most observations (N : P = 20-50+; e.g., Letelier andKarl, 1998 andSañudo-Wilhelmy et al., 2004), as well as an optimality-based growth model (Klausmeier et al., 2004).
This allows diazotrophs to more efficiently fix N 2 into the ocean when they are P-limited, but does not significantly change the pattern of N 2 fixation.Since zooplankton are capable of maintaining their own stoichiometry (Sterner and Elser, 2002), we set zooplankton N : P to 16:1 and assume that they excrete excess N when grazing diazotrophs.Detrital N and P are explicitly calculated as separate prognostic variables.The carbon to nitrogen (C : N) ratio for all compartments is held constant at C : N = 6.625.

N 2 fixation
Diazotrophs grow according to the same principles as the general phytoplankton class in MOBI, but we account for some of their different characteristics as follows.N 2 fixation is energetically more costly than assimilating fixed-N.Diazotrophs must break down the strong triple-N bond and undergo extra respiration to keep the N 2 -fixing compartment anoxic since O 2 inhibits the expression of the N 2 -fixing gene (nifH).Thus, the growth rate of diazotrophs is reduced compared to ordinary phytoplankton, which MOBI considers by using a handicap factor (c Diaz = 0.13).Diazotrophs are allowed to grow at low rates in cool waters according to Eq. (B4), consistent with culture experiments (Pandey et al., 2004;Le Quéré et al., 2005), whereas in previous model versions (Schmittner et al., 2008;Somes et al., 2010b) their growth rates were zero below 15 • C.
Diazotrophs are not limited by nitrate and can outcompete general phytoplankton in surface waters that are depleted in fixed-N, but still contain sufficient phosphate and iron.They will consume nitrate if available, consistent with culture experiments (Mulholland et al., 2001;Holl and Montoya, 2005), which is another factor that inhibits N 2 fixation in high-nutrient, low-chlorophyll regions in the model.The prey-capture rate for zooplankton feeding on diazotrophs is reduced relative to that for the general phytoplankton class in order to account for observations (Table 1) (Mulholland, 2007), in contrast to Somes et al. (2010b), who used equal rates.
Due to uncertainties in the iron cycle (Aumont et al., 2003;Moore et al., 2004;Galbraith et al., 2010), iron is currently not included as a prognostic tracer.Instead, iron limitation is considered using a monthly mask based on modeled aeolian dust deposition estimates (Mahowald et al., 2005a(Mahowald et al., , 2006(Mahowald et al., , 2009)).This iron limitation mask is multiplied by diazotrophs' maximum growth rate to account for iron limitation (Somes et al., 2010a).Note that this diazotroph iron limitation mask does not account for upwelled iron.Since upwelled iron will be accompanied by other macronutrients, we assume this source of iron is consumed by other phytoplankton.This simplification may lead to an overestimation of iron limitation of diazotrophs in and around some upwelling regions because complete utilization of upwelled iron by other phytoplankton will likely not occur in all of these areas (e.g., Equatorial Atlantic; Subramaniam et al., 2013).We scaled our iron limitation mask of diazotrophs to best reproduce large-scale meridional δ 15 N patterns across the Pacific and Atlantic oceans.

Water column denitrification
WCD occurs when organic matter is respired under suboxic conditions.It is parameterized in MOBI according to Eq. (B14), which determines the relative amount of nitrate consumption that takes place during respiration at low ambient oxygen concentrations.We use a threshold of 5 µM O 2 that sets the oxygen level at which respiration by denitrification overtakes aerobic respiration.This parameterization was designed to best capture the decreasing NO − 3 : PO 3− 4 ratios in suboxic zones.Anammox is also removing fixed-N in these areas of low oxygen and high organic matter recycling.Although the exact partitioning between WCD and anammox is not well known, anammox likely depends on the first step of denitrification (NO − 3 → NO − 2 ) to supply sufficient nitrite that typically exists in very low concentrations (Lam et al., 2009).Since MOBI does not differentiate between different species of dissolved inorganic nitrogen, this WCD parameterization is designed to implicitly capture total fixed-N loss from both canonical WCD and anammox.

Benthic denitrification
BD is included using a new empirical function deduced from benthic flux measurements (Bohlen et al., 2012) .This function estimates denitrification in the sediments based on organic carbon rain rate into the sediments, and bottom water O 2 and NO − 3 concentrations (Eq.B15).It provides an efficient alternative to coupling a full sediment model that would significantly reduce simulation speed.Note that all organic matter instantaneously remineralizes in the bottom water when it reaches the seafloor sediment interface.Nitrate is then removed from the bottom water according to this benthic denitrification function.While the organic carbon rain rate predicts benthic denitrification rates to firstorder, Bohlen et al. ( 2012) found a strong non-linear relationship to the empirical parameter O * 2 = O 2 − NO − 3 in the bottom water overlying the sediments.BD rates are significantly higher as O * 2 decreases (i.e., when O 2 is low and NO − 3 is high) for similar organic carbon rain rates.We use high resolution (1/5 • ) bathymetry to account for shallow continental shelves and other topographical features that are not fully resolved in MOBI's coarse-resolution grid in order to calculate BD (Eq.B15).This scheme considers particulate organic matter sinking and remineralization, but it does not influence the physical circulation.Since details of the circulation that are not resolved in our global model are responsible for nutrient fluxes in these areas (Fennel et al., 2006), our model will likely underestimate productivity, rain rate of carbon into the seafloor, and thus BD there if the BD rate is not further adjusted.This is one of the reasons why we provide sensitivity experiments that linearly increase BD rates by different factors (α BD = 1.35-3.2;Tables 1, 2).Note that water in contact with this subgrid-scale shelf scheme can still be influenced from water below through physical mixing and circulation that can result in too much coupling with the ocean interior.

Nitrogen isotope model
The nitrogen isotope model simulates the two stable nitrogen isotopes, 14 N and 15 N, in all N species included in MOBI.
Fractionation results in the preferential incorporation of the more reactive, thermodynamically preferred 14 N isotope into the product of each reaction by a process-specific fractionation factor, α (Mariotti et al., 1981).We report these values in the more commonly used "δ" notation, where the fractionation factor takes the form ε = (1 − α) × 1000.The processes in the model that fractionate nitrogen isotopes are phytoplankton nitrate assimilation (ε assim = 5 ‰), N 2 fixation (ε NFix = 1.5 ‰), zooplankton excretion (ε excr = 6 ‰), WCD (ε WCD = 25 ‰), and BD (ε BD = 0, 2, 4 ‰).For example, diazotroph biomass becomes 1.5 ‰ depleted in δ 15 N relative to the source (atmospheric δ 15 N 2 = 0 ‰) giving diazotroph biomass a δ 15 N signature of −1.5 ‰ when they fix atmospheric N 2 for growth.We refer to the fractionation factor as the ε value chosen for the fractionation equation and the "isotope effect" as the overall effect the process has on the δ 15 N distribution in the ocean.The total isotope effect also includes effects from source values and processes that alter the impact of the net fractionation, such as nitrate utilization and dilution.NFix, for example, has a low fractionation factor, but a strong isotope effect by introducing very 15 N-depleted nitrogen (δ 15 N NFix = −1.5 ‰) into the oceanic fixed-N pool relative to deep ocean mean δ 15 NO − 3 near 5 ‰.In MOBI, fractionation factors are constant in space and time, and chosen to best represent estimates from field observations (Somes et al., 2010b).

Sensitivity experiments
The model experiments were initialized with World Ocean Atlas 2009 (WOA09) observations (temperature, salinity, oxygen, nitrate, and phosphate) (Antonov et al., 2010;Garcia et al., 2010a, b;Locarnini et al., 2010) and integrated for over 10 000 yr with pre-industrial boundary conditions as the nitrogen cycle approached equilibrium.While the scarcity of water column observations makes it difficult to estimate global mean δ 15 NO − 3 , vertical δ 15 NO − 3 profiles throughout the global ocean converge to ∼ 5 ‰ at depths below 2000 m (Somes et al., 2010b).In our experiments we manually adjust the global NFix and BD rates to match this deep ocean δ 15 NO − 3 value and observed global mean NO − 3 ∼ = 30.8µM by modifying uncertain parameters (Tables 1, 2).Adjusting α BD determines the total BD rate.For NFix, we adjust the preycapture rate, ω Diaz , and mortality rate, ν Diaz , of diazotrophs, which regulate their net growth.In the following sections, annual mean results from the final 500 yr of the integrations are reported.

Water column denitrification experiments
This set of experiments was designed to test the importance of the isotope effect of WCD on the global ocean mean δ 15 NO − 3 .We follow a method introduced by Moore and Doney ( 2007) that limits the rate of WCD (limWCD) at nitrate thresholds of 20, 26, and 32 µM (Table 2, Fig. 2).These values were chosen to produce WCD rates that lie within the range of modern estimates between 50 and 150 Tg N yr −1 .Note that recent studies report a new low-end estimate for water column denitrification rate (∼ 50 Tg N yr −1 ) (Eugster and Gruber, 2013;DeVries et al., 2013), which is not explicitly tested in our model experiments that cover the range 75-140 Tg N yr −1 .
Only the highest threshold applied (limWCD = 32) captures the average nitrate concentration in the suboxic zones suggested by WOA09 (∼ 32 µM), while the other experiments underestimate nitrate there.Models with higher thresholds also have lower levels of nitrate utilization in the suboxic zones that influences the isotope effect of WCD.In order to exclude impacts on the distribution of other biogeochemical tracers (oxygen, phosphate), remineralization of organic matter was left unchanged in these sensitivity experiments.Only nitrate consumption during WCD was turned off below the respective nitrate threshold.

Benthic denitrification fractionation experiments
Field studies estimate the fractionation factor of BD to be much smaller than that of water column denitrification, but to what extent still remains uncertain.The small increase in bottom water δ 15 NO − 3 overlying BD zones suggests the fractionation factor should be in the range 0-3 ‰ (Brandes and Devol, 2002;Lehmann et al., 2004)   recent studies have suggested a much higher net fractionation (Lehmann et al., 2007;Granger et al., 2011;Alkhatib et al., 2012).They suggest the 15 N-enriched ammonium measured in the overlying bottom water, presumably released from the sediments, was due to fractionation during the nitrificationdenitrification loop in the sediments.These studies indicate a net fractionation factor for BD in the range 4-8 ‰.
However, Lehmann et al. (2007) show that shallow regions have a higher net fractionation compared to deep ocean sites and the global average net fractionation is likely closer to 4 ‰.Thus we performed experiments with BD net fractionation at 0, 2, and 4 ‰.Since these experiments nearly cover the full range of previous estimates for BD rates between ∼ 100-300 Tg N yr −1 and BD's net fractionation factor, they likely test the full range of uncertainty of BD's effect on the global nitrogen isotope inventory.

N 2 fixation
Diazotrophs perform N 2 fixation primarily in the tropics/subtropics where sufficient atmospheric Fe deposition occurs (e.g., North Atlantic, western Pacific, and north Indian Ocean; Fig. 2).Besides temperature and Fe availability, another important factor that determines where model diazotrophs are most abundant includes competition for resources with other phytoplankton (i.e., N vs. P limitation).Denitrification consumes nitrate and increases Nlimitation "downstream" of denitrification zones.In MOBI, this creates ecological niches for diazotrophs because they are not limited by nitrate.This applies not only to the three WCD zones, but also to the BD zones that occur with greater global rates and across all ocean basins.In fact, it is the higher BD rates in the Atlantic Ocean that stimulate higher N 2 fixation rates there in experiment #5 (Fig. 2, right panel) compared to experiment #3 (Fig. 2, center panel).Nitrate depletion at the surface thus stimulates N 2 fixation as long as temperatures are high enough and sufficient Fe and P is available.The spatial pattern of N 2 fixation is similar in all experiments, but with differing rates depending on the total denitrification rate in each experiment.

Water column denitrification
In MOBI WCD occurs when organic matter is respired in suboxic zones of the eastern tropical North/South Pacific, the Bay of Bengal, and the eastern tropical Atlantic (Fig. 2).Note that in the real ocean, large rates of WCD have been observed in the Arabian Sea, which is not reproduced in the model.Instead, the model displaces the suboxic zone to the Bay of Bengal, which is very close to the suboxic threshold in nature.Similar discrepancies between simulated and observed regions of WCD in the Indian Ocean have also been found in other models (Moore and Doney, 2007;Gnanadesikan et al., 2012).Its causes are not fully understood and may include a misrepresentation of coastal currents or precipitation in the Indian Ocean.Although denitrification is not regularly observed in the Atlantic Ocean, anammox has been measured in the eastern tropical South Atlantic suggesting that N-loss events can occur there (Kuypers et al., 2005).Nevertheless, the model likely significantly overestimates water column Nloss in the Atlantic Ocean.
Uncertainties in historical O 2 measurements and interpolation methods to a global scale make it difficult to assess the global volume of suboxic waters in MOBI.The global suboxic volume (O 2 < 10 µM) in our experiments is approx-imately 3 times larger than suggested by uncorrected annual WOA09 observations (Garcia et al., 2010a).However, Bianchi et al. (2012) found that after correcting for biases in historical O 2 measurements in suboxic zones and applying improved mapping/interpolation methods, this observational dataset is likely underestimating the global volume of suboxic waters by a factor of 3, which would be in general agreement with our model experiments.
Despite the uncertainty in global suboxic volume, it remains clear that the suboxic zones in the Pacific and Atlantic are displaced over the productive equatorial regions, where the large amounts of remineralization that occurs in the water column would be expected to result in too much WCD.This model bias can be attributed to the underestimation of equatorial and coastal undercurrents that ventilate the OMZs.
To account for this model deficiency, we provide sensitivity experiments that limit WCD in the upper part of the OMZ at given NO − 3 thresholds to produce WCD rates that are within previous estimates (see Sect. 2.4.1) and produce more realistic δ 15 NO − 3 values and NO − 3 utilization levels in the suboxic zones (see Sect. 3.2.2).Note that the suboxic zones in MOBI also extend too deep in the water column (to ∼ 1500 m), whereas according to WOA09 observations they do not extend below ∼ 1000 m.Since remineralization rates are much lower at these depths, only 15 % of total WCD occurs here which suggests that our best model may still be overestimating WCD by ∼ 15 %.

Benthic denitrification
BD occurs where large amounts of particulate organic carbon (POC) sink into the seafloor sediments (Middelburg et al., 1996;Seitzinger et al., 2006;Fennel et al., 2009;Bohlen et al., 2012).The largest BD rates in the model occur on the continental shelves when primary production is high (e.g., Bering Sea Shelf, Fig. 2).Although rates are much lower in the deep seafloor due to less sinking POC, it contains a much larger area and still significantly contributes to the global BD rate in the model.The percentage of total BD (including α BD ) simulated on the continental shelves (0-200 m), slopes (200-2000 m), and deep seafloor (2000-6000 m) is approximately 40 %, 40 %, and 20 %, respectively, which does not vary by more than ±5 % for the different model experiments.BD occurs across all ocean seafloor basins and thus has a more broad global distribution compared to WCD, which is confined to three distinct tropical regions (Fig. 2).

Seafloor δ 15 N
We compare sinking δ 15 N in the model with a new global database of seafloor δ 15 N (Tesdal et al., 2013) (Fig. 3).It is composed of 2347 bulk δ 15 N measurements covering all ocean basins.This seafloor δ 15 N database provides a more  (Robinson et al., 2012).See text for additional details.complete view of the global nitrogen isotope distribution compared to available water column δ 15 NO − 3 observations, which are sparse in space and time (Somes et al., 2010b).Since seafloor δ 15 N measurements represent the accumulation of material spanning the last hundreds to thousands of years, they remove seasonal, interannual, or anthropogenic variability that can impact any single water column observation, making seafloor δ 15 N an ideal dataset to constrain the long-term average of the pre-industrial ocean.
Diagenetic alteration of δ 15 N occurs as sinking δ 15 N becomes buried in seafloor sediments, which can potentially bias the interpretation of bulk sediment δ 15 N. A recent analysis shows a clear relationship of diagenetic alteration with water depth (Robinson et al., 2012).Data at > 50 sites where sinking δ 15 N from traps were compared to seafloor δ 15 N show a ∼ 0.8-1 ‰ δ 15 N diagenetic increase per kilometer of water depth.In our model-data comparison, we accordingly adjust the simulated sinking δ 15 N of the particulate organic nitrogen (PON) reaching the seafloor by increasing it by 0.9 ‰ km −1 to account for this effect that is not incorporated into MOBI.This diagenetic enrichment increases the globally averaged δ 15 N-PON in the sediments by 3.36 ‰ in the model experiments performed here and makes global mean values consistent with the observations (Table 3).Note that this seafloor diagenesis adjustment is only applied for the model-data comparison and does not affect isotope mass conservation in the model.The data-masked global seafloor δ 15 N average varies somewhat depending on rates of N 2 fixation and denitrification in the different model experiments, but all model experiments show similar values across the  ), the normalized (by the STD of the observations) standard deviation (NSTD), and the normalized (by the STD of the observations) root mean squared error (NRMSE).Observational estimates, where available, are given in parentheses.Note "full mean" refers to the entire ocean basin, whereas "data-masked mean" only include grid points where observations exist.
Southern Ocean (Table 3), which is not significantly affected by denitrification and N 2 fixation.
MOBI reproduces the major trends of the seafloor δ 15 N dataset (Fig. 3).The isotope effects of phytoplankton NO − 3 assimilation, WCD, and N 2 fixation are mainly responsible for these large-scale patterns of δ 15 N in the model.Relatively 15 N-depleted sinking nitrogen hitting the seafloor (δ 15 N = 0-4 ‰) occur in high-nutrient, low-chlorophyll regions of the Southern Ocean and eastern equatorial Pacific where NO − 3 utilization by phytoplankton is low.Here phytoplankton are able to preferentially incorporate 15 N-depleted nitrate into their biomass due to high availability of nitrate (ε assim = 5 ‰).The nitrate utilization isotope effect also produces more 15 N-enriched nitrate in surface waters as utilization increases.In the subtropical gyres where nitrate is nearly fully utilized, phytoplankton must consume this 15 Nenriched nitrate remaining in the surface water, causing much higher sinking δ 15 N values there (> 6 ‰), in the absence of N 2 fixation.All model experiments produce similar patterns and regional averages across the Southern Ocean where surface NO − 3 utilization dominates the δ 15 N trend (Table 3, Fig. 4).
Very high seafloor δ 15 N values (> 10 ‰) are observed near suboxic zones due to the large fractionation factor of WCD.Modeling the correct extent of the suboxic zones remains a challenge in global coarse-resolution models due to the limited spatial extent of suboxic zones.While the suboxic zones are generally simulated in the correct ocean basins (e.g., eastern tropical Pacific, northern Indian), they are all displaced in the model.Since WCD has a strong local effect on δ 15 N, this displacement causes rather poor model fits when comparing to the seafloor δ 15 N database (e.g., r < 0.6, Table 3).If all OMZ regions with less than 30 µM dissolved O 2 in the model and WOA09 are not included when calculating the global metrics, the correlation coefficient increases from 0.59 to 0.68 (in model experiment #4), an indication that the bias due to the displaced suboxic zones is one of the seafloor δ 15 N (Fig. 3), as well as water column observations (Somes et al., 2010a).Lowest values of seafloor δ 15 N occur in the North Atlantic, which is known to support high rates of N 2 fixation (Karl et al., 2002b).The model experiment that best reproduces these low values there has the highest rate of global N 2 fixation (Fig. 4).This model experiment also has the highest rate of BD, which removes fixed-N and stimulates additional N 2 fixation (Fig. 2).Note that atmospheric nitrogen deposition is not included in this model version.Since pre-industrial deposition rates are estimated to be approximately an order magnitude lower than N 2 fixation (Duce et al., 2008), it likely has a small effect in this pre-industrial scenario.However, atmospheric N deposition can introduce even more 15 N-depleted nitrogen in the North Atlantic (δ 15 N ∼ = −4 ‰) (Ryabenko et al., 2012), which could bias the model-data comparison.
The subarctic oceans contain large, shallow shelves where BD occurs.Seafloor δ 15 N show higher values (6-10 ‰) towards the shallow shelves in the Bering Sea, Sea of Okhotsk, Baffin Bay, and Banks of Newfoundland, despite the fact that less enrichment during burial occurs on these shallow shelves (Robinson et al., 2012).The model experiment with the largest net fractionation factor for BD (ε BD = 4 ‰) best reproduces the observational trends of high seafloor δ 15 N in these areas (Fig. 4), although values are still slightly underestimated.The other model experiments with smaller fractionation factors produced too low δ 15 N throughout this region.Some of the bias in this region may also be due to the shallow continental shelves that are not fully resolved in MOBI, nor the small-scale processes that take place on them.This could significantly affect surface NO − 3 utilization patterns, and thus seafloor δ 15 N. Future model versions with higher vertical resolution will evaluate this potential model bias.Nevertheless, this model-data comparison supports the view of net fractionation factors for BD to be ≥ 4 ‰, at least in the shallow subarctic ocean.

3
The global ocean mean δ 15 NO − 3 is determined by the rates and isotope effects of the source and sink terms of the fixed-N budget: N 2 fixation, WCD, and BD.N 2 fixation provides the ocean with 15 N-depleted nitrogen (δ 15 N NFix = −1.5 ‰).N-loss processes, on the other hand, preferentially consume this 15 N-depleted nitrogen, leaving the global mean nitrate pool enriched in 15 N (global mean δ 15 NO − 3 = ∼ 5.5 ‰).The average net fractionation that occurs during total N-loss determines how high the global mean δ 15 NO − 3 will be relative to the 15 N-depleted nitrogen source from N 2 fixation.We focus on two isotope effects with high uncertainty in this study: (i) OMZ nitrate utilization and dilution impacts on the isotope effect of WCD and (ii) the net fractionation factor associated with BD.

OMZ nitrate utilization and dilution isotope effect
The elevated δ 15 NO − 3 signature produced in suboxic zones depends on the level of nitrate utilization there.Utilization is determined by the balance between consumption by denitrification and replenishment by circulation and mixing.This balance determines the average δ 15 NO − 3 value that denitrifiers consume and convert to N 2 gas.The dilution effect takes into account that δ 15 NO − 3 will be weighted towards the water parcel with higher nitrate when mixing occurs (Deutsch et al., 2004).For example, if the nitrate concentration in the suboxic zone is only half of the nitrate concentration in surrounding oxic waters, the δ 15 NO − 3 signature of the oxic water will have twice the influence on total δ 15 NO − 3 of these water masses if they mix together.Note that the dilution effect is inherently simulated in the physical circulation model.
High levels of nitrate utilization reduce the influence of the isotope effect of WCD on global mean δ 15 NO − 3 .As average δ 15 NO − 3 increases in suboxic zones as denitrification occurs, the nitrate removed then becomes more 15 N-enriched.The influence of the isotope effect of WCD on global mean δ 15 NO − 3 is reduced when the removed nitrogen has a higher δ 15 N signature that is closer to global mean δ 15 NO − 3 .For example, imagine a situation with high nitrate utilization in which the average δ 15 NO − 3 value in the suboxic zone was 30.5 ‰ (instead of ∼ 12-15 ‰ in the real ocean).The nitrogen removed would then have an isotopic signature 25 ‰ depleted relative to this value, which would be equal to the global mean (δ 15 N removed = 5.5 ‰).In this case, it would have no influence on global mean δ 15 NO − 3 , even though denitrification is still fractionating the nitrogen isotopes.Thus, as nitrate utilization increases the δ 15 NO − 3 signature of the suboxic zone, it reduces the influence of WCD on global mean δ 15 NO − 3 .Experiments with high nitrate utilization in the suboxic zone require less NFix relative to WCD, and thus lower ratios of BD : WCD, to maintain global mean δ 15 NO − 3 with a balance fixed-N budget.The nitrate utilization effect alone (varying limWCD) required BD : WCD ratios that varied by nearly a factor of 2 in our range of sensitivity experiments (Table 2).BD is needed in the model to stimulate additional N 2 fixation to balance global mean δ 15 NO − 3 to the observed level.When nitrate utilization is high, the influence of the isotope effect of WCD is reduced and therefore less BD, and lower BD : WCD ratios, are required to balance global mean δ 15 NO − 3 .The model experiment that best reproduces nitrate and δ 15 NO − 3 observations in the suboxic zone is limWCD = 32 (Table 2, Fig. 5).It gives a good agreement with the amount of nitrate drawdown, as well as the slope of the increasing δ 15 NO − 3 as nitrate is consumed according to off-shelf observations in suboxic zones.The experiments with higher levels of nitrate utilization (limWCD = 20, 26) show too much nitrate consumption there.Due to deficiencies in the simulated suboxic zones, it still cannot be confirmed if the balance between nitrate consumption and replenishment is completely consistent with suboxic zones in the real ocean and the limWCD = 32 experiment.However, the high sensitivity of estimated rates of NFix and denitrification in our model experiments that test different levels of nitrate utilization suggests nitrate utilization plays an important role in global nitrogen isotope cycling.This highlights the need for higher resolution models that fully resolve all of the ventilation pathways (e.g., coastal undercurrents and eddies) of suboxic zones.

Isotope effect of benthic denitrification
Recent studies (Lehmann et al., 2007;Granger et al., 2011;Alkhatib et al., 2012) have suggested a higher net fractionation factor associated with BD (ε BD = 4-8 ‰) compared to previous estimates (ε BD = 1-3 ‰) (Brandes and Devol, 2002;Lehmann et al., 2004).They suggest BD should have a much higher net fractionation factor due to the measured high δ 15 NH + 4 that is presumably released from the sediments where BD occurs.They propose this signal is due to fractionation during the nitrification-denitrification loop in the sediments.If this high net fractionation factor is indeed cor-rect on a global scale, higher BD : WCD ratios would be required to balance the nitrogen isotope budget because additional N 2 fixation would be needed to balance this "extra" 15 N-enriched nitrogen produced in the sediments where BD occurs.However, ammonium efflux from the sediments is generally much higher on shallow shelves compared to deep ocean seafloor (Bohlen et al., 2012), suggesting that the global average fractionation of BD is likely lower than these estimates (Lehmann et al., 2007).We test the sensitivity of this effect by running experiments with ε BD set to 0, 2, and 4 ‰, while holding the limWCD parameter constant at 32, which best represented δ 15 NO − 3 observations in the suboxic zone.
The range of relative rates of BD : WCD required to closely reproduce observed global mean δ 15 NO − 3 for our sensitivity experiments (ε BD = 0-4 ‰) varied from 1.4 to 3.5 (Table 2).This large range is mostly caused by variations in BD rates and suggests that a misrepresentation of this isotope effect can significantly bias the estimate for the BD : WCD ratio.The lack of water column δ 15 NH + 4 measurements overlying sites of BD, most notably in the deep ocean, makes it difficult to constrain the global response at this time.We note our data-model analysis with the seafloor δ 15 N database supports the high estimates for the net fractionation factor of BD (ε BD ≥ 4 ‰) in the subarctic ocean that contains many shallow shelves where BD rates are high (Table 3 and Figs. 3,4).Riverine δ 15 N input is not included in this model, which can also influence some coastal settings with the input of δ 15 N between 1 and 5 ‰ (Brandes and Devol, 2002).Since riverine N input (∼ 25 Tg N yr −1 ) is relatively small compared to BD (≥ 150 Tg N yr −1 ) and introduces δ 15 N near the oceanic average, it is unlikely to have a large global impact on the ratio of BD : WCD in the pre-industrial ocean, but still may bias the model-data comparison at some locations.

Discussion
Our nitrogen isotope sensitivity experiments produce a large range of potential Nfix and denitrification rates that vary by over a factor of 2 (Table 2, Figs. 2, 6).We show that different nitrogen isotope parameters chosen for the isotope effects of WCD and BD significantly affect the estimates of the BD : WCD ratio needed to satisfy the global mean δ 15 NO − 3 constraint.This is in general agreement with previous box model studies that also estimate a large range from 1 to 3.7 (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007;Eugster and Gruber, 2012), when different isotope effects for denitrification were used in their respective models.These results show the importance of correctly modeling each isotope effect to simulate the balance of Nfix and denitrification that determines the observed global mean δ 15 NO −  , 2002;Deutsch et al., 2004;Altabet, 2007) and this study (Table 4).Each panel shows the sensitivity of one parameter on the ratio of benthic denitrification to water column denitrification (BD : WCD) needed to achieve observed global mean δ 15 NO − 3 in a steady-state scenario while all other parameters are held constant: to perform a more thorough sensitivity analysis of the key nitrogen isotope parameters.The 0D model is designed to calculate the required ratio of BD to WCD to maintain the observed global mean δ 15 NO − 3 in a steady-state pre-industrial ocean.The 0D model then calculates the BD to WCD ratio required for different nitrogen isotope parameters chosen from previous model studies and the sensitivity experiments from this study (Table 4) using Eq. ( A12).
The 0D model used here accurately reproduces the reported BD : WCD ratios of the various model configurations used in this study as well as previous studies despite the large range of model design and parameter selections (Table 4).This suggests our 0D model may be reliable to estimate the sensitivity of the different nitrogen isotope parameters in a steady-state scenario.It shows that the BD : WCD ratio needed to match global mean δ 15 NO − 3 is very sensitive to the level of OMZ nitrate utilization that determines WCD zone δ 15 NO − 3 and net fractionation factor chosen for BD, as well as other parameters (Fig. 6).The range of uncertainty for these two effects can alone account for a large range of estimates for BD (100-280 Tg N yr −1 ) from previous nitrogen isotope models (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007;Eugster and Gruber, 2012).
The large uncertainty associated with the net fractionation factor of BD adds further difficulties to constraining the BD : WCD ratio using global mean δ 15 NO − 3 .Our experiments show that increasing this factor from 0 to 4 ‰ requires almost triple the BD rate needed to maintain the global mean δ 15 NO 3 at observed levels.Recent estimates suggest that the net fractionation factor may be even higher (4-8 ‰) due to fractionation within the nitrification-denitrification loop in the sediments (Granger et al., 2011;Alkhatib et al., 2012).If these high-end estimates are validated on a global scale, this could require larger BD : WCD ratios than our largest value simulated here (> 3.5).An experiment testing ε BD = 6 ‰ in MOBI required too high BD rates to achieve a balanced global NO − 3 inventory at the modern level in the model and thus is not included here.However, Lehmann et al. (2007) show that shallow regions have a higher net fractionation compared to deep ocean sites and the global average net  Brandes and Devol (2002) reported ratio also included isotope effects from atmospheric N deposition, river input, and sediment burial, which are excluded in the 0D one-box model calculation to maintain consistency with the other model configurations.These processes slightly reduce the BD : WCD ratio and suggests the other model estimates may be slightly overestimating BD : WCD as well.
2 Altabet (2007) used a combination of reducing the fractionation factor of WCD (to account for circulation effects not included in the one-box model) and increasing WCD zone δ 15 NO − 3 so the δ 15 N value of nitrogen removed was −7 ‰.We leave ε WCD at 25 ‰ and increase the WCD zone δ 15 NO − 3 to 18 ‰ to achieve his suggested δ 15 N value of −7 ‰ for nitrogen removal.fractionation is likely closer to 4 ‰.Our 0D experiments also show that no model configuration is able to support a global net fractionation factor for BD greater than 6 ‰ and predict BD : WCD ratios in range of observational estimates (BD : WCD ≤ 4) (Fig. 6e).
MOBI experiments #4 and #5 predict BD rates on the continental shelves (60 and 108 Tg N yr −1 , respectively) that are on the low-end of most recent estimates (80-125 Tg N yr −1 ) (Bianchi et al., 2012;Bohlen et al., 2012;DeVries et al., 2013) and much lower than another estimate of 250 Tg N yr −1 (Seitzinger et al., 2006).If much more BD is occurring on the continental shelves than predicted by MOBI, this can impact the isotope effect of BD because δ 15 NO − 3 of bottom water on the shelves is on average ∼ 1.5 ‰ higher compared to the deep ocean in the model due to its close proximity to surface NO − 3 utilization.This would lead to higher values for the BD zone δ 15 NO − 3 parameter, which would require lower BD : WCD ratios (Fig. 6f).Note that this parameter has a smaller effect on BD : WCD ratios compared to the others (Fig. 6), suggesting that its uncertainty is likely lower than the other isotope parameters (e.g., net fractionation factor of BD).
The level of nitrate utilization in OMZs has a strong influence on the isotope effect of WCD.It determines the δ 15 NO − 3 value in WCD zones that is consumed by denitrifiers.Figure 6c shows the range of BD : WCD ratios required for given δ 15 NO − 3 signatures in the WCD zones for all experiments.These idealized experiments using the parameter settings of limWCD = 20, 26, and 32 (when ε BD = 0) show that the different level of nitrate utilization, and its effect on δ 15 NO − 3 in the WCD zone, is causing the range of BD : WCD ratios from 0.8 to 1.4 in these experiments.This demonstrates that if the isotope effect of WCD is not modeled accurately, it can lead to large biases of the estimates for BD : WCD.
For example, Brandes and Devol (2002) did not account for the locally high δ 15 NO − 3 of WCD zones in their one-box model (Fig. 6c).The δ 15 NO − 3 removed during WCD is thus much more 15 N-depleted compared to the other models that take into account nitrate utilization in the suboxic zone.This increases the isotope effect of WCD in the Brandes and Devel (2002) model configuration, and it thus needs more N 2 fixation to maintain global mean δ 15 NO − 3 , which is achieved by imposing a higher BD : WCD ratio.If they had accounted for a more realistic suboxic-zone δ 15 NO − 3 signature in the range of the other model configurations, our 0D model suggests their estimate for BD : WCD could have been nearly a factor of 2 lower (Fig. 6c).
Our MOBI experiment #5 (limWCD = 32, ε BD = 4) is in general agreement with the final results of Brandes and Devol (2002), suggesting a relatively high rate of BD, despite that the δ 15 N model configuration for the isotope effects of water column and BD is much different (Table 4).This suggests that even though the simple one-box model of Brandes and Devol (2002) was able to reach a similar result, it was due to a different combination of isotope effects that are difficult to constrain in a one-box model.The fact that MOBI results are directly comparable to δ 15 N observations in regions where denitrification occurs in the water column and sediments allows better validation of the various isotope effects.
The average level of nitrate utilization throughout the ocean's prominent suboxic zones remains difficult to assess.While studies in the eastern tropical North Pacific (ETNP) and Arabian Sea OMZs do not typically show nitrate depleted below half of its expected value (Brandes et al., 1998;Voss et al., 2001), recent results from the ETSP suboxic zone show more than two-thirds of the expected nitrate was consumed with a much smaller fractionation factor (Ryabenko et al., 2012).They note large rates of BD occurring within close proximity to WCD were likely contributing to this larger www.biogeosciences.net/10/5889/2013/Biogeosciences, 10, 5889-5910, 2013 nitrate deficit and reduced fractionation.Furthermore, evidence from an eddy entraining suboxic water from the ETSP OMZ showed an even larger level of nitrate utilization as it moved offshore (Altabet et al., 2012).These results suggests that an average level of nitrate utilization in the global suboxic zones may be higher than off-shore observations from the Arabian Sea and ETNP included in Fig. 4. Since MOBI does not incorporate these specific isotope effects inferred from observations in the ETSP OMZ, the model may overestimate the global net WCD isotope effect and thus BD : WCD ratios in the ocean (Altabet, 2007).
The 6-box model of Deutsch et al. (2004) accounts for local nitrate utilization with a designated suboxic box, but still estimates a higher BD : WCD ratio (2.69) compared to the results from our MOBI experiment in which the fractionation factor for BD was also set to 0 ‰ (BD : WCD = 1.4).The most significant difference between the Deutsch et al. (2004) and the other model configurations is the isotope effect of N 2 fixation.Deutsch et al. (2004) chose a δ 15 N NFix signature of 0 ‰, while all other models selected between −1 and −1.5 ‰.If the Deutsch et al. (2004) study would have chosen the same value as here (δ 15 N NFix = −1.5 ‰), our 0D model suggests this would decrease their estimated BD : WCD ratio from 2.69 to 1.83, which would then be more consistent with the results from MOBI.This shows that even small differences (< 2 ‰) for the isotope effect of N 2 fixation can alter the ratio of BD : WCD by 30 % or even more depending on the model configuration used (Fig. 6d).
The MOBI experiments that most closely reproduce seafloor δ 15 N observations are experiments #4 and #5 (limWCD = 32, ε BD = 2 and 4), with ε BD = 2 producing slightly better correlation and root mean squared error statistics.(Table 2).They predict a range of rates for N 2 fixation,WCD,76, Tg N yr −1 , respectively.These experiments produce a large range of BD : WCD ratios from 2 to 3.5 and highlight the high sensitivity of the BD : WCD ratio to the net fractionation factor of BD.Although the average level of nitrate utilization in the suboxic zones is uncertain, our experiments using limWCD = 32 best represent observations from the ETNP and Arabian Sea.Assuming this range for nitrogen isotope parameters, our model estimates a potential range for BD : WCD of 2.0-3.5.
Our model experiments are in general agreement with a recent 3D inverse model that included nitrogen isotopes to constrain marine denitrification rates (DeVries et al., 2013).They similarly show a high sensitivity to the NO − 3 utilization in the suboxic zone and fractionation factor assumed for BD.However, they estimate lower ratios of BD : WCD from 1.3 to 2.3 compared to our results with MOBI.The main reasons for this discrepancy are likely that DeVries et al. (2013) assumed a slightly higher level of nitrate utilization in the suboxic zones and lower values for ε BD = 0-3 ‰, whereas our high-end estimate for BD : WCD = 3.5 is due to using ε BD = 4 ‰.The high sensitivity to these parameters empha-sizes the need to better understand and quantify them in future studies.

Conclusions
Our study uses water column δ 15 NO − 3 and seafloor δ 15 N observations to constrain the rates of N 2 fixation, WCD, and BD in the global ocean.The uncertainty associated with isotope effects of denitrification in the water column and sediments makes it difficult to constrain N 2 fixation and total denitrification rates.Previous box model studies using δ 15 N have estimated a large range for the ratio of BD : WCD from 1 to 3.7 (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007;Eugster and Gruber, 2012).Here we used a set of experiments with a global coupled three-dimensional circulation-biogeochemistry-isotope model (MOBI) and a one-box model to show that nitrate utilization in the suboxic zone and the net fractionation factor of BD, both of which are not well constrained by observations, can lead to rates of BD that vary by over a factor of 2 if not modeled correctly.
With our global coupled three-dimensional model, we are able to compare δ 15 N observations in the water column and seafloor in the regions where denitrification occurs to constrain the nitrogen isotope parameters in the model.This highlights the importance of using models that can resolve all of the locally important nitrogen isotope effects that affect δ 15 N in denitrification zones.The model experiments that best reproduce δ 15 N observations in the water column and sediments estimate the rates of N 2 fixation, WCD, and BD in the ranges of 195-350, 65-80, and 130-270 Tg N yr −1 , respectively, assuming a balanced fixed-N budget in the preindustrial ocean.Although uncertainties still exist, this model result suggests that N 2 fixation is occurring at much greater rates than previously estimated, and the residence time for oceanic fixed nitrogen is between ∼ 1500 and 3000 yr.

Nitrogen isotope model description A1 Fractionation equation
Fractionation is calculated using kinetic fractionation (Mariotti et al., 1981): where α is the kinetic fractionation factor associated with the process and the N pro and N sub refer to the nitrogen of the product and substrate of the reaction, respectively.In the model, we include 15 N as the prognostic variable instead of the ratio 15 N/ 14 N.The 15 N equations are embedded within the marine ecosystem model that calculates total N ( T N = 15 N+ 14 N).Solving Eq. (A1) for 15 N pro with respect to T N pro yields where R sub is the isotopic ratio 15 N/ 14 N of the substrate of the reaction.This equation can be equivalently expressed in the commonly used delta ("δ") notation by applying the relation (Mariotti et al., 1981): which gives positive values for ε with this definition.Eq. (A2) then becomes where et al., 2000), which is the nitrogen isotope equation coded into MOBI.Note we use a R std value of 1 so that 14 N and 15 N have concentrations of the same order of magnitude.This reduces the impact of numerical noise caused by the advection scheme on the δ 15 N value.If the atmospheric N 2 ratio was used (R std = 0.0036765), the 15 N concentration would be over two orders of magnitude smaller and be more susceptible to numerical noise, which produces erroneous δ 15 N values.In the polar oceans (> 80 • N/S) where numerical noise is the highest in our global model, some model grid points still contain erroneous isotope values so this region is not included in the statistical analysis.

A2 Coupled model equations
The fractionation equation used for NO − 3 consumption during phytoplankton uptake and WCD follows Eq. (A4) where 1000] and u is the fraction of available total nitrate consumed during each time step.This fractionation equation is used to ensure that if a significant portion of the nitrate pool is consumed in one time step, mass balance of the different nitrogen isotope species is conserved.In the experiments here, we artificially limit WCD at high enough nitrate concentrations (26-32 µM) so this term [(1 − u)/u • ln(1 − u)] has a negligible effect for WCD in this study.Since zooplankton excretion and BD are parameterized in the model, the instantaneous fractionation equation is used (Eq.A4) with a given fractionation factor to mimic the net fractionation that occurs during the integrated reaction.
The full set of time-dependent equations for 15 N that are embedded into the marine ecosystem biogeochemical model are as follows: and where T R X = 15 N X /( 15 N X + 14 N X ).Here it suffices to note that the equations for total nitrogen ( 14 N + 15 N) are identical to the ones for 15 N, except that T R X = β X /(1 + β X ) = 1 in the total nitrogen equations.The parameter list is given in Table B1.where β = α • R sub = R sub (1 − ε/1000), consistent with Eq. (A4).Solving for benthic to WCD ratio yields These results are displayed in Table 4 and Fig. 6.

Appendix B Marine ecosystem-biogeochemistry model description
This appendix provides a description of the parameters used in the full set of time-dependent equations in the marine ecosystem model.It suffices to note that the equations for to-tal nitrogen ( 14 N + 15 N) ecosystem variables are identical to the ones of 15 N if R X = β X /(1 + β X ) = 1, which are located in Appendix A. The function J * O provides the growth rate of nondiazotrophic "ordinary" phytoplankton, determined from irradiance (I ), NO − 3 and PO 3− 4 .
The maximum growth rate is dependent only on temperature (T ), ( * ).We use a 0 = 0.35 d −1 for the maximum growth rate at 0 • C which was determined to optimize surface nutrient concentrations.Under nutrient-replete conditions, the lightlimited growth rate J OI is calculated according to where α is the initial slope of the photosynthesis vs. irradiance (P-I) curve.The calculation of the photosynthetically active shortwave radiation I and the method of averaging Eq. (B3) over one day is outlined in Schmittner et al. (2005).This version also includes the correction for the error in the calculation of light limitation in previous versions (Schmittner et al., 2008).Nutrient limitation is represented by the product of J * Omax and the nutrient uptake rates, u N = NO − 3 /(k N + NO − 3 ) and u P = PO 3− 4 /(k P + PO 3− 4 ), with k P = k N r P : N providing the respective nutrient uptake rates.
Diazotrophs grow according to the same principles as the ordinary phytoplankton class, but are disadvantaged in nitrate-bearing waters by a lower maximum growth rate, J * Diazmax : The coefficient c Diaz handicaps diazotrophs by dampening the increase of their maximal growth rate versus that of the general phytoplankton class with rising temperature.We use c Diaz = 0.13, such that the growth rate of diazotrophs is 13 % that of ordinary phytoplankton.This handicap is further decreased by the Fe limitation parameter, which is scaled between 0 and 1 by multiplying a monthly climatology of aeolian dust deposition (Mahowald et al., 2005b(Mahowald et al., , 2006(Mahowald et al., , 2009) by a constant factor and setting the maximum value to 1 (Somes et al., 2010a).However, diazotrophs have an advantage in that their growth rate is not limited by NO − 3 concentrations, J * Diaz (I, PO 4 )= min(J * DiazI , J * Diazmax u P ), (B5) although they do take up NO − 3 if it is available (see term 1 on the right hand side of Eq.A7).The N : P of model diazotrophs is set to 40 : 1.
The first-order mortality rates of phytoplankton and diazotrophs are linearly dependent on their biomass concentrations, P O and P Diaz .Dissolved organic matter and the microbial loop are folded into a single fast-remineralization process, which is the product of their biomass and the temperature dependent term.where g is grazing rate, ω is prey-capture rate, and P is phytoplankton concentration (Table B1).Note prey-capture rate is reduced for diazotrophs relative to ordinary phytoplankton in these experiments (Table 1).
Since diazotrophs have a higher N : P ratio Remineralization transforms the N and P content of detritus to NO − 3 and PO 3− 4 .Photosynthesis produces oxygen, while respiration consumes oxygen, at rates equal to the consumption and remineralization rates of PO 4 , respectively, multiplied by the constant ratio R O : P .Dissolved oxygen exchanges with the atmosphere in the surface layer (F sfc ) according to the Ocean Carbon-Cycle Model Intercomparison Project protocol (Orr, 1999).
Oxygen consumption in suboxic waters (O 2 < ∼ 5 µM) is inhibited, according to (B15) BD is the rate at which nitrate is removed from the bottom water.We assume that the rain rate of carbon into the sediments occurs at a ratio of R C : N = 6.625 of the nitrogen in the sinking organic detritus.
Since the continental shelves and other small-scale bathymetric features are not well resolved in the model, we use a subgrid-scale parameterization.The portion of each bottom ocean grid box that is deeper than the real sea floor is calculated at each location from high-resolution (1/5 • ) bathymetry.The rain rate of carbon that is included in the BD function in this shelf parameterization is the amount of particulate organic carbon that sinks into the portion of the grid box covered by a shallower continental shelf.In the model, ∼ 30 % of BD occurs within this shelf parameterization.The remaining particulate organic matter continues to sink to greater depths.The coarse-resolution physical circulation model's inability to fully resolve coastal systems generally underestimates primary production and sinking carbon fluxes on these continental shelves, which likely results in too-low BD rates there.To account for this deficiency, we multiply the BD transfer function by an arbitrary coefficient (α BD ).This parameter is tuned to set the global deep ocean δ 15 NO − 3 in the model to ∼ 5 ‰ for each experiment.Figure 2 shows the spatial distribution of BD.

Fig. 6 .
Fig. 6. Results of the one-box δ 15 NO − 3 model with different nitrogen isotope parameter settings according to previous configurations (Bran-des and Devol, 2002;Deutsch et al., 2004;Altabet, 2007) and this study (Table4).Each panel shows the sensitivity of one parameter on the ratio of benthic denitrification to water column denitrification (BD : WCD) needed to achieve observed global mean δ 15 NO − 3 in a steady-state scenario while all other parameters are held constant: (a) WCD rate, (b) WCD fractionation factor (ε WCD ), (c) average δ 15 NO − 3 value where WCD occurs, (d) N 2 fixation fractionation factor (ε NFix ), (e) benthic denitrification fractionation factor (ε BD ), and (f) average δ 15 NO − 3 value where benthic denitrification occurs.Note ε NFix makes diazotroph biomass depleted in δ 15 N by its associated value.The crosses (X) show the parameter chosen for each reported model configuration and denotes the BD to WCD ratio of the previous studies as well as coupled three-dimensional experiments from this study.
Fig. 6. Results of the one-box δ 15 NO − 3 model with different nitrogen isotope parameter settings according to previous configurations (Bran-des and Devol, 2002;Deutsch et al., 2004;Altabet, 2007) and this study (Table4).Each panel shows the sensitivity of one parameter on the ratio of benthic denitrification to water column denitrification (BD : WCD) needed to achieve observed global mean δ 15 NO − 3 in a steady-state scenario while all other parameters are held constant: (a) WCD rate, (b) WCD fractionation factor (ε WCD ), (c) average δ 15 NO − 3 value where WCD occurs, (d) N 2 fixation fractionation factor (ε NFix ), (e) benthic denitrification fractionation factor (ε BD ), and (f) average δ 15 NO − 3 value where benthic denitrification occurs.Note ε NFix makes diazotroph biomass depleted in δ 15 N by its associated value.The crosses (X) show the parameter chosen for each reported model configuration and denotes the BD to WCD ratio of the previous studies as well as coupled three-dimensional experiments from this study.
by the oxygen-equivalent oxidation of nitrate,r NO 3 sox = 0.5 [1 − tanh (O 2 −5)] .(B14)Denitrification consumes nitrate at a rate of 80 % of the oxygen equivalent rate, as NO − 3 is a more efficient oxidant on a mol per mol basis (i.e., one mol of NO − 3 can accept 5 e − while 1 mol of O 2 can accept only 4 e − ).We include the scheme of Bohlen et al. (2012), which parameterizes BD based on the rain rate of POC (RR POC ) into the seafloor and bottom water oxygen and nitrate: BD = α BD 0.09782 + 0.22944 × 0.9811 bwO 2 −NO − 3 ×RR POC .

Table 2 .
Model sensitivity experiment results.
Model results for the different sensitivity experiments in the global ocean and suboxic zone (O 2 < 10 µM) with the observational estimate given in parentheses.
* Statistical measures are the correlation coefficient (R