Marine denitriﬁcation rates determined from a global 3-D inverse model

. A major impediment to understanding long-term changes in the marine nitrogen (N) cycle is the persistent uncertainty about the rates, distribution, and sensitivity of its largest ﬂuxes in the modern ocean. We use a global ocean circulation model to obtain the ﬁrst 3-D estimate of marine denitriﬁcation rates that is maximally consistent with available observations of nitrate deﬁcits and the nitrogen isotopic ratio of oceanic nitrate. We ﬁnd a global rate of marine denitriﬁcation in suboxic waters and sediments of 120–240 Tg N yr − 1 , which is lower than many other recent estimates. The difference stems from the ability to represent the 3-D spatial structure of suboxic zones, where denitriﬁcation rates of 50–77 Tg N yr − 1 result in up to 50 % depletion of nitrate. This depletion reduces the effect of local isotopic enrichment on the rest of the ocean, allowing the N isotope ratio of oceanic nitrate to be achieved with a sedimentary denitriﬁcation rate about 1.3–2.3 times that of suboxic zones. This balance of N losses between sediments and suboxic zones is shown to obey a simple relationship between isotope fractionation and the degree of nitrate consumption in the core of the suboxic zones. The global denitriﬁcation rates derived here suggest that the marine nitrogen budget is likely close to balanced.


Introduction
Relative to the cycles of other biologically important nutrients, the marine nitrogen cycle is potentially highly dynamic, with large input and output rates and a relatively short turnover time. A question of central importance is whether the marine nitrogen budget can sustain long-term imbalances, or if the primary sources and sinks of nitrogen are closely coupled by self-stabilizing feedbacks, preventing significant variability in the ocean's nitrogen inventory. The answer remains unclear in large part because the nitrogen budget of the contemporary ocean is poorly constrained, primarily due to uncertainty in the rate of denitrification, by which fixed nitrogen is converted to N 2 gas and lost from the ocean. Estimates of denitrification rates in the contemporary ocean range from about 200 Tg N yr −1 Sarmiento, 1997, 2002; to over 400 Tg N yr −1 (Middelburg et al., 1996;Codispoti et al., 2001;Brandes and Devol, 2002;Codispoti, 2007).
Denitrification occurs both in small areas of the ocean where waters become suboxic (water-column denitrification), and in the pore-waters of sediments throughout the ocean (benthic denitrification). Recent work suggests that the rate of water-column denitrification is about 60-70 Tg N yr −1 (DeVries et al., 2012), but the rate of benthic denitrification remains poorly known. Direct measurements of denitrification rates in sediments have been made at a handful of sites (e.g., Devol and Christensen, 1993;Devol et al., 1997;Laursen and Seitzinger, 2002;Rao et al., 2007), but scaling these up to a global estimate is precluded by the sparsity of observations and the spatial and temporal variability in these rates.
One solution to the challenge of deriving a global estimate of marine denitrification rates is to make use of observed marine nitrate (NO 3 ) deficits and nitrogen isotopic ratios (Fig. 1). These quantities are more widely measured than denitrification rates, and tend to reflect the integrated effects of highly variable denitrification processes. Nitrate deficits can be measured by the N * tracer, which reflects the difference between the in situ NO 3 concentration and that expected due to the average nitrate to phosphate (PO 4 ) ratio of organic matter, N * = NO 3 − 16×PO 4 . The isotopic ratio of 15 N to 14 N in oceanic nitrate (R = 15 NO 3 / 14 NO 3 ) is commonly expressed as δ 15 NO 3 = (R/R std − 1) × 1000, where R std is the isotopic ratio of atmospheric N 2 .
In tandem, N * and δ 15 NO 3 provide powerful constraints on marine denitrification rates. Denitrification in both the water column and the sediments consumes nitrate but not phosphate, imparting a negative signature to N * . The influence of denitrification is clearly visible in the thermocline N * distribution, which shows strongly negative N * due to watercolumn denitrification in suboxic waters of the Arabian Sea and the Eastern Tropical Pacific, and negative N * due to benthic denitrification in the sub-Arctic Pacific and elsewhere (Fig. 1b). Water-column denitrification preferentially removes the lighter nitrogen isotope, and imparts a heavy isotopic signature of w ∼ 25‰ (Sigman et al., 2003) to surrounding waters, which explains the elevated δ 15 NO 3 values in the suboxic waters of the eastern tropical Pacific and the Arabian Sea (Fig. 1d). However, benthic denitrification has a much lighter isotopic signature of 0‰ ≤ b 3‰ Devol, 1997, 2002;Lehmann et al., 2004Lehmann et al., , 2007, as evidenced by the lack of isotopic enrichment in the sub-Arctic Pacific (e.g., Yoshikawa et al., 2006) (Fig. 1d).
The mean ocean δ 15 NO 3 of about 5‰ primarily reflects a balance between the input of isotopically light NO 3 by nitrogen fixers (with an isotopic signature of −2‰ fix ≤ 0‰ (Macko et al., 1987;Carpenter et al., 1997)) and a mixture of water-column and benthic denitrification, and therefore provides a strong constraint on the relative amounts of denitrifi-cation occurring in the water column (W ) and the sediments (B), which can be expressed compactly as the ratio B/W . Simple geochemical box models have been used to derive estimates of B/W , but these estimates still yield a large uncertainty of 1 < B/W < 4 (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007;Eugster and Gruber, 2012). The uncertainty in B/W partly reflects uncertainty in isotopic enrichment factors and the isotopic ratio of organic nitrogen (Brandes and Devol, 2002;Altabet, 2007). More importantly however, the wide range of estimates of B/W reflects inaccuracies associated with simple box models, which cannot resolve important spatial features of the ocean circulation and denitrification processes. In order to correctly account for the isotopic effects of nitrate consumption in the suboxic water column, a model should correctly simulate the degree of nitrate consumption in suboxic zones, and accurately simulate how the suboxic zones are ventilated and how tracers are exchanged between the suboxic and oxic ocean (e.g., Deutsch et al., 2004). These latter effects can only be accurately captured in a spatially explicit 3-D ocean circulation model. The only previous study to simulate nitrogen isotopes in a global ocean circulation model suffered from severe overconsumption of nitrate in the suboxic water column (Somes et al., 2010), rendering the B/W estimate from that model inaccurate.
Here we address these issues by coupling a simple nitrogen cycle model to a data-constrained ocean circulation model (Sect. 2.1). The parameters of the nitrogen cycle model are iteratively adjusted using an adjoint approach to achieve an optimal fit to the observed distributions of N * and δ 15 NO 3 (Sects. 2.2 and 2.3). The solution to this global inverse nitrogen cycle model is an estimate of the global rates and patterns of water-column and benthic denitrification (Sect. 3). The effects of denitrification on the N * distribution (Sect. 4.1) and on the mean ocean δ 15 NO 3 (Sect. 4.2) are discussed. Model parameters that are poorly constrained by the available data are varied by a Monte Carlo procedure in order to derive uncertainty estimates on denitrification rates (Sect. 4.3). We also discuss the implications of our findings for the global marine nitrogen budget (Sect. 5).

Circulation model and nitrogen cycle model
The physical component of the nitrogen cycle model is based on the data-constrained model of DeVries and Primeau (2011), which has been extended to increase the resolution of the model to 2 • in the horizontal, with 24 vertical levels. As in DeVries and Primeau (2011), the circulation of the model has been tuned to closely reproduce the observed temperature, salinity, and radiocarbon distributions in the ocean using an adjoint method. We also add a simple simulation (Najjar and Orr, 1998) of the phosphate and dissolved organic phosphorus (DOP) cycle to the model, with the observed PO 4 distribution (Garcia et al., 2010a) as an additional constraint on the circulation and biological fluxes. Unlike the study of DeVries et al. (2012), we do not assimilate  observations in this version of the model.
The annual mean steady-state circulation determined in the above step is then taken offline and used to compute the physical tracer transport in a simple marine nitrogen cycle model. The internal cycling of N is driven by restoring surface nitrate toward observations, as done for PO 4 . This ensures that the model reproduces the observed nutrient stoichiometry (N * ) of surface waters. One of the necessary chemical signatures of N 2 fixation -a non-Redfield uptake of NO 3 and PO 4 -is therefore already implicit in the model design. Thus the N * observations cannot simultaneously constrain the rates and spatial distribution of N 2 fixation. Rather, we model N 2 fixation according to a simple dependence on light, temperature, and nutrient availability (see Appendix A), and use it primarily to close the N budget. An inverse solution that solves simultaneously for N 2 fixation and denitrification rates requires an explicit treatment of non-Redfield stoichiometry of nutrient uptake, as well as riverine and atmospheric N inputs, and is left for a future study.
Water-column denitrification in the model occurs where observed oxygen concentrations fall below a critical threshold O 2, crit . Observed oxygen concentrations are taken from the 2009 World Ocean Atlas monthly climatology (Garcia et al., 2010b) after correcting for measurements in lowoxygen regions (Bianchi et al., 2012). Benthic denitrification occurs within grid cells that have some contact with the ocean floor. The grid cells having contact with the ocean floor are determined by interpolating observed bathymetry to the model grid, as described in Appendix A. Water-column denitrification rates are proportional to the rate of organic matter remineralization within suboxic zones, while benthic denitrification rates are proportional to the rate of organic matter supply to the sediments, and also depend on bottom-water nitrate and oxygen concentrations. Denitrification in the water column and sediments is balanced by nitrogen fixation, with a rate that depends on local surface NO 3 concentrations, temperature, light levels, and iron supply. Production of organic nitrogen is parameterized by restoring to observed NO 3 in the top two model layers (above 73 m depth). Organic nitrogen is exported out of the surface layers as either dissolved organic nitrogen (DON) which remineralizes with first-order kinetics, or as particulate organic nitrogen (PON), which remineralizes according to a power-law dependence on depth (Martin et al., 1987). See Appendix A for a full description of the nitrogen cycle model.

Inverse nitrogen cycle model
The parameters of the nitrogen cycle model include the critical oxygen threshold for water-column denitrification; the ratio of nitrate consumed to organic matter remineralized during water-column denitrification, and that during benthic denitrification; the oxygen and nitrate dependence of benthic denitrification; the isotopic enrichment factors for nitrogen fixation, water-column denitrification, benthic denitrification, and uptake of nitrate to form organic nitrogen; the maximum nitrogen fixation rate as well as its temperature, light, nitrate, iron, and depth dependence; the fraction of organic matter production routed to the dissolved organic nitrogen (DON) pool; and the decay timescale for DON. Most of the model parameters can be constrained by the N * and δ 15 NO 3 data, and these parameters (see Table B1) are iteratively adjusted using an adjoint method to find values that best fit observed N * and δ 15 NO 3 (Fig. 1). See Appendix B for a full description of the inverse model and the formulation of the cost function measuring the model-data misfit.
We withhold four of the model parameters from the inversion because they are likely to be unconstrained by the available N * and δ 15 NO 3 data. Primarily due to the sparsity of the δ 15 NO 3 observations, only one of the isotopic enrichment factors can be constrained independently of the others. We solve for w as part of the solution to the inverse model, since there is good data coverage in the suboxic zones where water-column denitrification occurs (Fig. 1c). We fix the isotopic enrichment factor for uptake of nitrate to form organic matter ( up ) at 5‰ as in previous studies (e.g., Somes et al., 2010;Eugster and Gruber, 2012). We account for uncertainty in the remaining isotopic fractionation factors by rerunning the inverse model with various combinations of b (0, 1, 2, or 3‰) and fix (−2, −1, or 0 ‰). The fraction of organic nitrogen routed to the DON pool (σ DON ) is also poorly constrained, since observations of DON concentrations are not used to constrain the model. So we also vary the value of σ DON over a wide range (σ DON =1/4, 1/3, 1/2, or 2/3) in different versions of the inverse model. In all, the different combinations of w , fix , and σ DON produce 48 configurations for the inverse model. The uncertainty ranges quoted below represent the full range of optimal model solutions under each of these 48 different model configurations.
Data constraints for the inverse model include N * concentrations derived from the 2009 World Ocean Atlas objectively mapped annual mean NO 3 and PO 4 data (Garcia et al., 2010a), which are interpolated to the model grid, and δ 15 NO 3 observations compiled from the literature (Somes et al., 2010;Rafter et al., 2012;De Pol Holz et al., 2009;Liu, 1979;DiFiore et al., 2009;Yoshikawa et al., 2005, 2006, P. Rafter andand unpublished data), which are binned to the model grid. We exclude δ 15 NO 3 observations above 200 m depth, due to numerical noise in the simulated δ 15 NO 3 fields that can occur near the surface where nitrate concentrations are very low (close to zero).

Model-data comparison
The optimization procedure produces a good fit to the observed distributions of N * (Fig. 2b) and δ 15 NO 3 (Fig. 2c).

Fig. 2. (a-c)
Joint distribution function for the gridbox-volumeweighted observed and modeled tracer concentrations. The joint distribution function was estimated using the kernel density estimation method described in Botev et al. (2010), with a modification to account for the volume of the model grid boxes (Primeau et al., 2013). Printed on each plot is the gridbox-volume-weighted mean model-data difference and the gridbox-volume-weighted root mean squared error. Although NO 3 is not included in the cost function measuring model-data misfit, a good fit to observed NO 3 (Fig. 2a) is achieved by virtue of the fact that both PO 4 and N * are included in the cost function. The mean modeled nitrate and N * concentrations are slightly lower than the observed values, and the mean ocean δ 15 NO 3 is around 5 ‰, in excellent agreement with the observations. The model also demonstrates good agreement with the observed degree of nitrate consumption (f c = 1-NO 3 /16PO 4 ) and δ 15 NO 3 in lowoxygen zones (Fig. 2d). The degree of nitrate consumption in waters with less than 20 µM O 2 in the model is between 0.1-0.5, in agreement with observations. At high degrees of nitrate consumption (f c ∼ 0.4), the modeled δ 15 NO 3 is about 12-16 ‰, depending on the oceanic region. This agrees with the mean observed δ 15 NO 3 in these regions, although the observations show more scatter than the modeled values. This could be due to spatial or temporal variability that is not captured by the coarse steady-state model. The model also matches the depth distribution of N * and δ 15 NO 3 in both the Atlantic and Indo-Pacific basins quite well (Fig. 3). Depth profiles of modeled and observed δ 15 NO 3 in the Atlantic show that the model does not pro- duce quite low enough δ 15 NO 3 near the surface (Fig. 3a). Atmospheric deposition, which is not accounted for in the model, may play an important role in determining the shallow nitrate δ 15 NO 3 in the Atlantic (Knapp et al., 2008). Deep ocean values in the Atlantic are slightly lower than observed, but the observations lie within the envelope of model uncertainty. The modeled depth profile of δ 15 NO 3 in the Indo-Pacific basin matches the observations fairly well throughout the water column (Fig. 3a). The greatest mismatch occurs above 500 m, where modeled δ 15 NO 3 is lower than observed. This could potentially indicate too little N fixation in the model in regions with δ 15 NO 3 observations, primarily suboxic zones. N * concentrations in the Atlantic show a general decrease with depth in both the model and the observations (Fig. 3b). The model does produce a thermocline maximum in N * in the Atlantic, although not as pronounced as observed. Except in the deepest Atlantic, where the model slightly underpredicts N * , the observed mean N * lies within the model-predicted uncertainty. The model matches the general magnitude and shape of the observed N * depth profile in the Indo-Pacific, with some distinctive differences most notably in the depth range 500-1000 m where the model has a stronger and shallower negative N * peak than the observations, and in the deep ocean where the model N * is more negative than observed (Fig. 3b). However, the largest model-data misfit (∼ 1 µM) is small compared to the N * signature of denitrification (see for example Fig. 5a and discussion below).

Denitrification rates and their global distribution
The global distribution of water-column and benthic denitrification ( Fig. 4) was determined from the optimal solution under each of the 48 different model configurations. Watercolumn denitrification best fits observations when confined primarily to low-oxygen zones (O 2 5 µM) in the eastern tropical North and South Pacific, and in the Arabian Sea (Fig. 4a). Small (negligible) amounts of water-column denitrification occur in the sub-Arctic Pacific and in the Bay of Bengal. Globally, the rate of water-column denitrification The N : P ratio of remineralized organic matter calculated according to the method of Anderson and Sarmiento (1994) using the modeled NO 3 and PO 4 fields (red circles with error bars), and the ratio after correcting for the amount of NO 3 that is lost due to denitrification (black circles with error bars). All calculations are performed on isopycnal surfaces and plotted at the mean depth of each isopycnal. predicted by the model is 50-77 Tg N yr −1 (Table 1). Of this total, 9-14 Tg N yr −1 is due to denitrification in the Arabian Sea, and 41-63 Tg N yr −1 is due to denitrification in the eastern tropical Pacific (Table 1). These rates agree, within their uncertainty, with an independent estimate of fixed nitrogen loss rates based on N 2 gas observations from within these same suboxic zones (DeVries et al., 2012). This agreement occurs despite the fact that the present study does not explicitly account for fixed nitrogen loss due to anaerobic ammonium oxidation (anammox), while the estimate of DeVries et al. (2012) does implicitly include the effects of annamox. This agreement indicates that either anammox is similar to traditional denitrification in its imprint on N * and δ 15 NO 3 , or that annamox represents a small sink of fixed nitrogen within suboxic zones. Benthic denitrification rates are highest in highly productive coastal upwelling regions (e.g., the eastern tropical Pacific, West African coast, and Arabian Sea), and in regions of shallow continental shelves (e.g., the sub-Arctic Pacific, northeast North America, Indonesian Archipelago, Arctic margin, and southeastern South America) (Fig. 4b). These areas experience benthic denitrification rates over 100 times greater than rates typical of deep ocean sediments. Globally, the model predicts that 70-170 Tg N yr −1 of benthic denitrification is needed to match the constraints provided by the N * and δ 15 NO 3 data, with the largest contribution from sediments in the Pacific Ocean (Table 1).
The vertical distribution of water-column denitrification is highly concentrated and confined above 1000 m depth, coincident with the depth of suboxic zones (Fig. 4c). Watercolumn denitrification rates have a shallower peak than the suboxic volume (Fig. 4c) due to the fact that organic matter remineralization rates decrease with depth. Benthic denitrification also has a shallow peak, with highest rates above 1000 m depth, although approximately half of benthic denitrification occurs below 1000 m (Fig. 4d). The primary factor controlling the distribution of benthic denitrification is the rate at which organic matter is delivered to the sediments. However, the data imply a mid-depth peak in benthic denitrification at about 1000-2000 m that is not found in the rate of organic matter supply to the sediments (Fig. 4d), which the model ascribes to enhanced benthic denitrification rates under low oxygen and high nitrate conditions (e.g., Middelburg et al., 1996;Bohlen et al., 2012).
The rates of benthic denitrification reported here are lower than some earlier estimates, which were on the order of ∼ 250-300 Tg N yr −1 (Middelburg et al., 1996;Seitzinger et al., 2006). However, more recent estimates are in better agreement with the benthic denitrification rates derived here. Bianchi et al. (2012) used the meta-model parameterization of Middelburg et al. (1996) and satellite-derived estimates of sinking organic matter flux to infer a benthic denitrification rate of 182 ± 55 Tg N yr −1 , while Bohlen et al. (2012) derived data-based transfer functions to scale satellite-derived estimates of sinking organic matter flux to estimates of benthic denitrification, arriving at a globally integrated benthic denitrification rate of ∼ 155 Tg N yr −1 . These more recent estimates agree within uncertainty with the benthic denitrification rate found in this study (70-170 Tg N yr −1 ), although our estimate is, on average, slightly lower.
A significant difference between this study and previous studies is that our model predicts that only about 20 % of benthic denitrification occurs in shelf sediments (< 160 m depth), in contrast with other estimates suggesting 35-70 % of benthic denitrification occurs on continental shelves (Middelburg et al., 1996;Bianchi et al., 2012;Bohlen et al., 2012). Although we have used a parameterization to account for the presence of shallower continental shelves than are resolved by the model's bottom topography (see Appendix A), the model resolution is still insufficient to resolve the enhanced biological productivity and high organic carbon export rates within shelf regions. Because the model underestimates benthic denitrification on continental shelves, and in order to maintain the correct ratio of benthic to water-column denitrification, the model probably slightly overpredicts the amount of denitrification in the deep ocean. That there is too much benthic denitrification in deep sediments can be seen in Fig. 3, which shows that model-predicted N * is about 20 % lower than observed in the deep ocean. If denitrification in shelf areas produces less of an imprint on the mean ocean δ 15 NO 3 than denitrification in deep sea sediments, this could produce a slightly low bias in the total benthic denitrification rate found by the model.
A further consideration is that riverine N sources are ignored in the model. An additional input of nitrate to the continental shelf areas could be supplied by rivers, which are estimated to deliver about 30 Tg N yr −1 into coastal waters (Gruber and Galloway, 2008), with an isotopic enrichment of about 4 ‰ (Brandes and Devol, 2002). If this were balanced entirely by denitrification on the shelf areas, this would also serve to increase the proportion of benthic denitrification occurring on continental shelves.

Effect of denitrification on vertical N * distribution
Water-column and benthic denitrification have distinct effects on the vertical distribution of N * in the ocean. Quantifying these effects can help to reveal the imprint of watercolumn and benthic denitrification on N * , and allow us to judge the relative misfit between the modeled and observed N * depth profiles (Fig. 3). To estimate the effect of denitrification on N * , we simulate an idealized denitrification tracer that is produced at a rate of 1 mol tracer per 1 mol NO 3 consumed by denitrification. The tracer is immediately removed from the ocean when it reaches the top model layer, where all nutrients are considered "preformed". We perform separate calculations for benthic and water-column denitrification. For comparison with the study of Anderson and Sarmiento (1994), we calculate the average amount of tracer on isopycnal horizons (below 400 m depth) and plot the results as a function of the average depth of each isopycnal horizon (Fig. 5a). The results show that N * due to watercolumn denitrification reaches a minimum in the thermocline at about 700-800 m depth, with a peak value of around −1 µM (Fig. 5a). Benthic denitrification produces an N * profile with a mid-depth peak of about −3 to −4 µM at 1500-3000 m depth (Fig. 5a). These are the globally averaged effects of denitrification; the effect is about 50 % larger in the Indo-Pacific Basin, where most denitrification occurs and where waters are older, allowing for accumulation of reaction products. From this calculation, we see that the misfit between modeled and observed N * depth profiles (Fig. 3b) is relatively small compared to the denitrification signal.
Denitrification also affects the nitrogen to phosphorus (N : P) ratio of subsurface waters, an effect which must be Table 1. Integrated denitrification rates by ocean basin. Median rate is given, with range in parentheses. The range includes uncertainty due to uncertainty in the isotopic enrichment factors of fixation and benthic denitrification, and in the parameterization of dissolved organic matter. Southern Ocean is south of 34 • S. Global rates also include contribution from benthic denitrification in the Arctic Ocean and the Mediterranean Sea. accounted for when determining N : P ratios of regenerated organic matter from subsurface nutrient concentrations. We calculated the N : P ratio of remineralized organic matter in the interior ocean in our model following a procedure similar to that of Anderson and Sarmiento (1994). The amount of "preformed" nitrate (or phosphate) in the interior ocean is calculated from an idealized tracer with a surface concentration equal to the modeled surface concentration, and no sources or sinks in the interior ocean. The amount of remineralized nitrate (or phosphate) is then determined by subtracting the preformed component from the total concentration. We determine the N : P ratio of remineralized organic matter from a linear regression of remineralized nitrate against remineralized phosphate along the same isopycnal horizons used in Fig. 5a. Similar to the results of Anderson and Sarmiento (1994), we see that the near-surface N : P ratio calculated in this way is close to that of "fresh" organic matter (∼ 16 : 1), but drops to a minimum in the depth range 1500-3000 m (Fig. 5b). Anderson and Sarmiento (1994) hypothesized that the actual N : P ratio of remineralized organic matter was approximately constant with depth, but that the mid-depth minimum may be an artifact of benthic denitrification. By adding back in the remineralized NO 3 that is lost due to denitrification (Fig. 5a) and repeating the calculation, we find that indeed the N : P ratio of remineralization is approximately constant with depth (Fig. 5b), which is consistent with the original hypothesis of Anderson and Sarmiento (1994). The deepest isopycnals in Fig. 5b are associated with deep Labrador Sea water in the western North Atlantic, a small region in which production is clearly non-Redfieldian.

Controls on the partitioning between benthic and water-column denitrification
We find a median value of B/W = 1.7 in our suite of optimized models, with a range of 1.3 < B/W < 2.3. This significantly reduces the uncertainty on B/W over that derived from geochemical box models, which gave a range of 1 < B/W < 4 (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007). The relatively low value of B/W determined in this study contrasts with results expected from a linear isotope mass balance model (Brandes and Devol, 2002), which predicts the following relationship for B/W , where δ ≈ 5‰ is the mean nitrogen isotopic ratio of oceanic nitrate, and a balance between inputs by nitrogen fixation and removal by benthic and water-column denitrification is assumed. If both fix = 0‰ and b = 0‰, Eq. (1) predicts B/W = 4. However, the global inverse model predicts that B/W is about 1.9 in this case. This is because the impact of isotopic fractionation associated with water-column denitrification is diminished by the degree of nitrate consumption in suboxic zones (Deutsch et al., 2004). We find that the modeled B/W ratio can be well predicted (Fig. 6) by applying a simple correction to Eq. (1) to account for the degree of nitrate consumption (f c ) in suboxic zones To apply Eq.
This predictive relationship suggests that if there is complete consumption of nitrate in suboxic zones (f c = 1), then the isotopic enrichment effect of water-column denitrification does not affect the mean ocean δ 15 NO 3 . In fact, the full effect of the isotopic enrichment due to water-column denitrification can never be achieved, because f c > 0 whenever there is any water-column denitrification. In the model, the average f c within suboxic zones of about 0.34 effectively reduces the enrichment factor for water-column denitrification from about 25 ‰ to about 17‰. The modeled f c in suboxic zones agrees well with the observed average f c = 0.31 within those same locations. Equation (2) slightly overpredicts B/W by about 0.15, on average, primarily because Eq. (2) does not take into account the isotopic enrichment during the assimilation of nitrate to form organic matter. To obtain the magnitude of this effect we re-ran the model with up = 0‰ and all other parameters fixed at their values determined by the inversion process. The results show that the fractionation associated with nitrate assimilation lowers the mean ocean δ 15 NO 3 by about 0.5 ‰, on average (minimum of 0.3 ‰ and maximum of 1.2‰ in all the model runs). Taking uptake fractionation into account, by replacing δ in Eq. (2) by the value of δ in the case that up = 0‰, leads to a prediction that slightly underestimates the modeled B/W by about 0.15, on average. Equation (2) is therefore accurate within about ±0.15, depending on the value of up and on the actual value of B/W . Remaining discrepancies between the value of B/W predicted by Eq. (2) and that determined by the model are likely due to the inadequacy of using a single value of f c to account for the spatially heterogeneous effects of water-column denitrification.
These results emphasize the importance of simultaneously achieving good fits to both nitrate deficits (to get f c correct) and nitrogen isotopes (to get δ correct) in order to derive a good estimate of marine denitrification rates. Using only one or the other constraint can produce misleading results. For example, a box model that was tuned to fit mean ocean δ 15 NO 3 but not N * found a value of B/W ∼ 4 (Brandes and Devol, 2002) because it did not take into account the effects of nitrate consumption in suboxic zones, while an ocean cir- culation model that was tuned to fit mean ocean δ 15 NO 3 but not N * found a value of B/W ∼ 0.5 (Somes et al., 2010), because the modeled f c was too large.
To reduce the uncertainty on the exact value of B/W will require more accurate knowledge of the isotopic enrichments factors for nitrogen fixation ( fix ), benthic denitrification ( b ), and water-column denitrification ( w ). Laboratory experiments with denitrifying bacteria have found values for w as high as 30‰ (Barford et al., 1999) or as low as 10-15 ‰ (Kritee et al., 2012). In our model, the optimal value of w depends on the values of fix and b , and ranges from 19-29 ‰ throughout our suite of model configurations. The lower values of w are associated with high values of b and fix . If b were as high as 4-5 ‰, as suggested by some measurements (Alkhatib et al., 2012), then the optimal value of w would approach the lower values suggested by Kritee et al. (2012).
It is interesting to compare our results to those from a conceptually similar study that systematically tuned the parameters of a multibox ocean model to fit observed N * and mean ocean δ 15 NO 3 (Eugster and Gruber, 2012). The two studies arrive at similar overall denitrification rates (107-188 Tg N yr −1 in Eugster and Gruber (2012) compared to 120-240 Tg N yr −1 in this study) and B/W (1.6-2.0 in Eugster and Gruber (2012)   ocean δ 15 NO 3 and the fractional consumption in suboxic zones (Eq. 2), and that details of the ocean circulation are probably a second order influence on B/W . However, the total rate of denitrification is likely to be more sensitive to details of the ocean circulation, which influences the pattern and magnitude of biological production, the extent of suboxic zones, and the spatial distribution of remineralization. Therefore, the agreement between the box model and our study in this regard may be partly fortuitous.

Uncertainty in denitrification rates
The globally integrated rate of marine denitrification predicted by the model ranges from about 120-240 Tg N yr −1 , with a median rate of 170 Tg N yr −1 (Table 1). We find that uncertainty in the isotopic enrichment factors ( fix and b ) and uncertainty in the fraction of organic matter production routed to the DON pool (σ DON ) contribute approximately equally to the uncertainty in the globally integrated denitrification rate (Fig. 7a and c). Denitrification rates generally increase with larger values of b and fix (Fig. 7a). This relationship follows because for larger values of b or fix , a larger B/W ratio is needed to achieve a mean ocean δ 15 NO 3 of ∼5‰ (Fig. 7b), which is achieved in the model by increasing benthic denitrification rates. It is also the case that B/W increases with smaller σ DON values (Fig. 7d). This is because with smaller σ DON values a larger fraction of particulate organic matter sinks out of the euphotic zone, yielding a larger supply of organic matter to the sediments.
The ratio B/W is relatively well constrained despite large uncertainties in the isotopic enrichment factors. This is because in the inverse model, the fractional consumption f c in suboxic zones generally increases slightly with increasing values of b or fix , while w decreases with increasing values of b or fix . Both of these effects reduce the sensitivity of B/W to the isotopic enrichment factors for fixation and benthic denitrification.
The effects of the N * and δ 15 NO 3 constraints on the ratio of B/W predicted by the model can be illustrated by comparing the relative model-data misfit for each observational constraint as a function of B/W . The results for one particular model configuration show that the value of B/W needed to optimally match only the observed N * is about 1.4, while that needed to optimally match only the observed δ 15 NO 3 is about 2.3 (Fig. 8). When both constraints are used, the model strikes a compromise between the two constraints such that the optimal value of B/W is about 1.8 (Fig. 8). In this particular model configuration, N * provides the stronger constraint on B/W , as evidenced by the deeper minimum associated with the relative model-data misfit. This is owing to both the larger number of N * observations compared to δ 15 NO 3 observations, and the fact that N * generally shows a stronger sensitivity to changes in B/W (holding all other model parameters fixed) than does δ 15 NO 3 .
The fact that the N * and δ 15 NO 3 constraints require different optimal B/W values is not surprising, given uncertainties in the model parameters and imperfections inherent in representing complex phenomena with simple parametric equations, as well as uncertainties in the data. Only with perfect data and a perfect model could we expect the model to match both data sets optimally with the same set of parameters. This further illustrates the importance of bringing together both the N * and δ 15 NO 3 data as constraints on denitrification rates, as long as one deals with imperfect models and imperfect data. Furthermore, the advantage of the inverse model is that uncertainties are explicitly coded into the model in terms of the adjustable control parameters, and the model is given freedom to choose between different parameter values in order to optimally match the observed N * and δ 15 NO 3 .

Implications and conclusions
The results of this study suggest that the marine nitrogen budget is unlikely to be strongly out of balance. Previous studies suggesting that marine denitrification rates are much higher than the rate of fixed nitrogen inputs depended on having either a very high rate of water-column denitrification, exceeding 90 Tg N yr −1 , and/or a large value of the ratio of benthic to water-column denitrification (B/W ), exceeding 3 (Codispoti et al., 2001;Codispoti, 2007). Here we have shown that neither of these possibilities is likely, given the constraints provided by observed nitrate deficits (N * ) and nitrogen isotopic ratios of oceanic nitrate (δ 15 NO 3 ).
Results of our global 3-D inverse model simulations suggest that the optimal rate of water-column denitrification needed to match the observed N * and δ 15 NO 3 is about 60 (range of 50-77) Tg N yr −1 , in good agreement with a previous estimate based on N 2 gas measurements (DeVries et al., 2012). Meanwhile, the optimal value of B/W is about 1.7 (range 1.3-2.3). These estimates represent a significant improvement over previous estimates from box models (Brandes and Devol, 2002;Deutsch et al., 2004;Altabet, 2007;Eugster and Gruber, 2012), which could not resolve the 3-D ocean circulation and the full spatial variability in denitrification rates.
While the denitrification rates estimated here significantly reduce the uncertainty on the global rate of fixed N loss from the ocean, some significant uncertainties in the marine N cycle remain that cannot be addressed using the present model. Perhaps most importantly, we have not addressed the magnitude and distribution of N 2 fixation rates. A modeling study that used surface N * distributions to estimate nitrogen fixation rates found a global N 2 fixation rate of about 140 Tg N yr −1 (Deutsch et al., 2007), while a recent observational study suggests that the global rate of N 2 fixation is about 180 Tg N yr −1 (Grosskopf et al., 2012), which is approximately equal to the mean global denitrification rate found in this study. However, it is not known precisely what amount of fixation is supported by the N * and δ 15 NO 3 data in our model. This is because the information contained in the surface N * distribution, which can in principle be used to constrain rates and patterns of nitrogen fixation (e.g., Deutsch et al., 2007), has already been absorbed by the surface restoring condition used to simulate the production of organic phosphorus and organic nitrogen.
Similarly, we have not explicitly considered sources of N due to riverine inputs and atmospheric deposition, which could have significant local impacts on surface N * and δ 15 NO 3 distributions (e.g., Hansell et al., 2007;Knapp et al., 2008). However, whatever spatial pattern these fluxes impart to N * in surface waters is achieved by restoring toward observed NO 3 and PO 4 independently, and their effect on the N reservoir will be implicitly included in our "N 2 fixation" term. In the case of isotopic constraints, the effect of these surface fluxes is not similarly accounted for, but if the isotopic signature of those fluxes is on balance not significantly different from that of N 2 fixation (e.g., Brandes and Devol, 2002) then these inputs could again be considered part of the "N fixation" term that closes the budget. Given the uncertainty surrounding the isotopic ratio of these N inputs, this assumption must be taken as provisional. Lastly, some estimates suggest that human activities may have more than doubled riverine and atmospheric N inputs from their preindustrial values (Gruber and Galloway, 2008). The effect of these anthropogenic perturbations on nutrient distributions and on denitrification rates is poorly known, but should be addressed in future studies.
for warm seawater temperatures (T ) and adequate light (I ) and iron (Fe) supply (e.g., Monteiro et al., 2011). T max is the maximum modeled surface temperature (about 31 • C), which is included to ensure that the maximum of the exponential expression in equation (A6) is 1. Iron is not modeled explicitly, but rather we use a modeled dust deposition field (Mahowald et al., 2006) as a proxy for its availability (e.g., Somes et al., 2010). The surface irradiance I is derived from International Satellite Cloud Climatology Project-C1 data (Zhang et al., 2004). The environmental controls in equation (A6) help the model to reduce N 2 fixation in places where rates are known to be low, such as the Southern Ocean.
Of all the newly fixed organic matter, a fraction φ e is routed to the PON pool and remineralizes in the watercolumn (J wc fix ) and sediments (J sed fix ), following the same formulation as regular organic matter (Eq. A5). Another fraction φ d = φ e × σ DON /(1 − σ DON ) is routed to the DON pool (J DON fix ) and remineralizes following the first-order kinetics for DON remineralization. The remaining fraction (1 − φ e − φ d ) remineralizes immediately in the surface layers where fixation occurs. Thus, the individual fixation terms in Eqs. (A1) and (A2) are Denitrification in the water column occurs wherever local observed O 2 levels fall below a critical level (O 2, crit ) representing the threshold at which denitrification replaces oxic respiration as the dominant pathway of organic matter degradation: crit where r N denit :N org represents the ratio of moles NO 3 used to respire 1 mol of organic nitrogen. The "if and only if" ⇐⇒ statement is handled by creating a mask from observed monthly climatology of oxygen concentrations (Garcia et al., 2010b) after applying the correction suggested by Bianchi et al. (2012) using the procedure described by DeVries et al. (2012). Benthic denitrification is parameterized as a function of the rate of organic matter respiration in the sediments: where F is a function that accounts for enhanced sedimentary denitrification rates under low-oxygen and high-nitrate conditions, where C O 2 , K O 2 and K NO 3 are parameters governing the oxygen and nitrate dependence of sedimentary denitrification. A hyperbolic tangent was chosen for the O 2 dependence of benthic denitrification to allow for uncertainty in how the transition to O 2 -inhibition takes place (i.e., it can be either an abrupt or a smooth transition). Ultimately, we found the parameter controlling the smoothness of the transition (K O 2 ) to be poorly constrained (see Table B1). The linear dependence of benthic denitrification on sediment organic matter flux (Eq. A11) was chosen because of its simplicity and ease of implementation in the model. By contrast, Middelburg et al. (1996) suggest a nonlinear relationship between benthic denitrification and organic matter fluxes. However, we found that our linear formulation produces similar overall rates and spatial patterns to the Middelburg et al. (1996) formulation.
We also tested a version in which benthic denitrification depended quadratically on sedimentary organic matter fluxes, and found no improvement over the linear model. The coupled system of nonlinear equations (Eqs. A1-A2) for NO 3 and DON are solved using Newton's method, which produces convergence to an equilibrium state orders of magnitude faster than traditional time-stepping techniques (Kwon and Primeau, 2006). Fast convergence to an equilibrium solution is necessary for application in the inverse model, which requires O (10 3 ) runs of the forward model to converge to a solution.
The governing equations for 15 NO 3 and DO 15 N are the same as Eqs. (A1) and (A2) except that a fractionation factor α representing the discrimination of chemical reactions toward the lighter isotope is introduced in each term that involves a chemical reaction (e.g., Deutsch et al., 2004). Generically, the reaction rate (J reac ) for 15 N is related to the reaction rate for 14 N by This fractionation effect is taken into account in the uptake of NO 3 to form organic N (J prod ), the remineralization of organic N by denitrifying bacteria in the water column (J wcd ) and the sediments (J sd ), and the fixation of atmospheric N 2 (J fix ). The isotopic enrichment factor for a reaction is given by = (1 − α) × 1000.
Given the steady-state solution for NO 3 and DON obtained from solving Eqs. (A1) and (A2), and making the approximation 14 NO 3 ≡ NO 3 and DO 14 N ≡ DON, the equations for 15 NO 3 and DO 15 N can be cast in terms of a coupled linear 2492 T. DeVries et al.: Marine denitrification rates from a global inverse model system of equations for the isotopic ratios: where R std is the isotopic ratio of atmospheric N 2 , which we take to be 1 for convenience. The resulting set of coupled linear equations can be solved by direct matrix inversion. In places where NO 3 concentrations are very low (close to zero) we find that the ratio R NO 3 can become ill-defined, leading to noise in the model-simulated R NO 3 distribution. For this reason, we neglect R NO 3 values above 200 m depth, where very low NO 3 can occur, when comparing the modeled and observed δ 15 NO 3 values in the inverse model (see below). We find that a similar problem occurs where DON concentrations are close to zero, which occurs in many places in the interior ocean. These points can cause an ill-conditioned (nearly singular) matrix when inverting for the modeled R NO 3 and R DON values. We find that this problem is eliminated when we set all values where DON< γ to γ , where γ is a small number (we used 10 −4 µM). This does not affect the modeled R NO 3 values. When comparing modeled to observed isotopic ratios for NO 3 , we convert observed δ 15 NO 3 values to R NO 3 values using the relationship δ 15 NO 3 = (R NO 3 − 1) × 1000.

Inverse model description
The procedure by which the model is fit to observed N * and δ 15 NO 3 involves two steps. In the first step, the model circulation and air-sea fluxes are adjusted to minimize the misfit between modeled and observed temperature, salinity, radiocarbon, and phosphate distributions. This procedure follows that outlined in DeVries and Primeau (2011), except that here we use a higher resolution model grid (2 • horizontal resolution with 24 unevenly spaced vertical levels) and we include phosphate observations from the 2009 World Ocean Atlas gridded database (Garcia et al., 2010a) in the set of observations constraining the model. The cycling of phosphate and dissolved organic phosphorus (DOP) are both explicitly modeled, following the same sets of equations described for the nitrogen cycle, except that the fixation and denitrification terms are of course not included. We include the depth attenuation coefficient b for particulate organic phosphate (POP) remineralization as an additional control parameter of the model to be determined as part of the optimization. Two additional parameters, σ DOP and τ DOP are needed for the model, but these cannot be determined as part of the optimization because DOP data is not included in the set of observational constraints. Rather, we specify σ DOP and τ DOP based on values determined in previous studies. In one model we specify σ DOP = 2/3 and τ DOP = 1/2 yr (Najjar and Orr, 1998), and in another model we specify σ DOP = 1/2 and τ DOP = 2 yr (Schlitzer, 2002).
The relative error of model state variables in this first step of the optimization is 0.8 for temperature, 0.8 for salinity, 0.65 for 14 C, and 0.75 for phosphate in the case that σ = 1/2 and τ DOP = 2 yr (0.78 in the case that σ DOP = 2/3 and τ DOP = 1/2 yr). Relative errors of about 1 indicate that the model-data residuals are distributed according to the prior estimated error covariance for the global gridded data sets (cf. DeVries and Primeau, 2011). Primarily due to computational restrictions, we did not use CFC observations to constrain the circulation, unlike a previous study (DeVries et al., 2012). When compared to observed CFC concentrations, the model does show some deficiencies in ventilating the suboxic zones. In particular, the Arabian Sea suboxic zone is apparently too weakly ventilated near the surface and too well ventilated at depth, while the eastern tropical South Pacific is too well ventilated throughout. In the future it will be important to assimilate CFC observations into the model to reduce these circulation errors. Despite this, we find that at the global scale the water-column denitrification rates derived here (50-77 Tg N yr −1 ) are similar to those derived by DeVries et al. (2012) (66±6 Tg N yr −1 ) using a model that was constrained by CFCs and that used N 2 /Ar observations to constrain denitrification rates. This indicates that circulation errors within the suboxic zones are small enough not to have a large impact on the globally integrated water-column denitrification rates.
The circulation found in step one of the optimization is then taken offline and used in the nitrogen cycle simulation. Most of the parameters of the nitrogen cycle model (Eqs. A1-A17) are included as control parameters in this second step of the inversion (Table B1). One important exception is that we fix the depth attenuation coefficient for PON remineralization (b) at the value found in step one of the inversion. The value of b determined from step one of the inversion is 0.77 (for the model in which σ DOP = 2/3) or 0.79 (for the model in which σ DOP = 1/2). This is within the range of estimates for the globally averaged value of b based on sediment trap data (e.g., Primeau, 2006), although we do not account for the fact that the value of b within suboxic zones may be lower than the global average value due to reduced respiration rates at low oxygen concentrations (e.g., Van Mooy et al., 2002). The choice of a constant and uniform value of b for both POP and PON remineralization is obviously a simplification, but one that in this case makes the problem more computationally tractable. Despite this simplification, we find that our model well matches the observed N deficits in suboxic zones (Fig. 2d), and that the water-column denitrification rates inferred by our model agree with recent estimates that did take into account the possibility of different b values within suboxic zones (Bianchi et al., 2012;DeVries et al., 2012). In the context of the inverse model, it may be that there are enough additional degrees of freedom (e.g., in Table B1. Prior (pre-optimization) and posterior (post-optimization) values of the control parameters, and associated uncertainties. Prior uncertainties are 1 standard deviation of a normal distribution.  Barford et al. (1999), e Depends on value of σ DON used and varies from 0.4 ± 0.1 yr for σ DON = 2/3, to 6.6 ± 0.8 yr for σ DON = 1/4. the ratio of NO 3 consumed to organic N remineralized, or in the critical O 2 threshold for denitrification) to effectively make up for any biases in the value of b within suboxic zones. Several other variables were held fixed (i.e., not included as control parameters of the inverse model) including the isotopic enrichment factors for nitrogen fixation ( fix ), benthic denitrification ( b ), and assimilation of nitrate to form organic matter ( up ), which are fixed at various values in different model configurations; and the fraction of organic nitrogen routed to the DON pool (σ DON ), which is fixed at either σ DOP or 0.5σ DOP , depending on the model configuration. The values of b , fix , up , and σ DON used in the different model configurations are given in Sect. 2.2.
In total, there are 21 parameters that are iteratively adjusted to find the optimal solution (Table B1). The adjustable parameters include 6 parameters controlling the rate and spatial pattern of N 2 fixation (F o , λ, T o , K I , K Fe , φ e ); 6 parameters controlling the rate and spatial pattern of water-column denitrification (O 2,,crit and r N denit :N org are allowed to vary separately in the Indian Ocean, the South Pacific, and the North Pacific); 7 parameters controlling the rate and spatial pattern of benthic denitrification (a 0 − a 3 , C O 2 , K O 2 , K NO 3 ); the timescale for remineralization of DON (τ DON ); and the enrichment factor for water-column denitrification ( wcd ).
The optimal solution is defined as the set of control parameters that minimize the following cost function: where n N * (n δ ) is the number of grid cells with N * (δ 15 NO 3 ) observations, and n p = 21 is the number of control parameters. p pos represents the final (posterior) value of the model parameters, and p pri their prior values, which were used as the initial guess (see Table B1). We chose σ 2 p to be very large for all parameters so that the model solution was not biased toward our initial guess. In all solutions, the value of the last term in Eq. (B1) is much smaller than the value of the first two terms.
The optimization procedure takes an initial guess at the set of control parameters and iteratively adjusts the parameters to find the minimum of the cost function using a quasi-Newton algorithm. The algorithm typically takes several hundred iterations to find a suitable minimum (defined where the gradient of the cost function with respect to the parameters is less than 10 −3 ). This requires O (10 3 ) simulations, each of which necessitates computing the steady-state solution to the model equations. This large number of simulations is made possible by applying Newton's method to the nonlinear nitrogen cycle equations, which allows rapid convergence to a steady-state solution. Table B1 lists the control parameters of the inverse model, their initial guesses and their final values. The initial guesses for parameter values were set by published estimates where available, and by hand tuning to achieve rough consistency with observations for parameters without a published estimate. The final (post-optimization) parameter values are the mean of the values at the end of each of the 48 different optimizations (using different values of fix , b , and σ DON ), and the uncertainty is the standard deviation of each parameter in this set of 48 optimal solutions. We should note that the quasi-Newton algorithm cannot distinguish between global and local minima. However, most parameters have a relatively small posterior uncertainty indicating a strong minimum (Table B1). Some parameters have a large posterior uncertainty (e.g., K O 2 ), which might reflect multiple local minima in that parameter, but more likely reflects a broad, weak minimum indicating that these parameters are poorly constrained.