Constraints on global oceanic emissions of N 2 O from observations and models

. We estimate the global ocean N 2 O ﬂux to the atmosphere and its conﬁdence interval using a statistical method based on model perturbation simulations and their ﬁt to a database of ∆ pN 2 O (n=6136). We evaluate two submodels of N 2 O production. The ﬁrst submodel splits N 2 O production into oxic and hypoxic pathways following previous publications. The second submodel explicitly represents the redox transformations of N that lead to N 2 O production (nitriﬁcation and hypoxic denitriﬁcation) and N 2 O consumption (suboxic denitriﬁcation), and is presented here for the ﬁrst time. We perturb 5 both submodels by modifying the key parameters of the N 2 O cycling pathways (nitriﬁcation rates, NH +4 uptake, N 2 O yields under oxic, hypoxic and suboxic conditions), and determine a set of optimal model parameters by minimisation of a cost function against 4 databases of N cycle observations derived from observed and model ∆ pN 2 O concentrations. Our estimate of the global oceanic N 2 O ﬂux resulting from this cost function minimisation is 2.4 ± 0.8 and 2.5 ± 0.8 Tg N y − 1 for the 2 N 2 O submodels. These estimates suggest that the currently available observational data of surface ∆ pN 2 O constrain the global 10 N 2 O ﬂux to a narrower range relative to the large range of results presented in the latest IPCC report.


Introduction
Nitrous oxide (N 2 O) is the third-most important contributor to anthropogenic radiative forcing, after carbon dioxide (CO 2 ) and methane (CH 4 ) (Myhre et al., 2013).It is also currently estimated as the dominant contributor to stratospheric ozone depletion (Portmann et al., 2012).Yet our quantita-tive understanding of the magnitude and processes controlling natural N 2 O emissions from the Earth surface to the atmosphere is very poor.A range of methods have been used to constrain total oceanic N 2 O emissions, including the combination of surface ocean N 2 O partial pressure anomalies with gas-exchange parameterisations (Nevison et al., 1995), empirically derived functional relationships applied to global ocean datasets (Nevison et al., 2003;Freing et al., 2012), and ocean biogeochemistry models (Suntharalingam and Sarmiento, 2000;Suntharalingam et al., 2000;Jin and Gruber, 2003;Martinez-Rey et al., 2015).In spite of the multiple methods used, the reported oceanic emissions of N 2 O are still poorly constrained, ranging from 1.9 to 9.4 Tg N yr −1 according to the latest report of the Intergovernmental Panel on Climate Change (IPCC; Ciais et al., 2013).The uncertainty in the oceanic emissions of N 2 O accounts for a large part of the total uncertainty in the natural N 2 O emissions, which is approximately 11 Tg N yr −1 (Ciais et al., 2013).Part of the uncertainty in the oceanic emissions is whether estuaries are included, which could emit as much as 2.3-3.6 Tg N yr −1 (Bange et al., 1996).
The large uncertainty in the oceanic emissions of N 2 O stems from the complexity of its production pathways.There are two main pathways of N 2 O production in the ocean, nitrification and denitrification, which both stem from redox reactions of nitrogen, under oxic and hypoxic conditions, respectively (Fig. 1).N 2 O is formed as a byproduct of marine nitrification of ammonium (NH + 4 ) to nitrate (NO − 3 ); N 2 O is also an intermediate product of denitrification, during the reduction of NO − 3 to nitrogen gas (N 2 ) (Frame and Casciotti, 2010;Loescher et al., 2012;Merbt et al., 2012).Denitrification can also consume N 2 O, using extracellular N 2 O, and reduce it to N 2 (Bange, 2008).In the oxic part of the ocean (i.e.most of the ocean, 97 %>34 µmol O 2 L −1 , using O 2 data taken from Bianchi et al., 2012) denitrification is suppressed, and the primary formation pathway is usually ascribed to nitrification (Cohen and Gordon, 1978), although denitrification may be significant in the anaerobic centres of large marine snow particles in oxic waters (Klawonn et al., 2015).Oceanic N 2 O production in oxic regions is often derived from the linear relationships observed between apparent oxygen utilisation (AOU) and apparent N 2 O production ( N 2 O) (e.g.Yoshinari, 1976;Cohen and Gordon, 1978).However, the N 2 O / AOU ratio varies in different water masses and oceanic regions (Suntharalingam and Sarmiento, 2000).Previous studies have suggested that differences in the N 2 O / AOU ratio could be driven by changing N 2 O yields under varying pressure and temperature (Butler et al., 1989) or varying O 2 concentration (Nevison et al., 2003).Additional mechanisms not yet quantified could include variations in the elemental stoichiometry of the organic matter that is being remineralised, and spatial separation of organic matter remineralisation and nitrification.Throughout the paper we will refer to N 2 O stoichiometries relative to O 2 , NH + 4 and NO − 3 as ratios, because they have been optimised against global databases of concentration measurements, rather than from microbiological yields.Using the latter would be more mechanistically satisfying, but the relevant yields are at present insufficiently constrained by observations.
Estimates of the contribution from suboxic regions of the ocean (about 3 %) to the global N 2 O flux vary from net depletion via denitrification (Cohen and Gordon, 1978) to 33 % for the total N 2 O production in the suboxic ocean (Suntharalingam et al., 2012) and to more than 50 % from denitrification alone (Yoshida et al., 1989).This ambiguity remains unresolved.Bottom-up microbial physiology data are relatively scarce (see Sect. 2.4-2.6), while top-down data need relatively complicated inverse methods to estimate the contribution from suboxic regions.These inverse methods are complicated both because of the variation in the N 2 O / AOU ratio, which is negative under suboxic conditions, maximal under hypoxic conditions and lower under oxic conditions (e.g.0.31-0.033mmol mol −1 , Law and Owens, 1990), and because the influence of mixing gradients makes in situ ratios an unreliable gauge to the biological yields under in situ conditions (Nevison et al., 2003).
Here, we estimate the global ocean N 2 O flux to the atmosphere and its confidence interval.First, we estimate N 2 O flux from observations only (Sect.2.1).This estimate has large uncertainty.We subsequently use a statistical approach introduced by Buitenhuis et al. (2013a) to estimate the global oceanic emissions of N 2 O and its confidence interval by combining ocean N 2 O model simulations with a global database of measurements of surface pN 2 O.This approach involves minimisation of a cost function that compares a series of model simulations with a global database of point measurements of surface pN 2 O.To achieve this, we use four observational databases of the N cycle (Sect.2.2) to extend the global ocean biogeochemistry model Plank-TOM10 (Le Quere et al., 2016b) with additional N cycle processes.We derive the biogeochemical parameters for nitrification rate and phytoplankton use of NH + 4 from the observational databases of nitrification rate and NH + 4 concentration (databases 1 and 2 and Sect.2.4-2.5).Then, we describe two separate submodels of different levels of complexity that represent N 2 O cycling pathways (Sect.2.6-2.7).Finally, we apply the statistical approach (Sect.2.8) to the two submodels to estimate the N 2 O production in the low-O 2 regions from the depth resolved N 2 O concentration database (database 3 and Sect.3.1), and the global oceanic N 2 O flux from the surface pN 2 O database (database 4 and Sect.3.2), followed by a discussion of the results (Sect.4).

Calculation of global ocean N 2 O production from N cycle observations
In this section we provide an initial estimate of global marine N 2 O production based on observationally derived quantities characterising marine productivity and the global ocean N cycle.This follows a similar method to Cohen and Gordon (1979), who estimated ocean N 2 O production using Redfield type ratios.N 2 O is produced either during production of NO − 3 in NH + 4 oxidation or during NO − + ρ urea )) as a function of latitude, from 15 N uptake experiments.Small dots were estimated without measuring NH + 4 or urea concentrations (Prakash et al., 2008(Prakash et al., , 2015;;Gandhi et al., 2010Gandhi et al., , 2012)).Large dots did not give a significant linear relationship with absolute value of latitude and were therefore averaged at 0.29 ± 0.18 (Wafar et al., 2004;Varela et al., 2005Varela et al., , 2013;;Joubert et al., 2011;Thomalla et al., 2011;Simpson et al., 2013).
N 2 O production on total NO − 3 turnover, calculated from primary production times the f ratio.The f ratio is the fraction of primary production that is supported by nitrate.Primary production (PP) was estimated at 58 ± 7 Pg C y −1 based on 14 C primary production measurements (n = 50 050), parameter perturbations of a previous version of the model used here, and Eq. ( 5) (Buitenhuis et al., 2013a).We compiled a database of uptake rates of NO − 3 , NH + 4 and urea, which gives an average f ratio of 0.29 ± 0.18 (Fig. 2, large symbols, n = 34).The globally averaged N 2 O / AOU ratio was calculated from the MEMENTO database (Bange et al., 2009) as 81.5 ± 1.4 µmol / mol (Fig. 3).Finally, since primary production is expressed in carbon terms, and N 2 O production was correlated with oxygen (O 2 ) utilisation, we need to include the −O 2 : C ratio (the − sign indicates the O 2 is consumed as CO 2 is produced), which was taken from Anderson and Sarmiento (1994) as 170 ± 10 / 117 ± 14, and the molar weights of C (12) and N in N 2 O (28).Here and in the rest of the paper, errors were propagated in the usual way: Thus N 2 O production was calculated as PP ×f ratio × −O 2 : C × N 2 O / AOU.Our best estimate of N 2 O production using this method is 58 ± 7 × 1000 × 0.29 ± 0.18 × 170 ± 10/117±14×81.5×10−6 ±1.4×10 −6 ×28/12 = 4.6 ± 3.1 Tg N yr −1 .This estimate lies in the middle of other reported estimates (Fig. 4), but the 68 % confidence interval is very large.We therefore investigate the N 2 O fluxes using a model optimised with observations in the rest of the paper.Nevison et al., 1995;Nevison et al., 2003;Freing et al., 2012Freing et al., (plotted in 2011)); Bianchi et al., 2012;this study -atmospheric inversions: Hirsch et al., 2006;Huang et al., 2008;Thompson et al., 2014Thompson et al., (plotted in 2013)); Saikawa et al., 2014), model estimates shown as crosses (Suntharalingam and Sarmiento, 2000;Jin and Gruber, 2003;Suntharalingam et al., 2012;Martinez-Rey et al., 2015).

Observational databases for model development
We used four databases to tune or optimise different aspects of the N cycle in the PlankTOM10 ocean biogeochemistry model.The number of data points reported for each database are after gridding to 1

Cost function formulation
To parameterise the model N cycle, we use a cost function to minimise the difference between model and observations, following the methods of Buitenhuis et al. (2013a): This formulation gives equal weight to the relative correspondence between model and observations at small and large observational values.A value of 2 means that, on average, the model deviates from the observations by a factor 2 in either direction.To calculate the cost function (and also to calculate MSE in Eq. 6), the model was regridded to the same grid as the observations, and residuals were calculated at months and places where there are observations.The cost function results for the optimised simulations are summarised in Table 1.

Nitrification
Our initial biogeochemical model configuration is Plank-TOM10 (Le Quere et al., 2016b), which represents growth and loss terms from 10 plankton functional types (PFTs), including N 2 fixers, picoheterotrophs (Bacteria plus Archaea) and denitrification rate, but not denitrifier biomass.A full model description and parameter values are provided in the Supplement.Here, we extend the model representation of redox reactions in the N cycle, to create the global biogeochemical model PlankTOM10.2.We describe the new N cycle components below.
In order to represent nitrification rate, the state variable for dissolved inorganic nitrogen was split into NO − 3 and NH + 4 .Respiration by all PFTs produces NH + 4 .The parameterisation for nitrification used in our model is based on the analysis of a database of NH + 4 -specific nitrification rates (Yool et al., 2007).Yool et al. (2007) found that observed nitrification rates are highly variable, with no obvious relationship with either latitude or depth.In their model they therefore used a constant rate of 0.2 d −1 throughout the ocean.Implementing this rate in our model resulted in a cost function relative to the nitrification rate observations of 4.22 (Table 1).We tested if including temperature, O 2 or light dependence improves the ability of the model to reproduce observed nitrification rates.Regarding the response of ammonia-oxidising Archaea (AOA), the main nitrifiers in the ocean (Francis et al., 2005;Wuchter et al., 2006;Loescher et al., 2012), to temperature, we are only aware of the measurements of Qin et al. (2014).These show a ∼ 4-fold variation in maximum growth rate between three strains, which poorly constrains the temperature dependence of AOA.We therefore first used a generic Q 10 of 2 and optimised the rate at 0 • C using the nitrification rate observations.This led to only a slightly improved representation of the observations (cost function = 4.18).Although the response of AOA and ammonia-oxidising Bacteria (AOB) to O 2 has only been measured at 21-25 • C (Frame and Casciotti, 2010;Loescher et al., 2012), which limits the range of O 2 concentrations, there was a significant logarithmic relationship between N 2 O yield and O 2 (Fig. 5).A logarithmic function fitted the data better than linear, exponential or power functions.Since nitrification consumes O 2 , in the model it decreases as remineralisation switches from O 2 to NO 3 (Eqs.61, 67 and 70 in the Supplement).Implementing this response to O 2 led to only a further small improvement of the model nitrification rate relative to the observations (cost = 4.16).This implies that nitrification never becomes O 2 limited, reflecting a lack of data to parameterise an expected decrease.As will be described more fully in Sect.3.1, we used observed O 2 concentrations in the simulations (Bianchi et al., 2012) rather than interactively modelled O 2 to minimise the impact of model biases in simulated O 2 fields (Suntharalingam et al., 2012).The response of AOA to light is estimated to be 50 % inhibited at 5 µmol photons m −2 s −1 .However, this estimate is not well constrained (Merbt et al., 2012).Implementing this light response did not improve the model, either in combination with the O 2 and temperature responses or with the temperature response only, and was subsequently omitted.The lack of improvement in nitrification rates by adding light inhibi-tion might reflect the lower sensitivity of AOA to light found by Qin et al. (2014).

Phytoplankton K 1/2 for NH + 4 uptake
We used the calculation of the preferential uptake of NH + 4 over NO − 3 by phytoplankton PFTs of Vallina and Le Quere (2008) (Eq. 9 in the Supplement).The K 1/2 of phytoplankton for NH + 4 has mostly been measured based on uptake rates (syntheses by Goldman and Glibert, 1983;Killberg-Thoreson et al., 2014).Aksnes and Egge (1991) have shown a theoretical expectation of a linear increase of K 1/2 with cell radius.The observations are so variable that they neither confirm nor contradict such an increase.The model uses a fixed C : N : O 2 ratio for all organic matter of 122 : 16 : −172, and Michaelis-Menten kinetics for growth based on inorganic N uptake by phytoplankton (Buitenhuis et al., 2013a; Eqs. 8 and 9 in the Supplement).We therefore need a K 1/2 for growth rather than for uptake to be consistent with the fixed C : N ratio (Morel, 1987).The available uptake rate data do not include the supporting data to allow conversion to the K 1/2 for growth.We are only aware of measurements of the K 1/2 for growth by Stawiarski (2014).Based on the latter values of 0.09 ± 0.15 µmol L −1 for picoeukaryotes, the K 1/2 of phytoplankton for NH + 4 was set to 0.1 to 5 µmol L −1 , increasing linearly with nominal size (Buitenhuis et al., 2013b).Due to the highly dynamic nature of NH + 4 turnover, the model produces a much smoother distribution of NH + 4 concentrations than the observations, but the large-scale pattern of surface NH + 4 concentration shows an increase with latitude, consistent with the observations (Fig. 6), which translates into a cost function of 3.0.

N 2 O production
N 2 O production is implemented as two distinct submodels.The diagnostic submodel is based on statistical relationships of N 2 O / AOU ratios taken from observations and has previously been published (Suntharalingam et al., 2000(Suntharalingam et al., , 2012)).In oxic waters it uses one ratio to estimate the open ocean source of N 2 O production.In hypoxic waters it uses a higher ratio to represent the increased yield of N 2 O from both nitrification and denitrification in oxygen minimum zones.The hypoxic N 2 O yield is maximal at 1 µmol O 2 L −1 , and decreases with an e-folding concentration of 10 µmol O 2 L −1 (Suntharalingam et al., 2000(Suntharalingam et al., , 2012;;Eqs. 35, 67 and 69 in the Supplement).Previous studies using regional databases have found different oxic ratios (Suntharalingam and Sarmiento, 2000, and references therein).Therefore, both the oxic and hypoxic ratios have been reoptimised to the global databases (Sect.3.1-3.2).
The prognostic submodel presented here is based on process understanding and explicitly represents the primary N 2 O formation and consumption pathways associated with the marine nitrogen cycle (Fig. 1).It includes the production of  (Frame and Casciotti, 2010;Loescher et al., 2012); crosses: marine AOB at high cell numbers (Goreau et al., 1980;Frame and Casciotti, 2010); plusses: soil AOB at high cell numbers (Lipschultz et al., 1981).Black line: logarithmic fit to AOA and low to medium cell number AOB (yield N 2 O during oxic nitrification (blue arrows in Fig. 1) and during hypoxic denitrification (red arrow in Fig. 1), as well as a consumption term during denitrification at even lower (suboxic) O 2 concentrations (yellow arrow in Fig. 1).The ratios of the three processes are globally invariant (Eqs.61, 63, 70 and 71 in the Supplement).The functional form of the O 2 dependence of N 2 O consumption (Eq.71 in the Supplement) was the same as that of denitrification (Eq.67 in the Supplement), and with an O 2 response function that is 1.5 µmol L −1 lower than that of denitrification, which is similar to that used by Babbin et al. (2015).We independently optimised the ratios of N 2 O production and consumption from denitrification (Sect.3.1), which controls the net N 2 O production as a function of O 2 concentration.There is not enough information at present to optimise the O 2 concentration parameters of denitrification and N 2 O consumption as well.The low O 2 ratios of both submodels (Supplement Sect.8.7) were optimised using the database of observed N 2 O concentration (Sect.3.1) and the oxic ratios of both submodels were optimised using the database of observed pN 2 O (Sect.3.2).The N 2 O concentrations from both the diagnostic and the prognostic submodels are transported in the same way by physical transport, and the formulation of their gas exchange is also identical.
in which K0 is the solubility (Weiss and Price, 1980), p watervapour is the water vapour pressure (Sarmiento et al., 1992), piston velocity = 0.27 × (windspeed) 2 (Sweeney et al., 2007), which is optimised for use with the NCEP reanalysis data used here, the Schmidt number for N 2 O was taken from Wanninkhof (1992), and the ice cover is calculated by the sea ice model LIM2.
In most of the simulations, atmospheric pN 2 O was calculated from Eq. (2).For the optimised low O 2 production we also ran a series of simulations with the NOAA pN 2 O atm observational data that included seasonal and latitudinal variations (see Sect. 2.2 for the ftp address where we downloaded the data, and Sect.3.2 for the results).Between 2000 and 2014, we used the monthly observations for the 12 available latitudes.Monthly anomalies relative to the global average were calculated at each available latitude from the 2000-2016 observations.These were added to Eq. ( 2) from 1965 and 1976, and to the global average observations between 1977 and 1999.In the model simulation, the data were linearly interpolated between the 12 latitudes and monthly observations.
The PlankTOM10.2 biogeochemical model coupled with the two N 2 O submodels is incorporated into the ocean general circulation model NEMO v3.1 (Madec, 2008).The model resolution is 2 • in longitude, on average 1.1 • in latitude and has 30 vertical layers, from 10 m in the top 100 m to 500 m at 5000 m.The model simulations were initialised in 1965 from observations (Le Quere et al., 2016b), with NH + 4 initialised as 0, and N 2 O initialised from a horizontal interpolation of the MEMENTO observations (see Sect. 2.2).Simulations were run to 2014, forced with daily atmospheric conditions from the NCEP reanalysis (Kalnay et al., 1996; for details see Buitenhuis et al., 2013a).Results are reported averaged over the last 5 years.

Estimation of global N 2 O flux from point measurements of pN 2 O
In previous versions of the PlankTOM model (Buitenhuis et al., 2006(Buitenhuis et al., , 2010(Buitenhuis et al., , 2013a) )  found to be more appropriate when the observations span several orders of magnitude.Unfortunately, statistical confidence intervals have only been defined for χ 2 statistics such as Eqs.( 5) and ( 6), which minimise absolute error, so that we end up with two cost functions (Eqs.3 and 5), depending on the application.To estimate the global air-sea flux of N 2 O that best fits the pN 2 O data, and its ±1 − sigma (68 %) confidence interval, we use the formula described in Buitenhuis et al. (2013a): in which MSE is mean square error: MSE min is the MSE of the model simulation that is closest to the observations, and n is the number of gridded observations.
In addition to the uncertainty that arises from the modelobservations mismatch, uncertainty is contributed by the uncertainties in the N 2 O solubility and the piston velocity, the two quantities that connect the measured pN 2 O to the estimated air-sea flux.The uncertainty in the solubility has been estimated as 3 % (Cohen and Gordon, 1978).The uncertainty in the piston velocity has been estimated at 32 % (Sweeney et al., 2007).Uncertainties in the solubility and piston velocity are proportional to uncertainty in the optimised N 2 O airsea exchange because the optimised N 2 O production needs to change proportionally with solubility and piston velocity to achieve the same pN 2 O.

N 2 O production at low O 2
The global N 2 O production rate in oxygen minimum zones (OMZs) was optimised using the depth-resolved N 2 O data of the MEMENTO database.As noted in previous model studies of ocean O 2 , global models do not well represent the extent and intensity of OMZ regions (Bopp et al., 2013;Cocco et al., 2013).The modelled OMZs in PlankTOM10 occur at greater depths than observed, resulting in unrealistic vertical distributions of N 2 O (results not shown).Therefore, following Suntharalingam et al. (2012), the model was run using fixed observed O 2 concentrations (Bianchi et al., 2012), which corrected, in part, the vertical distribution of N 2 O production from the two submodels, though it still occurred at too great depths (Fig. 7).In the equatorial regions and in the Pacific Ocean the N 2 O concentrations are underestimated between ∼ 200 and ∼ 1500 m depth, and overestimated below that.This shortcoming is not significantly improved in the prognostic model (Fig. 7), even though the prognostic model is more detailed, separately representing the processes of N 2 O production and consumption at low O 2 concentrations.The depth of maximum N 2 O in the model is generally deeper than observed, suggesting that organic matter remineralisation may be too low at shallow depths.This is confirmed by the depth profile of NO − 3 , which is underestimated relative to the WOA2009 observations between 100 and 1500 m, and overestimated at greater depths (Fig. 8).In both submodels, the N 2 O concentrations in the deep sea are also too high, but since only 5 % of N 2 O production occurs below 1600 m this does not have a big impact on the global N 2 O fluxes.The addition of N 2 O consumption in the prognostic N 2 O model does result in improvement of the N 2 O depth profiles in the Indian Ocean.
In order to find the optimal N 2 O production that minimises the MSE (Eq.5), we ran a range of simulations in which the low-O 2 N 2 O production was varied in the diagnostic model (Fig. 9a), and a range of simulations in which both the hypoxic N 2 O production and the suboxic N 2 O consumption were varied in the prognostic model (Fig. 9b).The optimum solution for the prognostic model was found at a gross production of 0.33 Tg N yr −1 .The optimised (net) N 2 O production in low-O 2 regions and its confidence interval were 0.16 ± 0.13 Tg N yr −1 for the diagnostic model, and 0.12 ± 0.07 Tg N yr −1 for the prognostic model.In the optimised diagnostic model the hypoxic N 2 O ratio (i.e.net production) is 1.7 mmol N 2 O (mol O 2 ) −1 .In the optimised prognostic model the maximum N 2 O production ratio (i.e.gross production from hypoxic denitrification) is 15.4 mmol N 2 O (mol NO − 3 ) −1 , decreasing to 0 above 34 µmol O 2 L −1 .The maximum N 2 O consumption ratio (from suboxic denitrification) is 15 mmol N 2 O (mol NO − 3 ) −1 , decreasing to 0 above 28 µmol O 2 L −1 .This leads to net production that is always positive and has a maximal ratio of 183 µmol N 2 O (mol NO − 3 ) −1 at 10 µmol O 2 L −1 .

N 2 O flux
We used the surface pN 2 O distribution to constrain the total global N 2 O flux, and the uncertainty arising from the model-data mismatch (the uncertainties arising from solubility and piston velocity are added at the end).We ran a range of simulations in which both the (net) low O 2 and the oxic N 2 O production rates were optimised in both submodels (Fig. 10).pN 2 O provided a better constraint than the N 2 O concentration distribution, since more N 2 O production mostly leads to more N 2 O outgassing to the atmosphere rather than a significant increase in shallow N 2 O concentrations (data not shown).This is because outgassing is proportional to pN 2 O, but N 2 O concentration is proportional to pN 2 O, and pN 2 O / pN 2 O is small in most of the surface ocean.The zonal average surface pN 2 O distribution was well simulated by both submodels (Fig. 11d), and the model ensemble covered a wide range of global N 2 O fluxes (Fig. 10).The total N 2 O flux that best reproduced the pN 2 O distribution was 2.4 ± 0.3 Tg N yr −1 for the diagnostic submodel and 2.5 ± 0.3 Tg N yr −1 for the prognostic submodel (Fig. 10).In the diagnostic model, the optimised oxic N 2 O / AOU ratio was 10.6 µmol N 2 O (mol O 2 ) −1 .In the prognostic model, the optimised oxic nitrification ratio was 123 µmol N 2 O (mol NH + 4 ) −1 .The results were the same in both diagnostic and prognostic submodels for the 2000-2004 and 2005-2009 averages, showing that the model was sufficiently spun up.
High N 2 O fluxes have been reported for the coastal ocean (Bange et al., 1996), and near-shore upwelling regions (e.g.Arevalo-Martinez et al., 2015).To test whether these regions contribute more to the global N 2 O flux than their surface area would suggest, we did the optimisation separately for the coastal ocean (≤ 200 m bottom depth) for the near-shore noncoastal ocean (≤ 2 • from land, > 200 m bottom depth) for the eastern tropical Pacific (180-70 • W, 5 • S-5 • N, > 2 • from land) and the rest of the open ocean (Table 2).The results show that the coastal ocean contributes only 2 % of the global N 2 O flux, less than would be expected from its surface area, although there are also fewer observations on the coast (2 % of the total), so that the relative error is slightly higher.The near-shore non-coastal ocean contributes 14 % of the global N 2 O flux both submodels, hardly more than its areal percentage (13 %), and it is also fairly well sampled (12 % of the observations).The eastern equatorial Pacific Ocean contributes 27 % in the diagnostic submodel and 25 % in the prognostic model, more than its areal percentage (22 %), and it is undersampled (17 %).The open ocean contributes 57-59 %, slightly less than its areal percentage (61 %).This is as expected, because we have separated out the main N 2 O hotspots, but the differences are quite small.When we used observed atmospheric pN 2 O that varied with latitude and month (see Sect. 2.2) the results were essentially the same, with an N 2 O flux of 2.4 ± 0.3 Tg N yr −1 for the diagnostic submodel and 2.6 ± 0.3 Tg N yr −1 for the prognostic submodel (data not shown).
Finally, we add the uncertainties in the solubility and the piston velocity to the total N 2 O flux through error propagation.This gives a total uncertainty of 2.4 ± 0.8 Tg N yr −1 for the diagnostic submodel and 2.5 ± 0.8 Tg N yr −1 for the prognostic submodel.

Discussion
Cohen and Gordon (1979) estimated global N 2 O production directly from N cycle observations as 4-10 Tg N yr −1 .However, they did not have information on the f ratio, so their estimate was based on total N assimilation in primary production.We use an updated estimate of primary production and its error (Buitenhuis et al., 2013a), and compile a database of the f ratio (Fig. 2).We also use a much larger database of the N 2 O / AOU ratio (Fig. 3).We recalculate the N-cycle-based N 2 O production based on these extended databases.We find that we can estimate all the relevant steps in the N cycle with observational data, including their uncertainty (Sect.2.1).At present this uncertainty is still fairly large, at 4.6 ± 3.1 Tg N yr −1 .The uncertainty in this estimate is similar to that in Cohen and Gordon (1979), but our uncertainty is based on the uncertainty in all components of the calculation, while their uncertainty was based only on the uncertainty in the N 2 O / AOU ratio.The upper 60 % of our estimate overlaps with the lower 62 % of the Cohen and Gordon (1979) estimate.The biggest contributor to our uncertainty is the f ratio, especially in the tropics, which constitute 44 % of the ocean surface area, and additional measurements and/or data synthesis could help constrain the N 2 O budget.The f -ratio data are only based on uptake of NO − 3 , NH + 4 and urea, whereas phytoplankton can also take up NO − 2 and organic N (other than urea).One of the major sources of uncertainty in using the N 2 O / AOU ratio is that it is conceptually based on the N 2 O production during nitrification, which uses O 2 as the electron acceptor.N 2 O production during denitrification is spatially separated from the associated O 2 use that is needed to nitrify the NH + 4 to NO − 3 , the electron acceptor in denitrification.This NO − 3 is produced by nitrification, so in terms of mass balance our calculation is still valid, but this N 2 O production would show up as a vertical increase in N 2 O without associated increase in AOU at low O 2 concentrations (high AOU) in Fig. 4.This estimate of global marine N 2 O production derived from analysing the N cycle (4.6± 3.1 Tg N yr −1 ) is statistically indistinguishable from the N 2 O flux derived from pN 2 O observations (2.4-2.5 ± 0.8 Tg N yr −1 ) but has a much larger error.However, further observational constraints could not only reduce the error but also extend our understanding of the whole N cycle, including the option of evaluating the model representation of these N cycle processes against observations, and not just the part that N 2 O plays in them.Such further constraints are also likely to provide the most productive way to reduce unexplained variability that is found in the observations but not in the present models.For example, we have shown that both the N 2 O and NO 3 are underestimated at ∼ 300-1500 m depth and overestimated below ∼ 2000 m (Figs. 6 and 7).Thus, improved representation of mesopelagic remineralisation might lead in improved representation of the N 2 O depth distribution.However, this falls outside the scope of this study.
Models of the global marine C cycle have been in use for decades, and a lot of the available information has been synthesised, cross-correlated and interpreted in detail (Le Quere et al., 2016a;Buitenhuis et al., 2013b).While actual measurements of N utilisation and transformation have also been made in abundance (Figs. 2,3,4,5a,6,7 and 9a), the synthesis and global modelling of these data are less advanced.In addition, N occurs in many different oxidation states in the marine environment (e.g.organic matter and NH + 4 as −3, N 2 as 0, N 2 O as 0 and +2, NO − 2 as +3, and NO − 3 as +5).Therefore, redox reactions complicate the representation of the N cycle a good deal.This lack of data synthesis and of identification of the most important controls in a complex system is reflected in a relatively low ability of the model to model observed nitrification rates and to a lesser extent NH + 4 concentrations (Table 1).
This lack of knowledge also means that partitioning the global marine N 2 O production over the nitrification and denitrification pathways is poorly constrained.Both the diagnostic and the prognostic models assign a small percentage of the total N 2 O production to the denitrification pathway, 6 and 4 % respectively.However, because of the large bias between the observed and modelled N 2 O concentration depth profiles (Fig. 7) these may be underestimates (Suntharalingam et al., 2012;Arevalo-Martinez et al., 2015).Possibly because of the model bias (Figs. 7 and 8), the addition of N 2 O consumption in the prognostic submodel does not lead to a significantly better distribution of N 2 O across depth or between different basins (Fig. 8).As a result, the pN 2 O distributions are also quite similar (Figs. 11 and 12) and the optimised N 2 O flux and confidence intervals of the two submodels are also quite similar (Fig. 10).However, it should also be noted, first, that the optimisation using surface pN 2 O agrees with the optimisation using N 2 O concentration that the contribution of the low-O 2 N 2 O production needs to be low (Fig. 10).Second, the error contribution from the model vs. observed pN 2 O comparison is low, with confidence intervals of 0.3 Tg N yr −1 for both submodels.Third, pN 2 O is equally well modelled above the low-O 2 regions as in the rest of the ocean (Figs. 11 and 12), and the contributions of the coastal and near-shore non-coastal ocean are nearly proportional to their surface areas (Table 2).These three features are supporting evidence for our results that suggest that  Despite these shortcomings, the global marine N 2 O flux is well constrained to 2.4-2.5 ± 0.8 Tg N yr −1 by both submodels (Fig. 10).This constraint reflects the fact that the integrated effect of the different physical and biogeochemical processes determines the surface pN 2 O distribution (Fig. 11), so that the integrated total can be well constrained even if the individual processes are not.The N 2 O flux is at the lower end of previous estimates, and with a similar confidence interval to other recent estimates (Fig. 4).The confidence interval is dominated by uncertainty in the piston velocity (32 %) rather than model-observation mismatches (12 %).Because of differences in methodology it is not possible to provide reasons for why our estimate is lower than the more recent estimates.We can, however, compare our estimate to that of Nevison et al. (1995), because it is also based on a database of pN 2 O. Compared  to their high-end estimate using the piston velocity of Wanninkhof of 5.2 ± 3.6 Tg N yr −1 , our estimate is lower because we use the more recent 13 % lower estimate of piston velocity of Sweeney et al. (2007), and because our pN 2 O of 7.6 ± 18.1 ppb is 25-28 % lower compared to 10.55 natm in Nevison et al. (1995) (the range is calculated based on the water vapour correction for conversion between ppb and natm, which increases from 0.6 to 4.1 % at temperatures from 0 to 30 • C, which brings the values slightly closer together).We also tested how much influence sampling biases of very high supersaturation values might have on the estimated air-sea exchange.If the 40 pN 2 O measurements in the gridded database that are higher than 100 ppb (Fig. 12) are doubled, the optimised N 2 O air-sea exchange becomes 2.8 ± 0.5 Tg N yr −1 for the diagnostic model and 3.1 ± 0.5 Tg N yr −1 for the prognostic model.If the 24 pN 2 O measurements in the gridded database that are higher than 152 ppm are excluded, to decrease the frequency of the highly oversaturated observations down to what both submodels simulate (Fig. 12), the optimised N 2 O flux becomes 2.0 ± 0.2 for the diagnostic model and 2.3 ± 0.2 Tg N yr −1 for the prognostic model.These results still fall within the confidence intervals of the results using the complete database.
Possible biases in ocean physical transport could in theory affect N 2 O production in low-O 2 regions.The indirect impact of ocean physics on low N 2 O production through its impact on the distribution of O 2 , which Zamora and Oschlies (2014) have shown to be substantial, is not quantified here because we used observed O 2 (Bianchi et al., 2012) instead of modelled O 2 .Our model results suggest that the model representation of ocean physics is adequate for the purpose of estimating N 2 O flux from biogeochemical model perturbations.On the one hand, if the model had too much ventilation in the OMZs, shallow N 2 O concentrations would be underestimated, as they are in the model (Fig. 7), but this would also lead to pN 2 O overestimation in the surface areas above the OMZs, which is not the case.The high pN 2 O values are generally lower but spread over a larger area than in the observations (Fig. 11), with a good frequency distri- bution of high pN 2 O (Fig. 12).On the other hand, if the model had too little ventilation in the OMZs, the optimisation would reduce N 2 O production in the OMZs in compensation, but the optimisation to pN 2 O would then estimate a higher OMZ N 2 O production than the optimisation to the N 2 O depth profiles to compensate for the low transport, and this is also not the case.Therefore we conclude that potential biases in ocean physical transport do not appear to have a large direct impact on low N 2 O production.
Global oceanic N 2 O emissions estimated using atmospheric inversion methods based on atmospheric N 2 O concentrations tend to be higher than our results (Fig. 4).However, N 2 O emissions from inversions in the Southern Ocean are lower than the priors (Hirsch et al., 2006;Huang et al., 2008;Thompson et al., 2014;Saikawa et al., 2014).These low Southern Ocean emissions (0.02-0.72 Tg N yr −1 ) are consistent with our results (0.68-0.79 Tg N yr −1 ).South of 30 • S, 88 % of the Earth surface is ocean, resulting in a clearer attribution in the inversions of the atmospheric N 2 O anomalies to ocean fluxes.We suggest that the higher emissions estimates from inversions for the global ocean could be due to a combination of overestimated priors of ocean fluxes in combination with insufficient observational constraints at latitudes north of 30 • S to allow correct partitioning between land and ocean fluxes.Results presented here are for the open and coastal ocean.The largest coastal seas are resolved in our model, although specific coastal processes, such as the interactions with sediments and tides, are not.Our results do not include emissions from estuaries.Fluxes from these could be as large as 2.3-3.6 Tg N yr −1 according to one estimate (Bange et al., 1996), and could be another contributing factor to the difference between our results and those of atmospheric inversions.
To improve the estimate of the ocean N 2 O flux, first, the uncertainty in the piston velocity would need to be reduced.Once that is achieved, further improvements might be possible by a more accurate model representations of the remineralisation length scale and of the physiology of N 2 O producing picoheterotrophs (nitrifying and denitrifying Archaea and Bacteria).
Figure 1.Primary biological pathways of the oceanic nitrogen N cycle represented in the model simulations, along with redox states of N. Nitrification occurs in the oxic ocean (blue arrow).Denitrification yields net N 2 O production in hypoxic conditions (red arrow) and net N 2 O consumption in suboxic conditions (yellow arrow).Only organic nitrogen (N org ), NH + 4 , NO − 3 and N 2 O are represented as model state variables.

Figure 6 .
Figure 6.Surface NH + 4 concentration (µmol L −1 ).(a) Observations (symbol size is 5 • × 5 • ).(b) Model results are for the same months where there are observations, and annual averages everywhere else.(c) Zonal average.Black: observations; red: model results.Model results are for the same months and longitudes as the observations.Latitude y axis to the left of panel (a).

Figure 8 .
Figure 8. Depth (m) profile of average NO − 3 concentration (µmol L −1 ).Black line: WOA2009 synthesis of observations, not interpolated.Red line: model results sampled at the places where there are observations.

Figure 9 .
Figure 9. MSE 0.5 for the two N 2 O submodels compared to the N 2 O concentration database as a function N 2 O production in the low-O 2 regions.MSE min was obtained as the minimum of a second-order polynomial fit (black lines).The 1σ confidence interval, where MSE equals the value calculated from Eq. (5), is indicated by the horizontal lines.(a) Diagnostic submodel; each point represents a simulation with a different low O 2 ratio.(b) Prognostic model; "no c" is with no N 2 O consumption, i.e.net production = gross production.All other lines have a constant gross production (see legend for a description of the symbols, Tg N yr −1 ), and net production varies with different N 2 O consumption rates.Range of parameter values is given in Supplement Sect.8.7.

Figure 10 .
Figure 10.MSE 0.5 for the two N 2 O submodels compared to the pN 2 O database as a function of global N 2 O flux at different (net) N 2 O production rates in the low-O 2 regions.MSE min and confidence intervals as in Fig. 9. (a) Diagnostic submodel.The four lines represent the four best low O 2 production rates from Fig. 9a; each point represents a simulation, different symbols indicate different low O 2 ratios, and points with the same symbols have different oxic N 2 O production ratios.(b) Prognostic submodel.The four lines represent the optimised net production rates at the four best gross production rates from Fig. 9b; points with the same symbols have different N 2 O ratios for nitrification.

Figure 11 .
Figure 11.Surface pN 2 O (ppb).(a) Observations (symbol size is 5 • × 5 • ).(b) Optimised diagnostic model.(c) Optimised prognostic model.Model results are for the same months where there are observations, and annual averages everywhere else.(d) Zonal average.Black line: observations; green dashed: diagnostic model; red dotted: prognostic model.Model results are for the same months and longitudes as the observations.Latitude y axis to the left of panel (a).

Figure 12 .
Figure 12.Frequency distribution of pN 2 O in the observations (solid black), and the optimised simulations of the diagnostic submodel (green squares) and the prognostic submodel (red lines).

Table 1 .
Cost function (Eq. 3) for the optimisation simulations of Sect.2.2-2.4,relative to the respective observational databases.The nitrification rate in bold was used in this study.