MOPS-1.0: towards a model for the regulation of the global oceanic nitrogen budget by marine biogeochemical processes

<jats:p>Abstract. Global models of the oceanic nitrogen cycle are subject to many uncertainties regarding the representation of the relevant biogeochemical processes and of the feedbacks between nitrogen sources and sinks that determine space- and timescales on which the global nitrogen budget is regulated. We investigate these aspects using a global model of ocean biogeochemistry that explicitly considers phosphorus and nitrogen, including pelagic denitrification and nitrogen fixation as sink and source terms of fixed nitrogen, respectively. The model explores different parameterizations of organic matter sinking speed, oxidant affinity of oxic and suboxic remineralization, and regulation of nitrogen fixation by temperature and different stoichiometric ratios. Examination of the initial transient behavior of different model setups initialized from observed biogeochemical tracer distributions reveal changes in simulated nitrogen inventories and fluxes particularly during the first centuries. Millennial timescales have to be resolved in order to bring all biogeochemical and physical processes into a dynamically consistent steady state. Analysis of global properties suggests that not only particularly particle sinking speed but also the parameterization of denitrification determine the extent of oxygen minimum zones, global nitrogen fluxes, and hence the oceanic nitrogen inventory. However, the ways and directions in which different parameterizations of particle sinking, nitrogen fixation, and denitrification affect the global diagnostics are different suggesting that these may, in principle, be constrained independently from each other. Analysis of the model misfit with respect to observed biogeochemical tracer distributions and fluxes suggests a particle flux profile close to the one suggested by Martin et al. (1987). Simulated pelagic denitrification best agrees with the lower values between 59 and 84 Tg N yr−1 recently estimated by other authors.                     </jats:p>

section_content: The balance of fixed, i. e., biotically available, nitrogen in the global ocean is determined by processes that either remove it (denitrification, anammox, burial) from or add it (nitrogen fixation, atmospheric, riverine supply) to the ocean. The magnitude of these biotic and abiotic fluxes, and therefore their combined effects on the oceanic fixed nitrogen inventory, is currently not well constrained. A decade ago, some geochemical and model-based studies suggested rather high fluxes (Codispoti et al., 2001; Gruber, 2004), sometimes with a high imbalance between sources and sinks (Codispoti, 2007), but more recent studies point towards lower, and balanced fluxes (Eugster and Gruber, 2012; DeVries et al., 2013; Somes et al., 2013), in accordance with earlier geochemical estimates (Gruber and Sarmiento, 1997). 

Water-column denitrification and anammox are restricted to suboxic zones, i. e., regions with low oxygen, most notable the Arabian Sea, the eastern tropical North Pacific (ETNP) and the eastern tropical South Pacific (ETSP). Because nitrogen fixers experience their optimum growth in warm (surface) waters they are not generally expected to thrive in waters colder than 18 • C (Breitbarth and LaRoche, 2005; Breitbarth et al., 2007), and are therefore thought to be limited to low latitudes. So far, it is not entirely clear, whether areas of denitrification and nitrogen fixation are tightly coupled in space (Deutsch et al., 2007) or not (Landolfi et al., 2013). Some spatial decoupling could be deduced from the distribution of reported direct measurements of nitrogen fixation (Luo et al., 2012) and also from the temperature limitation of nitrogen fixers and generally cold surface waters associated with eastern boundary upwelling regions overlying regions of nitrogen loss. In case of spatial segregation of nitrogen loss processes and nitrogen fixation, not only the processes them-Published by Copernicus Publications on behalf of the European Geosciences Union. 

selves (and their representation in models) but also the feedback processes and the physical transport processes linking the respective regions are of importance. Relatively slow biogeochemical turnover rates may also be the reason for rather long residence times of nitrogen in the ocean, which are estimated to range between ≈ 1000-4000 years (see Eugster and Gruber, 2012, and citations therein). 

Attempts to further constrain the residence time of marine nitrogen will benefit from a better understanding of both the spatial relation of nitrogen loss processes and nitrogen fixation, and the individual biogeochemical processes themselves. For example, nitrogen loss processes are associated with, and sensitive to, low (and therefore difficult to measure) concentrations of organic substrates and oxidants. Direct incubation measurements have led to different interpretations of the relevance and magnitude of various processes that determine the loss of fixed nitrogen in different ocean regions (e. g., Kuypers et al., 2005; Ward et al., 2009; Bulow et al., 2010). A possible explanation for the apparent discrepancies is the dependence on substrate availability, which is difficult to conserve in incubation experiments (Ward et al., 2008; Galan et al., 2009; Kalvelage et al., 2013). Geochemical estimates based on nutrient ratios and/or the distribution of nitrogen isotopes, on the other hand, integrate over space and time, and therefore depend on our knowledge and assumptions of underlying physics. 

Model-based studies, especially when combined with observations, may provide some insight into the associated processes, and help to integrate over space and time. As noted above, one of the two main drivers in setting the global budget and distribution of nitrogen is denitrification, a process confined to suboxic zones. Unfortunately, many models suffer from systematic deficiencies in the spatial representation of low oxygen areas, with often too large and too intense oxygen minimum zones (OMZs). Possible reasons include deficiencies in the description of diapycnal mixing (Duteil and Oschlies, 2011) and an insufficient representation of the equatorial current system and equatorial deep jets (Dietze and Loeptien, 2013). Applying strongly increased zonal isopycnal diffusivities in the equatorial current band, Getzlaff and Dietze (2013) could improve the performance of coarseresolution models with respect to oxygen and temperature in the eastern equatorial Pacific (EEP). Duteil et al. (2014) recently showed that a very fine (1/10 • ) spatial resolution can significantly improve the representation of the eastern tropical Atlantic OMZs by lateral ventilation via the betterresolved equatorial current system. 

Another possible cause of model deficiencies is the representation of the sensitivity of nitrogen loss processes to ambient oxygen as a rather abrupt switch that does not seem to match recent observations made in the Peruvian upwelling zone (Kalvelage et al., 2011). However, global models that parameterize combined cycles of oxygen and nitrate often represent oxidant dependency in a rather coarse manner, as a hyperbolic or step function of oxygen only (e. g., Moore and Doney, 2007; Schmittner et al., 2008), sometimes involving switches to prevent nitrate from depletion (e. g., Moore and Doney, 2007; Assmann et al., 2010; Ilyina et al., 2013). A recent version of PISCES involves a complex parameterization for denitrification of dissolved organic carbon, depending on oxygen and different forms of dissolved inorganic nitrogen (Aumont et al., 2015). Otherwise, more detailed models of oxidant and substrate cycles are usually run on smaller spatial and shorter temporal scales (e. g., Gutknecht et al., 2013). On the other hand, given the long residence times of nitrogen mentioned above, model integration times of decades to centuries (e. g., experiments carried out by Moore and Doney, 2007; Landolfi et al., 2013) may not be sufficient to fully examine effects of different parameterizations of nitrogen fixation and/or denitrification and their mutual interactions and feedbacks mediated by oceanic transport and mixing processes. 

Therefore, a correct representation of nitrogen fluxes in global biogeochemical ocean models is challenging as it entails a large variety of spatial and temporal scales, from cellscale biological-chemical interactions up to global circulation time and space scales. In order to investigate combined effects of small scale, biogeochemical processes such as oxidant affinities, and large-scale, long-term nitrogen budgets, we here present results from an efficient, global offline model for passive tracers transport coupled to a detailed pelagic model of combined phosphorus, nitrogen, and oxygen cycles, simulated over millennia. Spin-up times of the model experiments of several thousand years are sufficient to draw conclusions about feedbacks and the net impact of these processes on the global nitrogen distribution, budget, and fluxes. 

We note that the model is meant as a first step towards a global model that includes aerobic and anaerobic remineralization in the oceanic (both pelagic and benthic) environment, including all elemental cycles associated with these processes. To address the abovementioned problems in a stepwise manner, we here first focus on the sensitivity of nitrogen budgets to remineralization length scale, oxidant affinity of remineralization under oxic and suboxic conditions, temperature dependency, and stoichiometry of nitrogen fixation. Benthic denitrification is not explicitly included in the current model that focuses on pelagic processes. We also at this point neglect the iron cycle, and its effects on autotrophic biomass production, particularly on that of cyanobacteria. Both will be addressed in follow-up versions of this model. 

The paper is organized as follows: after an introduction into the model structure and experimental setup, we present model results against the background of observed concentrations, inventories, and fluxes. We finally discuss some characteristic features of the model that may shed some light on the dynamics of other model simulations as well as oceanic processes. Nitrogen is coupled to the biogeochemical, phosphorus-based core via fixed stoichiometric ratios, which is indicated by purple borders around these compartments. Fixed nitrogen loss and gain are indicated by minus and plus signs. Their magnitude for different model setups, is shown in Figs. 10 and 9, as well as Table 2.
section_name: Introduction

section_content: The model simulates the cycling of nutrients (N), phytoplankton (P), zooplankton (Z), detritus (D) and dissolved organic matter (DOM), as described for phosphorus by Kriest et al. (2012). It further parameterizes burial of organic matter (as phosphorus and associated elements) at the sea floor, and, to close the mass budget, its resupply via river runoff (Kriest and Oschlies, 2013). Hereafter, we refer to these previous model types as CTL (no burial) and BUR (burial at the sea floor). In Sect. 2. 1, we only present a brief overview of the P-based pelagic core of the model, which is common to CTL and BUR, and refer the reader to Kriest et al. (2012) and Kriest and Oschlies (2013) for a detailed presentation and analysis of these models. 

In these previous, phosphorus-only model versions, remineralization of organic matter back to inorganic nutrients (i. e., phosphate) continued even in the absence of oxygen, thereby parameterizing some implicit, non-oxygen oxidants. We have now extended this model to include oxidantdependent remineralization, together with an explicit representation of the nitrogen cycle, where both oxygen and nitrate may act as final electron acceptors for oxidation of organic matter. Consistent stoichiometric balance of all organic substrates, products, and oxidants is achieved following Paulmier et al. (2009). In Sect. 2. 2 we describe in detail this Model of Oceanic Pelagic Stoichiometry (hereafter called MOPS; see also Fig. 1 for a sketch illustrating the model structure, and the different elemental cycles). 

The biogeochemical model description will be supplemented by a presentation of the different experiments carried out with the biogeochemical model, and by a brief introduction into the Transport Matrix Method (TMM), the underlying physical forcing, and some technical aspects.
section_name: Model description and experiments

section_content: We assume that different biogeochemical processes operate in different domains of the water column, with fast and dynamic turnover of phosphorus in the upper ocean layers, and a slow turnover of phosphorus below. To specify processes operating only in the euphotic zone (0-100 m, or k ≤ 6 for the chosen z level circulation model with level index k), we use the symbol H e (k) ≡ H (k e -k), where H (k) is the Heaviside step function. H e (k) is 1 in the euphotic layers, and 0 outside.
section_name: The biogeochemical phosphorus core

section_content: Phytoplankton (P) light limitation f (I ) is parameterized following Evans and Parslow (1985), using a globally uniform initial slope of the P-I curve of 0. 025 (W m -2 ) -1 d -1 (see also Kriest et al., 2010). Its maximum growth rate µ PHY depends on temperature, with µ PHY (T ) = 0. 6 e T 15. 65 (following Eppley, 1972, in the notation by Schmittner et al., 2008) and T given in degrees Celsius. We assume that the most limiting resource determines the phytoplankton growth rate. Thus, phytoplankton growth is parameterized as PP = µ PHY PHY min (f (I ), g(X 1, X 2,. . )) where g(X 1, X 2,. . ) is a Monod function of only phosphate in models CTL and BUR, and a function of both phosphate and nitrate in MOPS (see Sect. 2. 2 below). Phytoplankton experience a linear loss term of λ PHY = 0. 03 d -1, and are grazed by zooplankton. Zooplankton grazing (G) is described by a Holling III function, i. e., via a quadratic dependence on phytoplankton, a maximum grazing rate of µ ZOO = 2 d -1, and half-saturation constant of K ZOO = 0. 088 mmol P m -3. Only a fraction of grazing, ZOO = 0. 75, is effectively ingested, the rest is released again via egestion. Zooplankton experience a quadratic mortality κ = 3. 2 (mmol P m -3 ) -1 d -1. We assume that a fraction σ DOP = 0. 15 of egestion, zooplankton mortality, and phytoplankton loss is released as dissolved organic phosphorus (DOP), the rest becomes detritus. Zooplankton further experience a linear loss term of λ ZOO = 0. 03 d -1.
section_name: Euphotic zone

section_content: Phytoplankton and zooplankton die with a constant mortality rate of λ PHY = λ ZOO = 0. 01 d -1, again only when present above the lower concentration threshold P min = 10 -6 mmol P m -3. The dead organisms immediately disintegrate to DOP. 

DOP in all layers remineralizes with a constant rate λ DOP = 0. 17/360 d -1, but only when present above lower limit P min. In models CTL and BUR this is simply a linear function of organic substrate:

Likewise, modeled detritus remineralizes with a fixed rate λ DET = 0. 05 d -1 :

In MOPS, remineralization of organic matter depends on the availability of oxidants, as detailed below. Remineralized DOP and detritus directly feed into the phosphate pool; therefore, phosphate gain due to remineralization is

We assume that the sinking speed of detritus increases linearly with depth, according to w(z) = a z, where z is the center of a layer. We note that in steady state, with constant λ DET and in the absence of any other processes, this parameterization can be regarded as equivalent to the so-called Martin (power law) curve of particle flux: F (z) ∝ z -b, with the exponent b given by b = λ DET /a (see Kriest and Oschlies, 2008, for a detailed discussion). Previous model experiments have shown a large sensitivity of global tracer distributions to variations in the particle flux profile (Kriest and Oschlies, 2008; Kriest et al., 2012; Kriest and Oschlies, 2013). We address the potential sensitivity of model results to changes in particle sinking speed via a set of experiments, as detailed further below.
section_name: All layers

section_content: The air-sea gas exchange (top layer only) is parameterized following the OCMIP-2 protocol, with piston velocity and saturation computed from a monthly mean wind speed, temperature, and salinity derived from the MIT ocean model, and interpolated linearly onto the current time step. Oxygen also changes due to photosynthesis and remineralization, using a fixed stoichiometric ratio of R -O 2 :P = 170 mmol O 2 /mmol P. While in earlier models CTL and BUR, for O 2 ≥ 4 mmol O 2 m -3 oxygen decreased proportionally to the concentration of organic matter:

in MOPS we will consider the dependence of remineralization on oxygen explicitly, as described on in Sect. 2. 2 below.
section_name: Oxygen and air-sea gas exchange

section_content: A fraction of detritus deposited at the sea floor (at the bottom of the deepest vertical box) is buried instantaneously in some hypothetical sediment. Non-buried detritus is resuspended into the water column, where it is treated as regular detritus again. The phosphorus budget is closed on an annual timescale through resupply via river runoff. For more details about burial and its budget closure see Kriest and Oschlies (2013).
section_name: Benthic exchange

section_content: The following equations describe the source-minus-sink terms for the earlier, phosphorus-based models presented in Kriest et al. (2010), Kriest et al. (2012), and Kriest and Oschlies (2013), including oxygen, coupled via a constant stoichiometry. Remineralization under oxic conditions is indicated by terms S R X, where X is the affected tracer. Additional terms that account for processes in suboxic environments in model MOPS are termed S D X, and explained below. Superscript M denotes mortality terms of phyto-and zooplankton.
section_name: Source-minus-sink terms

section_content: While in the earlier models, CTL and BUR described above aerobic remineralization was not limited by oxygen (see Eqs. 3 and 4) 7 ) and ( 10 )-( 12 ). Beside the addition of nitrate as seventh state variable, in MOPS we further account for multiple nutrient limitation in term PP of Eqs. ( 7 ) and (8). These changes will be explained in detail in the following section.
section_name: S(PO

section_content: At this stage, we implement N in the simplest possible way, by considering only nitrate, but neither nitrite nor ammonium, as additional nutrient. Further, we assume that all biological components (i. e., phytoplankton, zooplankton, detritus and DOM) have a constant stoichiometry, given by d = N : P = 16. Thus, for MOPS we only add one additional state variable to the phosphorus-based core model, namely nitrate. Adding nitrogen in this way requires the parameterization of three different processes: multiple nutrient limitation of phytoplankton growth, nitrogen fixation, and heterotrophic nitrate reduction under suboxic conditions (hereafter loosely termed "denitrification"):
section_name: MOPS: oxidant-dependent remineralization and the nitrogen cycle

section_content: We assume a minimum function for the co-limitation of phytoplankton growth by phosphate and nitrate. First, we define the limiting nutrient L via

If L > 10 -6, we then evaluate the combined light and nutrient limitation function in analogy to Kriest et al. (2012) ; i. e., we define total phytoplankton production as

where µ PHY (T ) and f (I ) are the maximum growth rate and light limitation as defined above, and K PHY = 0. 03125 mmol P m -3 is the half-saturation constant for nutrient uptake. This changes term PP in Eqs. ( 7 ) and (8).
section_name: Multiple nutrient limitation

section_content: Remineralization of organic matter by (implicit) bacteria depends on organic substrate concentration and the availability of oxidants. In MOPS, dependency of both aerobic and anaerobic remineralization on organic substrates (detritus, dissolved organic matter) is parameterized as a first-order process. We further assume that these processes depend on oxidant availability following a saturation curve:

where we only consider oxygen above a certain threshold

are the minimum concentration and half-saturation constant for the heterotroph's uptake of oxygen in setup BASE, respectively (see below and Table 1, for sensitivity experiments related to oxidant affinity). To restrict oxygen consumption per time step, we first calculate the theoretical oxygen demand for respiration u T O 2 :

where λ DET = 0. 05 [d -1 ] and λ DOP = 0. 17/360 [d -1 ] are the remineralization rates of detritus and dissolved organic phosphorus, respectively, as defined above (Sect. 2. 1). t = 1/16 d is the time step length of the biogeochemical model. R -O 2 :P = 170 again denotes mole oxygen required per mole phosphorus remineralized (see Table 1 ). As above, we restrict the minimum detritus and DOP concentration for the onset of remineralization: DET * = max(0, DET -P min ) and DOP * = max(0, DOP-P min ). The aerobic decay rate limitation is then

Instead of applying Eqs. ( 3 ) and ( 4 ), the decay of DOP and detritus under oxic conditions is now defined by

As above, decay of organic matter results in a gain of phosphate, as defined by Eq. ( 5 ), and a loss of oxygen as in Eq. ( 6 ). 

In addition to aerobic decay, we now assume that under suboxic conditions (here defined as regions where O 2 * < 36 mmol O 2 m -3 ) denitrification replaces aerobic remineralization with R -NO 3 :P = 0. 8×R -O 2 :P -d being the nitrate demand of denitrification (Paulmier et al., 2009, their Table 1). 

Stoichiometric budget considerations suggest that on long timescales, and over large spatial areas both anaerobic ammonium oxidation and canonical denitrification (here: heterotrophic reduction of nitrate to dinitrogen) result in the same nitrogen loss, and are undistinguishable in the model context (see Appendix A). We therefore refrain from resolving these processes explicitly, but refer to denitrification as oxidation of organic matter via reduction of nitrate to dinitrogen. 

The rate of fixed nitrogen removal ultimately depends on the availability of organic matter, nitrate, and the absence of oxygen. Motivated by recent observations of anammox and nitrate reduction under rather high (up to 25 mmol O 2 m -3 ) ambient oxygen concentrations, and by the potential cooccurrence of both aerobic (ammonium oxidation) and anaerobic processes (Kalvelage et al., 2011), we parameterize a gradual increase of denitrification with decreasing oxygen concentrations. By using Eq. ( 15 ), we ensure a smooth transition between regimes of low and high oxidant concentrations; again we define a quadratic rate limitation of this process, but reduce it by the inverse oxygen consumption rate:

where NO 3 * = max(NO 3 -min NO 3, 0); min NO 3 = 4 mmol N m -3 and K NO 3 = 32 mmol N m -3 are the minimum concentration and half-saturation constant for the denitrifiers' uptake of nitrate in setup BASE, respectively. 

We note that the choice of the half-saturation constants of aerobic and anaerobic processes is rather arbitrary; however, as shown above the only gradual decrease of nitrate reduction (likewise for anammox) under increasing oxygen concentrations is supported by observations. To account for the uncertainty in oxidant sensitivity of oxic and suboxic remineralization, we have carried out experiments with different half-saturation constants, which are explained below (see also Table 1 ). 

As for oxygen, we restrict the use of nitrate to the amount available:

The rate limitation of anaerobic decay is then

Organic matter decay under suboxic conditions (denitrification) can thus be described by

The associated phosphate changes are described in analogy to Eq. (5):

Aerobic decay of organic matter increases nitrate according to stoichiometric ratio d, denoted as S R NO 3. Under suboxic conditions, there is further a decrease of nitrate due to denitrification, S D NO 3 :

with stoichiometric ratio R -NO 3 :P as defined above.
section_name: Oxidant-dependent remineralization

section_content: The simulated removal of fixed nitrogen via denitrification is counteracted by supply of fixed nitrogen through pelagic cyanobacteria in the euphotic zone. Unfortunately, sampling species such as Trichodesmium presents methodological difficulties (Breitbarth and LaRoche, 2005). For example, assuming one nifH gene copy per cell results in calculated abundances up to 10 5 cells L -1 for station ALOHA near Hawaii (Goebel et al., 2007), corresponding to concentrations up to about 0. 35 mmol C m -3. However, this is an order of magnitude higher than previous estimates from microscopic counts of Trichodesmium in this region (Letelier and Karl, 1996). Despite considerable recent efforts, global data sets for biomass remain relatively sparse, with only 2280 biomass estimates when gridded onto a 1 • × 1 • grid with 33 vertical layers (Luo et al., 2012). 

Given these uncertainties, and the sparse observations of cyanobacteria biomass, we refrained from explicit simulation of cyanobacteria, and instead assumed zero net growth, with immediate release of fixed nitrogen as nitrate (which also assumes immediate nitrification in our model that does not resolve different inorganic nitrogen species). By implicitly assuming constant cyanobacteria biomass, and a relaxation of the nitrate : phosphate ratio via immediate release of fixed nitrogen, our approach is similar to the one described by Maier-Reimer et al. (2005) and Ilyina et al. (2013). 

Temperature dependence is parameterized based on observations of the diazotrophic filamentous cyanobacterium Trichodesmium, which is responsible for a large fraction of global nitrogen fixation (Breitbarth et al., 2007). Instead of taking the fourth-order polynomial fit of maximum growth rate to temperature presented in the study by Breitbarth et al. (2007), we have approximated their function for maximum growth by a second-order polynomial, fitted over 20-34 • C. The resulting function has a maximum rate of 0. 2395 d -1 at T = 26. 82 • C. We then normalized the temperature-dependent growth rate by the maximum rate. By doing so, we obtain a T -limitation curve that varies between 0 and 1 for a temperature range of 20-34 • C, and is zero elsewhere:

To examine the effect of temperature limitation, we carried out an experiment where we skipped the temperature dependency for nitrogen fixation (denoted as NFixNoTemp; see also Table 1 and below for more details). 

Using a rather geochemical approach to restore fixed nitrogen towards the Redfield stoichiometry, we further assume that -in the presence of phosphate -nitrogen fixation is regulated by the nitrate : phosphate ratio:

where d * is the ratio that sets the stoichiometry of this process (see below). Having defined the temperature and nutrient regulation of nitrogen fixation, both of which are allowed to vary between 0 and 1, we finally assign a rate of maximum nitrogen uptake by cyanobacteria: under optimum conditions of T = 26. 82 • C and nitrate phosphate, we set

. This value is within the range of many observed rates: oceanic rates usually range between 0 and ≈ 0. 2 nmol N L -1 h -1, with some higher values up to 3. 1 nmol N L -1 h -1 (Mulholland, 2007, we exclude one extremely high value from the Arafura Sea). The compilation by Staal et al. (2007, their Table 5) gives values between 0. 07 to 17. 3 nmol N L -1 d -1, with most values in the lower range. 

A value of 2 nmol N L -1 d -1 is also encompassed by the recent, comprehensive data compilation by Luo et al. (2012). However, in blooms or in certain incubations, rather high values were found (up to 30 nmol N L -1 d -1 ; Goebel et al., 2007; Staal et al., 2007; Kitajima et al., 2009). We thus consider our maximum value of 2 nmol N L -1 d -1 a more conservative estimate. 

Note that combining our maximum uptake of 2 nmol N L -1 d -1 with the maximum normalized growth rate of 0. 2395 d -1 implies a constant cyanobacteria concentration of 8. 4 nmol N L -1, equivalent to about 0. 05 mmol C m -3, when using a C : N ratio of 6. 3. This concentration is lower than an estimate of 0. 35 mmol C m -3 from nifH gene copies (Goebel et al., 2007), but close to the microscopic estimates of Trichodesmium cells (Letelier and Karl, 1996). The recent compilation by Luo et al. (2012) suggests very high biomass (up to ≈ 10 mmol C m -3 ) in the Caribbean Sea, and at the surface of the tropical Atlantic and the Arabian Sea. However, many open-ocean values especially in the Pacific Ocean and/or deeper layers are rather low (< 0. 1 mmol C m -3 ). Thus, our implicit biomass estimates are in line with the observed estimates. 

Choosing constant cyanobacteria concentrations, and the abovementioned dependence on nitrate : phosphate ratios, we therefore follow a rather pragmatic approach that parameterizes nitrogen fixation with immediate release of fixed nitrogen as nitrate as antagonist to the fixed N loss during denitrification (in mmol N m -3 d -1 ):

Note that nitrogen fixation with immediate remineralization of the fixed N to nitrate theoretically requires 1. 25 moles oxygen = 2. 5 oxygen atoms per mole nitrate produced. So far, this is not included explicitly in the model; instead, we assume that this loss will be compensated for immediately through air-sea flux of oxygen.
section_name: Nitrogen fixation

section_content: Combining the abovementioned processes and interactions, with changes of nitrate due to basic biological processes in the euphotic zone (phytoplankton growth; zooplankton grazing), the time rate of change for nitrate is
section_name: Nitrate as a new state variable

section_content: Many of these newly added processes occur only in certain areas, sporadically, and at low rates. Their parameterizations are therefore subject to uncertainties. Here we try to address some of the uncertainties by carrying out a set of sensitivity experiments, in which we simulate the different model versions with changed parameters, over the same spatial and temporal domain, and use the same physical forcing and setup. These sensitivity experiments address two alternative parameterizations of nitrogen fixation, and two experiments with changed oxidant affinities. As earlier studies have shown an impact of particle sinking speed on large-scale tracer distributions, we carry out these experiments before the background of different particle flux profiles.
section_name: Sensitivity experiments

section_content: In the standard setup BASE, nitrogen fixation relaxes the nutrient ratio to the stoichiometric relation used for the other biogeochemical processes, namely d * = d = 16 (see Table 1 ). We further carried out a sensitivity experiment NFixStoich, where this process depends on a ratio d * = 14. 28, which is the observed globally averaged nitrate : phosphate ratio. Another sensitivity experiment (NFixNoTemp) employs temperature-independent nitrogen fixation.
section_name: Nitrogen fixation

section_content: In the standard model setup BASE, aerobic remineralization depends on the half-saturation constant (affinity) for oxygen uptake of K O 2 = 8 mmol O 2 m -3, while denitrification is inhibited by oxygen, and is further determined by a half-saturation constant for nitrate uptake, K NO 3 = 32 mmol N m -3 (see Table 1 ). The latter value is at the upper end of experimentally derived estimates (Jensen et al., 2009). A sensitivity experiment DenHigh investigates the model's response to a higher nitrate affinity of denitrifying bacteria, as represented by a lower half-saturation constant for nitrate (K NO 3 = 8 mmol N m -3 ), while another experiment RemHigh examines the model's sensitivity to concomitant increase in nitrate and oxygen affinity of aerobic and anaerobic remineralization (K NO 3 = 8 mmol N m -3 and K O 2 = 2 mmol O 2 m -3, together with a lower oxygen level for the onset of denitrification; see Table 1 ).
section_name: Oxidant affinity of remineralization

section_content: As noted above, in our model we assume that particle sinking speed w increases linearly with depth z: w = w(z) = a z, which, under equilibrium conditions and absence of other processes beside sinking and remineralization, would result in a particle flux profile defined by F (z) ∝ z -b, where b = λ DET /a (see also Kriest and Oschlies, 2008). For all five model setups (BASE, DenHigh, RemHigh, NFixNoTemp and NFixStoich), we carried out experiments with varying sinking speed parameter a. To facilitate direct comparison to the particle flux profiles determined by Martin et al. (1987) or Van Mooy et al. (2002), we keep λ DET constant and express changes in parameter a in terms of b. The changes we applied would correspond to a variation of b between 0. 6435 (fast sinking) over 0. 858 (medium) to 1. 0725 (slow sinking) for all model experiments, and additionally to b = 0. 429 (very fast sinking) and b = 1. 287 (very slow sinking) for the refer-ence setup BASE. We note that in MOPS, due to reduction of accomplished remineralization by lack of oxidants (factors s O 2 and s NO 3 in Eqs. 19 and 24, respectively), the local, effective Martin exponent b may be smaller than its nominal value.
section_name: Particle sinking speed

section_content: Global model simulations were carried out using the TMM method described in detail by Khatiwala (2007). 

In the TMM method the three-dimensional (3-D), advective-diffusive transport (including all sub-grid scale parameterizations) of an ocean circulation model is represented in the form of sparse transport matrices (TM; see also Khatiwala, 2007, for more details on TMs extraction), which are then used to move any passive tracer. This method allows for the efficient and convenient testing of many different biogeochemical models in an offline mode. We used 12 monthly mean TMs derived from the Estimating the Circulation and Climate of the Ocean (ECCO) project, which provides circulation fields that yield a best fit to hydrographic and remote sensing observations over a 10-year period, on a spatial grid of 1 • × 1 • horizontal resolution with 23 vertical levels (Stammer et al., 2004). We note that this model may exhibit too efficient ventilation (Graven et al., 2012), a feature that may be common for many coarse-scale circulation models. TMs derived from ECCO were also applied by Kriest and Oschlies (2013), albeit with a different temporal resolution. We note that tests with the reference biogeochemical model and a shorter time step showed only small differences, particularly when compared to the impact of biogeochemical parameters. 

The nitrogen cycle model MOPS presented here operates on many timescales, determined by the relatively fast biological surface processes, slower deep remineralization, and the global ocean circulation. In addition, if areas of fixed nitrogen gain and loss are spatially segregated, timescales of physical transport between these regions are of importance. To ensure full equilibration of all processes involved, we spun-up the coupled system over 9000 years, using two time steps per day for tracer transport, and 16 time steps per day for the calculation of biogeochemical source-minussink terms (Kriest and Oschlies, 2013, used 1/8 d for tracer transport, and 1/64 d for biogeochemical processes, simulated over 3000 years). A spin-up time of 9000 years should be sufficient to ensure almost complete pelagic tracer equilibration, which may be as long as ≈ 10 000 years (see also Wunsch and Heimbach, 2008), but can depend on the surface boundary condition (Wunsch and Heimbach, 2008; Primeau and Deleersnijder, 2009; Siberlin and Wunsch, 2011). In our case, near-equilibrium conditions are also reflected by a very small disequilibrium of total fixed nitrogen inventory, which, in our model experiments, changes less than 0. 003 % between year 8900 and 9000. 

We initialized all experiments from observed distributions of phosphate, oxygen, and nitrate (monthly mean values for January above 500 m, and annual-mean values below), as provided by Garcia et al. (2006a, b). Initialization of other tracers has been carried out as in Kriest and Oschlies (2013). By adding a nitrogen cycle to the prior phosphorusbased model, nonlinear switches in oxidants are introduced. 

Though it has, to our knowledge, never been shown that marine biogeochemical models can exhibit multiple steady states, other components of the Earth system, such as ocean circulation, atmospheric circulation, land ice and terrestrial vegetation all have been found capable of displaying situations with multiple steady states. To test whether indeed the final model state is independent of its initial conditions, we carried out a number of test simulations with very different initial biogeochemical tracer distributions with model setup BASE, with a spatial resolution and circulation as in Kriest et al. (2012). These did always reach the same steady state, even if the model was started from near-zero oxygen and nitrate, indicating the robustness of the model results shown below. The final tracer distribution of the model being independent of its initialization, the data sets of observations (that went into initialization) may be used for model evaluation. 

Model simulations were mostly carried out on an SGI Altix ICE 8200 Linux cluster at the North-German Supercomputing Alliance (HLRN, www. hlrn. de), using different versions 3. x of PETSc (Portable, Extensible Toolkit for Scientific Computation; www. mcs. anl. gov/petsc/). Using a total of 256 cores, 3000 years of simulation took between 10 and 12 h; i. e., a full spin-up could be performed in less than 36 h. We note that experiments on other hardware/platforms (e. g., Intel Sandybridge; CrayXC30) and with different versions of PETSc (up to version 3. 5) resulted in very small -if anydifferences in model results.
section_name: Circulation model, spin-up, and initialization

section_name: Results

section_content: Figure 2 indicates that, depending on particle sinking speed, model spin-up times of at least a few millennia are necessary in order to fully equilibrate the different processes in the model. Within the first few decades, all model setups started from observed tracer concentrations initially lose some oxygen. The loss continues for experiments with moderate or fast-sinking speed, until oxygen approaches a global average equilibrium value that is about 16 to 18 mmol O 2 m -3 lower than that observed for the fast-sinking scenario, and between 8 and 12 mmol O 2 m -3 lower for moderate sinking speed. Model experiments with slower sinking speed show a different transient response: after the initial decline of oxygen, global average oxygen content increases again, until it almost approaches its initial value (setup BASE), or even exceeds it (setups DenHigh and RemHigh). Results of setup NFixStoich are very similar to those of BASE, and not shown here. 

A change in sign is also exhibited by the transient behavior of the globally averaged nitrate concentration in all slow scenarios of setups BASE, NFixStoich, DenHigh, and RemHigh (Fig. 2, lower panels). All model setups investigated exhibit an initial, very rapid loss of nitrate within the first few years (see also Fig. 3, left panels). This initial loss of nitrate can be attributed to surface processes, converting the inorganic dissolved tracers into organic forms (dissolved and particulate organic matter) as can be deduced from the almost constant total N found in almost all model configurations (Fig. 2, lower panels). In addition to this process, the loss of nitrate by denitrification, and its supply via nitrogen fixation affects the inventory of this tracer on longer timescales. MOPS' nitrate inventory starts to increase within the first few centuries, until it finally approaches an equilibrium value which is slightly (slow) or up to 2. 5 mmol N m -3 (fast) higher than at the beginning. An exception to this is experiment DenHigh, where, due to the higher nitrate affinity of denitrification, equilibrium nitrate in the slow sinking scenarios is slightly less than initially prescribed. Because in DenHigh more nitrate is used for oxidation, oxygen is higher than in setup BASE. 

Total fixed nitrogen is not affected by the abovementioned conversion of inorganic to organic nitrogen; therefore it lacks the initial rapid decline exhibited by the nitrate inventory (Fig. 2, lower panels). Total nitrogen first declines slightly in all slow model scenarios, except setup NFixNoTemp, which, due to the widespread occurrence of nitrogen fixation, may quickly compensate the loss through denitrification. After this initial decline, fixed nitrogen, like nitrate, equilibrates to values higher than initially prescribed -again with the exception of the slow scenario of model setup DenHigh. 

Most important for the transient behavior of fixed nitrogen are fluxes in the eastern tropical and subtropical Pacific, as illustrated in Fig. 3. Phosphate averaged over two regions EEP (here east of 140 • W, ±10 • ) and LLP (low latitudes of the tropical and subtropical Pacific; ±40 • latitude) shows only small initial variations (less then 0. 1 mmol P m -3 variation all experiments; no figure). In contrast, nitrate in the EEP decreases strongly within the first few hundred years, most strongly for the slow sinking scenario, down to a deficit of 3. 7 mmol N m -3 (Fig. 3, upper left panel). Nitrate then approaches a minimum, first in scenario fast (around year 200). The slow sinking scenario takes much longer to reach that minimum (about 400 years). The initial decline of nitrate is subsequently followed by a slight increase during the next centuries in all scenarios. In the LLP, nitrate exhibits a far less pronounced transient (Fig. 3, mid left panel), comparable to that of phosphate. The different transients of both dissolved tracers phosphate and nitrate are mirrored in the nitrate : phosphate ratio, which shows a strong decline (by about 0. 5-1 units) in the EEP within the first centuries, but 0 ), where x is either global average oxygen or nitrate). Average tracers are calculated from snapshots of day 360 every 10th year within the first 200 years, and every 100th year thereafter. Line colors denote different sinking speeds: red -fast, black -medium, blue -slow. Model identifier is shown on top of each column. For better visibility on all timescales, we present the time axis on a logarithmic scale. Note that as models were started from observed global average tracers, any deviation from the dashed (zero) line also denotes the model bias. 

only small variations in the LLP (Fig. 3, second from left panels). 

The strong variation of nitrate in the EEP is caused by vigorous denitrification, from initially ≈ 40 Tg N yr -1 to ≈ 10-50 Tg N yr -1 by year 900, while nitrogen fixation proceeds at a constant, low level of < 10 Tg N yr -1 (Fig. 3, third left from panels). The difference between the two fluxes in this region results in a net loss of nitrogen through local biogeochemical processes (Fig. 3, right panels). Despite this loss, fixed nitrogen in the EEP remains relatively constant after some centuries. The almost balanced fixed nitrogen budget, in the presence of high (biogeochemical) nitrogen loss, can be attributed to supply from the adjacent area LLP, which shows high nitrogen fixation, but relatively low denitrification. The resulting positive (between ≈ 10-40 Tg N yr -1 ) net N flux into the LLP region matches the loss in the EEP. Thus, after the first millennium, the more or less stable fixed nitrogen inventories in both EEP and LLP regions can be explained by transport between the two regions, which each display local disequilibria between dominant N loss (EEP) and N gain (LLP). The interactions between these two regions are also reflected in the initial transient pattern of global fluxes, which more or less mirror the combined fluxes of both regions (bottom panels of Fig. 3 ).
section_name: Transient concentrations, inventories, and fluxes

section_content: After the 9000 year spin-up, the final (near steady state) model solution is independent of its initialization, and solely reflects the combined effects of biogeochemistry and circulation. We can therefore use a comparison to observed tracers and fluxes in order to assess model skill and performance. In the next sections, we first examine model fit against observations of dissolved inorganic tracers. Examination of nitrogen fluxes against the quite sparse observational data sets can provide a first insight into the adequacy of some model parameterizations. Comparison to more robust, bulk diagnostics helps to examine the general model behavior. We finally combine some of the model-data comparisons to a total, global misfit function that should help to decide between the different model setups and scenarios.
section_name: Steady-state concentrations, inventories, and fluxes

section_content: In steady state all models exhibit similar volume distributions of phosphate (see left panels of Fig. 4 for model setups BUR and BASE; the other model setups show similar results, see auxiliary Fig. 1 ), which can be attributed to the fact that they are based on the same phosphorus core. Variations in sinking speed play only a small role for the volume distribution of phosphate. Surprisingly, the models also show very little dif- ference in the volume distribution of oxygen, which matches observations quite well particularly for the slow sinking scenarios. Even the introduction of the nitrogen cycle, together with oxidant-dependent remineralization, does not strongly affect the distribution of this tracer. Nitrate is simulated quite well by the slow sinking scenarios of model setups that explicitly include this tracer (MOPS). The explicit simulation of nitrate in MOPS results in an even better representation than nitrate diagnosed a posteriori from phosphate simulated by BUR, times a constant stoichiometric ratio of 16 (see straight lines in upper right panel of Fig. 4 ). Replacing this ratio for conversion by the observed global mean ratio of 14. 28 (dashed lines in right panels of Fig. 4 ) results in a better fit for all model experiments, yet this latter, empirically derived nitrate diagnostic provides only a weak constraint on model performance, because of its dependency on observations. 

A common way to look at model performance with respect to observed tracers, is to combine information about simulated and observed standard deviations, correlation co-efficient (R), and centered (unbiased) root mean square error (hereafter referred to as RMSE ) in a so-called Taylor plot (Taylor, 2001). Figure 5 shows these diagnostics for phosphate, oxygen, and nitrate of model setups CTL, BUR, and BASE. Obviously, all models deteriorate with respect to phosphate with increasing sinking speed, as indicated by too high a standard deviation, decreasing correlation coefficient, and RMSE. However, differences for slow settling speeds are rather small. Likewise, except for extreme sinking velocities, model results are quite similar when examining the fit to observed oxygen. However, results are more variable (between model types, and with respect to different metrics) when nitrate is considered, either diagnosed from phosphate times 16 for models CTL and BUR, or simulated explicitly for model BASE. First, for all model setups we find a quite strong overestimate of variance, and decrease in fit (RMSE and RMSE ) especially for fast-sinking velocities. Further, model BASE at first sight seems to exhibit a far worse fit to observations (with respect to correlation coefficient R and RMSE ) than BUR. This is in striking contrast to the right panels d-f ). Grey bars denote the corresponding observations (Garcia et al., 2006a, b). Colors denotes different sinking speeds. Blue: slow sinking; black: medium sinking; red: fast sinking. Nitrate for model BUR (without nitrogen cycle) has been calculated from phosphate ×16. In addition, we present nitrate calculated from phosphate ×14. 28 for the slow sinking scenario as dashed blue lines. panels of Fig. 4, which indicate a better fit of nitrate simulated by BASE than nitrate diagnosed from BUR's phosphate. However, it is important to note that RMSE does not account for the bias in total nitrate concentration (see also Taylor, 2001; Jolliff et al., 2009). Therefore, although BUR matches the general pattern of nitrate distribution (via phosphate; panel A of Fig. 4 ), its average concentration does not match the observed average of ≈ 31 mmol NO 3 m -3, as indicated by its overestimate of ocean volume containing nitrate concentrations > 40 mmol NO 3 m -3 (see panel c of Fig. 4 ). As a result, model BASE for each sinking speed shows a better fit to observations with respect to the total RMSE (color scale of Fig. 5 ). In contrast to the normally used Taylor plot RMSE, which would favor model BUR over BASE, the RMSE includes both the match to the pattern and to total tracer inventory, which are, in our case, best reproduced by model BASE. Despite the overall good match of the global distribution of dissolved tracers to observations, models may differ in regions which are particularly sensitive to the combined effects of oxygen supply, sinking, and remineralization. For example, as shown above (Fig. 3 ), the eastern tropical Pacific seems to play a large role for of global fluxes, and thus global tracer inventories. To investigate this region further, in Fig. 6 we have a closer look at nutrients and oxygen averaged over ±5 • in the eastern Pacific. The analysis is similar to the one presented in Dietze and Loeptien (2013), but integrates over the upper 6500 m, and thus disregards mismatches in the vertical distribution of tracers. 

Simulated regional phosphate varies only slightly (< 0. 2 mmol P m -3 ) between the different model experiments (Fig. 6, left panel). On the other hand, simulated nitrate shows large variations of up to almost 20 mmol N m -3 towards the American coast, much larger than would be expected from the variations in phosphate and some typical stoichiometry (middle panel of Fig. 6 ). In this area, particularly the slow sinking scenarios of all model setups that include nitrogen strongly underestimate observed nitrate. At the same time these model experiments are in quite good agreement with observed oxygen, especially when simulated with a high affinity of denitrification for nitrate (Fig. 6, right panel). The fast-sinking model experiments that yield a somewhat better fit to observed equatorial nitrate content, however, systematically underestimate the equatorial oxygen inventory. Thus, all configurations of MOPS show mismatches for either oxygen or nitrate in this region, and no experiment is able to sufficiently represent both oxidants at the same time. 

Diagnosing nitrate from phosphate in the phosphorus-only model BUR yields a quite good agreement to observations between 150 and 110 • W when applying the global observed nitrate : phosphate ratio of 14. 28. East of 110 • W the model with slow sinking speed overestimates nitrate. Using a stoichiometric ratio of 16 as typical for marine phytoplankton composition (Anderson, 1995), and also typically used in numerical models, results in a strong overestimate of observed nitrate by model BUR over the entire transect. Hence, only with the help of observed nutrient ratios will this model agree with observed nitrate, impeding the use of this tracer for model evaluation.
section_name: Patterns of dissolved inorganic tracers

section_content: Simulated vertically integrated nitrogen fixation in steady state (year 9001) is ≤ 100 µmol N m -2 d -1 for large parts of the subtropical ocean. Higher values of up to 200 µmol N m -2 d -1 occur mostly in the Pacific, the western Atlantic Ocean, occasionally in the Caribbean Sea, and in the Arabian Sea and Bay of Bengal (upper panels of Figs. 7 and 8 ). Depending on the parameterization of sinking speed and biogeochemistry, the central Pacific Ocean is characterized by large areas with fluxes > 160 µmol N m -2 d -1. Slow sinking speed, especially when combined with high nitrate affinity of denitrification (setup DenHigh) increases steady-state nitrogen fixation, which can be attributed to the compensation of an enhanced fixed nitrogen loss (Fig. 7 ; see also below). Simulated nitrogen fixation rates mostly lie well within the range of earlier estimates for the open ocean (e. g., Mahaffey et al., 2005; Staal et al., 2007; Kitajima et al., 2009). Note that the model setups due to the here-employed maximum fixation rate of 2 nmol N m -3 d -1 cannot not reach some high values observed by Kitajima et al. (2009) and Staal et al. (2007). The comprehensive data set by Luo et al. (2012, their Fig. 6a ) shows enhanced integrated nitrogen fixation of Because denitrification is restricted to regions with low oxygen, it is not as widely distributed as nitrogen fixation. Areas of simulated denitrification are the Arabian Sea and Bay of Bengal, the eastern tropical and subtropical Pacific extending north-and southwards to latitudes of about 30 •, and the upwelling off Namibia and Angola (lower panels of Figs. 7 and 8 ). The model experiments simulate the highest vertically integrated rates in the latter two regions, where loss of fixed nitrogen can be as high as ≈ 10 mmol N m -2 d -1 (slow sinking scenario of model setup RemHigh). Because of the longer residence time of particles in midwater depths, simulated nitrogen loss increases with decreasing sinking speed. It further increases with nitrate affinity in setups Den-High and RemHigh. 

Maximum volumetric denitrification mirrors that of its vertical integral, and can be as high as 118 nmol L -1 d -1 (setup RemHigh with slow sinking). Highest modeled values occur in the eastern tropical and subtropical Pacific, followed by the upwelling off Namibia and Angola, and the Arabian Sea and Bay of Bengal, especially for slow sinking speed and/or high nitrate affinity. Simulated maximum values up to 43 nmol L -1 d -1 in the Arabian Sea/Bay of Bengal are higher than maximum observed rates of ≈ 25 nmol L -1 d -1 (Ward et al., 2009; Bulow et al., 2010), but most simulated values are in the range 2-20 nmol L -1 d -1, and thus quite in agreement with the observations. High rates of nitrogen loss (up to ≈ 150 nmol L -1 d -1 ) have been observed in the Benguela upwelling by Kuypers et al. (2005), and Kalvelage et al. (2011) even reported values of almost 500 nmol L -1 d -1. The model only simulates annual-mean rates up to 40 nmol L -1 d -1 in that region. 

Off Chile and Peru maximum simulated rates of nitrogen loss range between ≈ 30 and 120 nmol L -1 d -1 which is within observed rates of fixed nitrogen loss (Hamersley et al., 2007; Galan et al., 2009; Kalvelage et al., 2011). However, most observations in this region are related to the anammox process, with little or no indication for denitrification off Peru or northern Chile. Some recent work suggests sporadic, yet very high rates of denitrification in this region (Dalsgaard et al., 2012), perhaps due to the local, sudden input of organic matter. Comparison of our model results to these ob- servations does not seem straightforward, as our model does not explicitly distinguish between different oxidation states of nitrogen, or between the mechanisms of fixed nitrogen loss, neither does it resolve the fine temporal-spatial scales that characterize the nitrogen cycle along the coasts of Chile and Peru. In the discussion section we will examine this alleged model deficiency, in light of the various observations made in this region.
section_name: Patterns of fixed nitrogen sources and sinks

section_content: When ordered according to the global nitrogen throughput, global primary production is the largest of all global fluxes, followed by export production and particle flux at increasing depths (Fig. 9 ). Generally, models agree quite well with the observed range; there is some tendency by some models to overestimate particle flux in 1000 and 2000 m, when compared to observations. N 2 -fixation and pelagic denitrification are more than an order of magnitude smaller than fluxes into and out of the euphotic zone. As will be shown below, N 2 -fixation and pelagic denitrification depend strongly on the complex interplay between particle sinking, remineralization, and the formation and extent of OMZs. 

Particle sinking speed has a strong impact on the volume of OMZs. In models with slow sinking speed, organic matter prevails in midwater depths for a long time, causing a strong depletion in oxygen due to remineralization, and thus a large suboxic volume (Fig. 10, left panel). If particles sink faster, they may be buried in the sediment before being remineralized -therefore, these model experiments exhibit the smallest suboxic volume. Increasing nitrate affinity (setup DenHigh) causes less depletion in oxygen, because more nitrate is used for oxidation of organic matter. Increasing both nitrate and oxygen affinity (setup RemHigh) increases the OMZ volume. The criterion for the definition of suboxia (here: < 8 mmol O 2 m -3 or < 4 mmol O 2 m -3 ) has a strong impact on the evaluation of model misfit: for the lower criterion, the fast-sinking model experiments show the best fit to observations, whereas for the higher criterion the best sinking speed depends on the parameterization of oxidant affinity (e. g., slow sinking for model setup DenHigh vs. fast sinking for model setup RemHigh). Thus, a model evaluation based on the suboxic volume can yield very different results for different oxygen thresholds used to define suboxia (see also Gnanadesikan et al., 2013). 3, incl. citations therein; export production and particle flux). Fluxes have been converted using a C : P ratio of 117, and a N : P ratio of 16. Global fixed nitrogen fluxes (N 2 -Fixation, pelagic denitrification) as in Table 2. 

The effect of particle sinking on OMZ volume propagates into simulated global nitrogen fluxes. The larger suboxic volume of model scenarios with slow sinking speed triggers high rates of denitrification (Fig. 10, middle panel). In contrast, in the scenario with very fast sinking of model setup BASE no denitrification occurs, because of the lack of suboxic zones (see above). Moving to high nitrate affinity in setup DenHigh exhibits the highest global integrated nitrogen loss and balancing nitrogen fixation for any given sinking speed (see also Table 2 ). The effect of higher nitrate affinity is less pronounced in setup RemHigh. Depending on sinking speed, global nitrogen loss and balancing nitrogen fixation increase by some 20 % with respect to experiment BASE in sensitivity experiment RemHigh with higher oxidant affinity, and by 42 % when only nitrate affinity, but not oxygen affinity, is increased (setup DenHigh; see Table 2 ). Turning off the temperature dependence of the model's diazotrophs that constrains nitrogen fixation to warm surface waters (setup NFixNoTemp) results in simulated fixed nitrogen gain in high latitudes, far away from suboxic regions of nitrogen loss (see also panel a of Fig. 8 ), and therefore increases globally integrated nitrogen fluxes. 

Because MOPS simulates pelagic denitrification as only nitrogen loss process, this loss has, in steady state, to be matched by nitrogen gain via nitrogen fixation. In the following paragraphs we will therefore focus mostly on pelagic denitrification and its comparison to other observed or simulated estimates. We refrain from a detailed comparison to global nitrogen fixation (which, in reality, and under the assumption of homeostasis, would have to balance benthic denitrification as well). 

The overall magnitude of steady-state global integrals of pelagic nitrogen loss of 59-84 Tg N yr -1 diagnosed from our model simulations (Table 2 and Fig 2. Note that in the models, due to their intrinsic assumptions pelagic nitrogen loss equates with nitrogen gain through nitrogen fixation. Right: global nitrate : phosphate ratio of different models. Left and right: observations according to Garcia et al. (2006a, b). 

roughly with the model estimates of Moore and Doney (2007), but is considerably lower than the 140 Tg N yr -1 simulated by Oschlies et al. (2008). The range of estimates of nitrogen loss based on observations is even larger: Codispoti (2007) suggested that water-column denitrification should be even higher than 150 Tg per year, whereas substantially lower estimates of 52 to 81 Tg per year were obtained by Gruber and Sarmiento (1997), Galloway et al. (2004), Deutsch et al. (2004) and, more recently, Bianchi et al. (2012), Eugster and Gruber (2012 ), and DeVries et al. (2012, 2013). 

With global oceanic phosphorus being conserved in all model configurations, the simulated global nitrate : phosphate ratio depends on the variable global nitrogen inventory and hence on the simulated fixed nitrogen losses and gains, which in turn are sensitive to particle sinking velocities. The establishment of suboxic zones, and thus areas of enhanced denitrification, results in regions with a lowered nitrate : phosphate ratio (see also Fig. 3 ), finally with an effect on the steady-state global nitrate : phosphate ratio (right panel of Fig. 10 right panel ). An extreme case is experiment BASE with very fast sinking, in which no suboxia develops and consequently no denitrification occurs. In this case the assumptions implicit in the model description of nitrogen fixation, increase the global-ocean molar nitrate-to-phosphate ratio from its observed initial value of 14. 28 to 16, the stoichiometric ratio for aerobic processes (Fig. 10, right panel). 

The relation of higher nitrogen losses and gains corresponding to a lower global nitrogen inventory generally holds also between different model configurations, the exception being the sensitivity experiment NFixNoTemp with temperature-independent nitrogen fixation, which shows both higher denitrification and fixation and higher stoichiometric ratios when compared to the reference model setup BASE. Generally, simulated global nitrate-to-phosphate ratios of the slow to very slow sinking scenarios are closest to the observed ratio. Thus, a high pelagic turnover of nitrogen in suboxic areas, as mediated via slow sinking or high affinity of denitrification towards nitrate, pushes the simulated nitrate-to-phosphate ratio towards the observed ratio of global average nutrients.
section_name: Global fluxes and inventories

section_content: So far, we have analyzed the models' fit to different observed quantities separately. While most of the models discussed in this study fit spatial patterns of dissolved tracers about equally well, some differences arise for different particle sinking speeds. Also, as shown above (Fig. 5 ) and discussed by Jolliff et al. (2009), accounting for the model bias of non-conservative tracers such as nitrogen may serve as an important additional information on model skill, and help to discriminate between the different model types. Investigating the fit to observed global inventories of non-conservative tracers, to OMZ volume, and to global pelagic nitrogen loss as a function of different model parameters further indicates some mutually independent controls exerted by the different global diagnostics (see Fig. 10 ). 

In this subsection we examine to what extent the additional consideration of the nitrogen cycle can constrain individual process descriptions better than the standard metric of global distributions of phosphate and oxygen, pelagic and benthic fluxes of organic matter (in terms of phosphorus) and global oxygen inventory, applied previously to model BUR by Kriest and Oschlies (2013). Specifically, we add the spatial distribution of nitrate, the global nitrate inventory, and global pelagic nitrogen loss as three additional diagnostics. Individual observations of nitrogen fixation or nitrogen loss are not considered because of the sparsity of data sets, their bias towards certain regions, and differences between the model's intrinsic assumptions and observations (see above) that may complicate a direct comparison. Also not considered in our metric is the volume of suboxic regions, because of its pronounced sensitivity to the oxygen threshold used in the definition of suboxia (see above). 

Most phosphorus-based metric components of our new model of the pelagic nitrogen cycle MOPS are similar to those for the phosphorus-only model BUR presented by Kriest and Oschlies (2013); i. e., for dissolved tracer distribution, oxygen inventory, and particle flux (panels b, d, h-j of Fig. 11 ), we find the best fit for a model with a sinking speed reduced relative to the classical Martin parameters. Similarly, in MOPS the global nitrate inventory is matched best for very slow sinking speeds (Fig. 11g ). Simulated global pelagic fixed nitrogen loss, on the other hand, shows a best fit to observations for model simulations with medium sinking speed (panel a of Fig. 11 ), as do benthic burial and remineralization (Fig. 11, panels e, f, k, l). Therefore, with the given weights of global integrated or average properties (see also Kriest and Oschlies, 2013) the overall misfit function (Fig. 11b ) now favors model experiments with a power law flux exponent of 0. 858, as initially suggested by Martin et al. (1987). Unfortunately, it is not possible to distinguish between the different model setups BASE, NFixNoTemp, NFixStoich and RemHigh, as all of them perform more or less equally well with respect to the combined metric. Thus, with a misfit function that gives equal weights to all types of observations, and for the parameter ranges investigated, the parameterization of oxidant affinity or nitrogen fixation may be less important for model performance than the parameterization of particle flux. It remains to be investigated whether this still holds if additional data, for example, for nitrogen fixation in the (currently unexplored) eastern tropical Pacific, become available. Considering MOPS' skill it is encouraging, however, that this model, which confines remineralization to the oxic or nitratebearing regions, performs as well as model BUR with its assumed infinite supply of oxidants.
section_name: Added value of the nitrogen cycle for model assessment

section_content: Our model simulations indicate adjustment times of the order of millennia. Substantial imbalances of fixed nitrogen sources and sinks exist over the first few hundred years after initialization with observed tracer concentrations, reaching some 10 to 20 Tg N yr -1 after 200 years (lower right panel of Fig. 3 ). This can be explained by the slow coupling between circulation, deep remineralization, and spatial separation of regions of fixed nitrogen loss and gain, particularly in the EEP. After 3000 years, we still find a gain of fixed nitrogen of 0. 6-7. 2 Tg N yr -1 (depending on particle sinking speed) for model setup BASE, which then slowly declines to < 0. 1-0. 4 Tg N yr -1 in year 9001. Earlier model studies (Moore and Doney, 2007; Schmittner et al., 2008; Ilyina et al., 2013) with generally shorter spin-up periods reported imbalances between nitrogen fixation and denitrification of up to several Tg N yr -1. According to our model results, this may be indicative of those models not having spun-up for long enough to be in equilibrium. The above results indicate the need to spin-up models over a long enough time in order to evaluate the full response of the coupled physico-biogeochemical system. Analysis of model results from shorter spin-ups may result in nitrogen fluxes that reflect the transient (but not steady-state) characteristics of the system. 

After several thousand years the model solution is independent of its initialization, and solely reflects the combined effects of biogeochemistry and circulation, allowing the use of observed dissolved tracer distributions (that also went into initialization) for model assessment. Model results in steady state agree reasonably well with respect to observed tracer concentrations and fluxes, both with respect to local as well as to global quantities, giving some confidence in the model's representation of the large-scale dynamics of marine nitrogen. The good match of the model MOPS that explicitly calculates nitrogen is promising, as this model agrees with observations even in the presence of more constraints (oxidant affinity of remineralization; no implicit oxidants) than in previous, phosphorus-only model versions. 

However, in steady state all model experiments exhibit some deficiency in nitrate or oxygen in the EEP, which may be either due to an ill-defined biogeochemical model, deficient representation of the physics, or both. Beside the lack of so far undefined processes (iron cycle; benthic denitrification; variable stoichiometry), possible reasons for this mismatch are (1) lack of distinction between different oxidation states of nitrogen, (2) inappropriate, rather weakly constrained parameters for oxidant affinity, and (3) equatorial dynamics that is notoriously difficult to represent by medium-to coarse-resolution ocean circulation models (Dietze and Loeptien, 2013). In the following subsections we will have a closer look at these model features, which should also provide some insight into the dynamics of models outside this study. Given the different tracers and diagnostics simulated by the model, we will finally discuss briefly if, and how, these can be useful for constraining the global marine nitrogen (and oxygen) budget.
section_name: Discussion

section_content: As noted above, in some areas the model shows a quite good agreement to observed rates of denitrification. However, a di- rect comparison between simulated and observed rates of nitrogen loss is not always straightforward, as many observations refer to anammox, which is not explicitly resolved by the model. Further, some of the experiments measuring nitrogen loss have been carried out under high levels of nutrient additions, and may therefore be regarded as potential rates. For example, Kalvelage et al. (2011) found maximum anammox rates up to 108 nmol L -1 d -1 off Peru, but many rates were much lower. Potential anammox rates below 30 nmol L -1 d -1 have been observed by Galan et al. (2009), Hamersley et al. (2007), and Thamdrup et al. (2006). Most of these works found little or no indication for denitrifica-tion in the Pacific off Peru or northern Chile, but recent work by Dalsgaard et al. (2012) suggested sporadic, yet very high rates (up to 190 nmol L -1 d -1 ) of denitrification in this region, at that time much higher than the maximum observed anammox rate of 21 nmol L -1 d -1. Depending on the method of calculation, they arrived at a mean contribution of denitrification of 65-77 % to total mean removal of fixed nitrogen (2. 1 mmol N m -2 d -1 ), which is in striking contrast to previous studies indicating a pronounced dominance of the anammox process (Hamersley et al., 2007; Thamdrup et al., 2006; Galan et al., 2009; Kalvelage et al., 2011). As suggested by Dalsgaard et al. (2012) a possible explanation for the very different contributions of denitrification and anammox found in the different studies could be the spatial patchiness of the relevant processes in this region. 

In model MOPS presented here, we do not distinguish between denitrification and anammox, but assume only one process for the reduction of nitrate, denitrification, where both steps of nitrogen reduction take place at the same rate, in conjunction with immediate, complete oxidation of ammonium released during both steps (as in Eqs. 13-18 of Paulmier et al., 2009). Given the controversy related to the importance of anammox and canonical denitrification on global and regional scales and the rather simplistic way in which we implemented fixed nitrogen loss in the model, it seems worthwhile to have a closer look at the model's intrinsic assumptions and their validity with respect to net nitrogen fluxes. 

To examine the stoichiometric consequences of different pathways of fixed nitrogen loss, we first assume the complete denitrification pathway, i. e., complete oxidation of ammonium by nitrate, as suggested by Richards (1965), and applied by Paulmier et al. (2009, their Eq. 17). As an alternative pathway, we assume that ammonium released during the first step of denitrification is oxidized by anammox. As demonstrated in detail in Appendix A, with the assumed model stoichiometry both nitrogen loss pathways require 120 moles nitrate per 16 moles of ammonium (or 1 mole of organic phosphorus) oxidized. Thus, our current model stoichiometry can be regarded to represent either denitrification plus ammonium oxidation by nitrate, with both steps of denitrification proceeding at the same rate, or a combination of denitrification and anammox. 

Note that in the absence of any nitrite accumulation, the combination of denitrification and anammox implies that the second step of denitrification happens 1. 3 times faster than the first step (see Appendix A). In this case, the contribution of anammox to total dinitrogen production would only amount to ≈ 24 % (see also Fig. 13 ). If both steps of denitrification proceed at the same rate, we would again arrive at a low contribution of anammox of ≈ 26 %. A contribution of 24-26 % for anammox is close to the values found by Dalsgaard et al. (2012), and discussed by Koeve and Kähler (2010) and Ward (2013), but is far lower than suggested by some observations made in the Peruvian upwelling region. A recent analysis confirms a rather low contribution of anammox, possibly depending on the carbon : nitrogen stoichiometry of organic matter, its oxidation state, and the magnitude of organic matter supply (Babbin et al., 2014). 

Our simple theoretical framework suggests that anammox may, on average, only play a secondary role in determining the nitrogen loss in suboxic open-ocean areas. A more detailed model of different oxidation states of nitrogen, and other potential sources for ammonium, such as dissimilatory nitrate/nitrite reduction to ammonium (DNRA) (see also the extensive analysis and discussion by Koeve and Kähler, 2010) or zooplankton excretion, would be required in order to investigate these processes in conjunction with diffusive transport processes across the oxycline more closely. However, as shown above for biogeochemical tracer distributions simulated by a relatively coarse-resolution global model, it will not matter much whether the loss of fixed nitrogen is caused by denitrification or anammox, as both processes are ultimately fueled by organic matter and its remineralization products, with very small differences (or none at all) in the net stoichiometry and end products.
section_name: The importance of resolving different nitrogen species

section_content: Motivated by the study by Kalvelage et al. (2011), we assumed for our model setups BASE and DenHigh a wide tolerance of denitrification/anammox towards high levels of oxygen. However, recent studies by Dalsgaard et al. (2012) and De Brabandere et al. (2013) indicate that much lower oxygen concentrations are required for these processes to operate. The different observational estimates of oxygen thresholds of anammox were explained with regional differences between the study areas (De Brabandere et al., 2013). In an attempt to account for these uncertainties, we increased the oxygen affinity of aerobic remineralization and reduced the tolerance of denitrification to low oxygen in sensitivity experiment RemHigh. 

Considering the slow growth of anammox bacteria (e. g., Dalsgaard et al., 2012; Ward, 2013) and the consecutive steps of nitrogen reduction in denitrification, which in our model are all parameterized implicitly, together with their slightly reduced energy yield, we simulated inorganic nitrogen uptake via an increased half-saturation constant. For the same oxidant concentration, this parameterization will result in reduced rates of anaerobic remineralization compared to aerobic remineralization. Unfortunately, the affinities of denitrifiers to low nitrate concentrations are not well constrained. Low half-saturation constants for the nitrate uptake during denitrification of 2. 9 and 2. 5 mmol N m -3 were measured in the Mariager Fjord in northern Denmark (Jensen et al., 2009) and in the Gotland Basin (Dalsgaard et al., 2013), respectively. However, direct comparison of these observed values to our model setup and results is complicated for two reasons: First, environmental conditions in both studies were characterized by very low nitrite and low (usually < 5 mmol N m -3 ) nitrate concentrations. Further, the electron donor for denitrification was usually sulfide instead of organic matter. These conditions differ from many openocean or even coastal-ocean environments considered here. Second, with the addition of labeled nitrate, the observed half-saturation constants were as high as 31 mmol N m -3 for nitrate reduction, and at least 15 mmol N m -3 for denitrification measured via addition of labeled nitrite (Jensen et al., 2009), which has been explained with different substrate concentrations within and around the bacterial cells (Jensen et al., 2009). Given these methodological complications, and the regional differences between observations and our model setup, we are thus left with an uncertainty of an order of magnitude for the half-saturation constant for nitrate uptake during denitrification, ranging from 2. 5 to 31 mmol N m -3. This range of variation is to some extent addressed via our model experiments BASE, DenHigh, and RemHigh. It remains to be investigated, which half-saturation constant would be most appropriate for global simulations, and how far these would have to be changed when addressing more regional questions with more finely resolved models. Although the effects of these parameters seem to be negligible for the overall, global volume distribution of dissolved tracers (Fig. 4 ), or for global metrics (Fig. 11 ), they influence the simulation of OMZ volume, global nitrogen flux, and inventory (Fig. 10 ). Depending on questions posed to the model, it may therefore be important to constrain these parameters, either inversely, or via direct measurements particularly in open-ocean areas of denitrification.
section_name: Constraints for oxidant affinity of suboxic processes

section_content: Although some uncertainties in the parameterization of nitrogen losses and gains remain, on long timescales model MOPS generally match the global distribution of observed dissolved tracers and associated fluxes quite well. However, despite the quite different parameterizations of oxic and suboxic remineralization tested in this study, model experiments fail to represent all dissolved tracers simultaneously in the EEP (Fig. 6 ). The similar response of all model setups to changes in sinking speed suggests that processes other than biogeochemistry -most likely, the physical exchange between the different regions, as suggested by Dietze and Loeptien (2013) -play a role in determining the nitrogen budget in this region. Despite using circulation fields derived from a data-assimilative optimization of an ocean circulation model (Stammer et al., 2004), it is likely that our model circulation suffers from an imperfect representation of physical processes in this region, as has been hypothesized for many models by Dietze and Loeptien (2013). So far, neither parameterization of oxidant-sensitivity, nitrogen fixation, or sinking speed has helped to relieve the models from these errors in the EEP (likewise, in the Arabian Sea or Bay of Bengal). Local increases of the zonal isopycnal diffusivity can help to emulate the yet unresolved equatorial intermediate current system, and thereby improve the models (Getzlaff and Dietze, 2013). Given these possibly systematic model deficiencies, the spatial representation of suboxic zones and associated processes presented here must be viewed with caution. Based on the encouraging results of our model assessment against observed biogeochemical tracer distributions, we cautiously assume that the connection between areas of nitrogen fixation and nitrogen loss processes, and the associated transport timescales linking these regions in our model is representative of the real ocean.
section_name: The representation of nitrate and oxygen in the eastern equatorial Pacific

section_content: The above examples show that nitrate may, in models that explicitly account for denitrification, act as a kind of scapegoat, that incurs the results of deficiencies built into the physical or biogeochemical model formerly represented by oxygen, or by unspecified oxidants. Nitrate diagnosed from simulated phosphate in model BUR is slightly less sensitive to variations in sinking speed (Figs. 4, 6) ; however, its good fit to observed nitrate is only achieved with help of a strong dependency on observations. Additionally, given the infinite supply of (undefined) oxidants embedded implicitly in model BUR, this phosphorus-only model may not be well-suited to discriminate between the effects of biogeochemistry and physics on oxygen and other biogeochemical tracer distributions in this area. Therefore, the explicit simulation of the nitrogen cycle in global models may help to better constrain the dynamics of remineralization and its role in establishing and maintaining oxygen minimum zones. Variation of parameters that govern remineralization, nitrogen fixation and denitrification can have quite different effects on properties such as the extension of oxygen minimum zones, model nitrogen inventory, or global pelagic nitrogen fluxes. The different responses of these diagnostics to changes in either of these parameters are summarized in Fig. 12. 

The temperature constraint of nitrogen fixation has only a small effect on OMZ volume, and therefore, denitrification. It decreases nitrate inventory slightly, because of a smaller area in which fixed nitrogen loss can be compensated. A high particle sinking speed reduces the residence time of organic matter in the water column, thereby reducing aerobic remineralization, and the extent of OMZs. This then reduces global pelagic denitrification (because of a lower suboxic volume). Because of a lower importance of denitrification at higher sinking speeds, and because of the prescribed stoichiometry of N : P = 16 under oxic conditions, the global nitrogen inventory increases. A high nitrate affinity of denitrification reduces the size of the OMZ (because more nitrate, and less oxygen, is used for oxidation of organic material). However, the preference for nitrate also has the effect of increasing denitrification, thereby reducing the nitrate inventory. If additionally the affinity for oxygen is increased, both pathways (the aerobic and anaerobic remineralization) would be enhanced, which to some extent cancels out any effect of global nitrogen fluxes and inventory, despite a larger OMZ volume. 

To summarize, the extent of OMZ, nitrogen fluxes, and the nitrate inventory can be altered quite independently in the model (as is illustrated in Fig. 11 ), and therefore nitrogen and oxygen-based diagnostics provide useful additional constraints for the model. Together with the examination of the a priori assumptions of the model (see above subsections) we hope that in the future we will be able to make further progress towards a better constrained model of combined phosphorus, oxygen, and nitrogen cycles.
section_name: Can we simulate and constrain the nitrogen cycle in a global biogeochemical model?

section_content: Our model results are within the range of local, in situ observations of nitrogen fixation and denitrification. So far, we have refrained from using local observations of nitrogen fixation for model calibration, the reason being a quite sparse database, which exhibits a strong bias to certain oceanic regions (Luo et al., 2012). A closer examination of the effects of this bias on global and regional estimates, and the consequences for model calibration is beyond the scope of this study, and will be carried out elsewhere. As new data become available particularly in the so far undersampled, yet sensitive eastern tropical Pacific, these may provide very valuable benchmarks for model performance and skill. So far, we have not explicitly accounted for the iron cycle and its potential influence on cyanobacteria abundance and production. Beside observations of different iron species, its boundary exchanges through atmospheric dust deposition and benthic exchange (e. g., Nickelsen et al., 2015), improved data sets may also help to better constrain a more detailed model of interactions between these organisms, and different macro-and micronutrients. Likewise, we have, so far, refrained from simulating spatiotemporal variation in the N : P ratio of other, non-diazotrophic organisms. As shown by Martiny et al. (2013), the elemental stoichiometry of particulate organic matter can vary significantly between the different regions; it remains to be investigated how this will affect nutrient signals in the deep ocean, and the global nitrogen inventory on long temporal and large spatial scales, compared to the impact of nitrogen fixation and denitrification presented in this study. Data sets as presented in Martiny et al. (2014) will provide very useful constraints on future global model developments that include more flexible stoichiometries. 

Based on his high estimated nitrogen loss rates, Codispoti (2007) suggested that the, until then, generally lower estimates of oceanic nitrogen fixation might have been too low. More recent estimates of nitrogen fixation are often higher (Deutsch et al., 2007; Eugster and Gruber, 2012, 130-175 Tg N yr -1 ), while the maximum fixation rate simulated with MOPS is only 117 Tg N yr -1 (see Table 2 ). Because our model MOPS does not include benthic denitrification, and because nitrogen fixation is parameterized to balance nitrogen loss, its global integral is necessarily at the lower end of global estimates. For example, model results by Somes et al. (2013) point towards high rates of nitrogen fixation between 195 and 350 Tg N yr -1, sufficient to balance combined pelagic and sedimentary nitrogen loss, the latter being about twice as high as the pelagic loss. Including benthic denitrification in MOPS would most likely increase the global nitrogen fluxes, but also result in a different dependency of these on particle sinking speed. 

The sensitivity experiments with different remineralization kinetics have not revealed any dramatic changes in the simulated biogeochemical tracer distributions. On the other hand, experiments with different particle sinking speed indicate that not only the transient behavior of the models depends on this parameter, but also the evolution of suboxic zones, and thus steady-state fluxes of fixed nitrogen (nitrogen fixation, pelagic denitrification) as well as the simulated global nitrogen inventory and the nitrate-to-phosphate ratio. However, this is probably partly due to the fact that in some model regions fast-sinking organic matter quickly reaches the sediment, where it is ultimately buried. Considering benthic denitrification might shift this pattern, because in that case the oxidant deficiency at any given location cannot be neglected anymore. As noted above, this will, under the assumption of homeostasis, also affect the model's estimate of global nitrogen fixation. Including more flexible, benthicpelagic exchange processes, together with a dynamic iron cycle will therefore be useful and illustrative to the current model.
section_name: Further constraints and improvements of the simulated nitrogen cycle

section_content: We have carried out model simulations using global coupled biogeochemical ocean models that simulate phosphorus, oxygen, and nitrogen fluxes. Starting from global ob-served distributions of these tracers, our long-term simulations indicate that model inventories and fluxes exhibit considerable changes within the first few decades to centuries, particularly in the eastern tropical Pacific, but also globally. Global integrated fixed nitrogen sources and sinks converge to a steady state only slowly, on millennial timescales, and suggest that model results and trends achieved after spin-up periods shorter than a few thousand years should be viewed with caution. 

Compared to a model without nitrogen cycle, but with some form of implicit oxidants, dissolved nutrients and oxygen simulated by our new Model of Oceanic Pelagic Stoichiometry (MOPS) do not look very different, despite the fact that the latter imposes many more functional controls on biogeochemical fluxes. Although in all models we can produce an equally good fit to observed nitrate by multiplication of simulated phosphate with the observed global stoichiometric ratio, only the model with explicit nitrogen MOPS can predict this tracer prognostically, i. e., in the presence of more mechanistic, a priori assumptions, quite well. For this model observations of nitrate, its inventory and global flux can serve as a useful additional constraint. 

In MOPS nitrate replaces oxygen as oxidant in certain regions. Especially in the EEP nitrate exhibits the mismatch that phosphorus-only models show with respect to observed oxygen distributions. In our model simulations the EEP plays a large role not only for the initial transient of the model, but also for steady-state nitrogen fluxes. However, it is not clear how much of this response can be attributed to a deficient representation of the equatorial current system in this region. 

Stoichiometric considerations indicate that for global models on relatively coarse spatial grids and simulated over long timescales, it might not be necessary to differentiate between the various processes of fixed nitrogen loss, in particular denitrification and anammox. Correctly representing these processes is difficult without the explicit consideration of the various inorganic nitrogen species involved, and also depends on the availability of experimental data that determine and constrain the kinetics, substrates, and oxidants of the different processes. A review of the few observations of oxidant affinities reveals a wide range of these parameters, indicating the need for further research, especially given the difficulty of the current model metric to constrain these from other observations. 

With the given model setup, the effects of parameter variations on the extent of suboxic zones, nitrogen fluxes, and inventories differ between the different parameters, and suggest that these model diagnostics should be used for model skill assessment. Including these new constraints in the overall misfit, our results point towards a Martin exponent of 0. 86, and relatively low pelagic loss of fixed nitrogen between 59 and 84 Tg N yr -1, the latter supporting more recent, observation-based estimates. 

-kiel_biogeochem_model. F maps the 1-D tracer fields used by the TMM onto 1-D arrays used by the biogeochemical core routine BGC_MODEL. F, and back again afterwards. It calls the following:

-BGC_MODEL. F carries out the actual computation of biogeochemical sources and sinks presented here, including organic matter sinking and remineralization, air-sea gas exchange, and computation of carbon chemistry (option -DCARBON, see below). Thus, it is the heart of biogeochemistry. This routine requires daily average photosynthetically active solar radiation below sea surface, and day length. For this, we use a routine insolation. F provided by the MIT (http: //mitgcm. org/public/source_code. html), with some minor modifications by us. Any other forcing field for light can be provided. 

-kiel_biogeochem_diagnostics. F maps the diagnostic output (production, sedimentation,. . . ) computed in BGC_MODEL onto arrays to be passed to external_forcing_kiel_biogeochem -kiel_biogeochem_set_params. F is a dummy that may serve as a future module for changing parameters during optimization. 

Communication between the different modules is carried out mainly via header files:

-BGC_PARAMS. h is a header file that passes biogeochemical parameters between the different model pieces (from BGC_INI to BGC_MODEL). It also contains the biogeochemical tracer fields (bgc_tracer). 

-BGC_DIAGNOSTICS. h contains arrays for diagnostic output. 

-BGC_CONTROL. h is a header file that passes more technical runtime parameters to biogeochemistry, e. g., time step length, and vertical geometry. 

-kiel_biogeochem. h make subroutines known to external_forcing_kiel_biogeochem. c

A rather simple carbon module may be coupled to the Pcore via compile option -DCARBON. Note that these modules (CAR_CHEM. F, CAR_INI. F, CAR_PARAMS. h) are still somewhat preliminary, and will be presented in a later publication. They are largely based upon the routines developed and provided by MIT (http://mitgcm. org/public/ source_code. html) Appendix C: MOPS-1. 2 biogeochemical subroutines update MOPS-1. 0 has been updated with respect to code naming, structure, and user friendliness, resulting in a new version MOPS-1. 2. We note that this does not change the model results; however, we have experienced small (≈ O(10 -10 ) -O(10 -9 ) maximum difference for phosphate and oxygen, respectively, over the entire, spatiotemporal domain) changes because of transition to new hardware and a different PETSc-Version. Below is a brief documentation of the main changes made to the source code. The source code, together with the TMM driver code, forcing and scripts can be downloaded from the main TMM-website: https://github. com/ samarkhatiwala/tmm. We would like to refer the reader to that website, where the most recent version of the code, updates, as well as further documentation are available, and will be stored in the future. 

Main changes to the biogeochemical source code are as follows:

-For consistency, we have changed module names: all former name components kiel were replaced by mops. 

-external_forcing_mops_biogeochem. c now has a much more generic, user-friendly interface to add tracers during runtime. This also affects the modules called by this routine. 

-Default budget closure of phosphate and other tracers is now via supply at the sea surface; a switch for riverine supply (as used in this paper) is available via compile option -DRUNOFF. 

-BGC_INI. F now contains an additional compile option -DIMPRO for adding the parameterization of implicit profiles used by Kriest and Oschlies (2011). 

In the Supplement we provide a detailed list of changes, as well as some documentation about tracer changes due to different hardware, PETSc versions, and changes due to this new code version. 

The Supplement related to this article is available online at doi:10. 5194/gmd-8-2929-2015-supplement.
section_name: Conclusions

section_content: Acknowledgements. This work is a contribution to the DFGsupported project SFB754. Parallel supercomputing resources have been provided by the North-German Supercomputing Alliance (HLRN). The authors wish to acknowledge use of the Ferret program of NOAA's Pacific Marine Environmental Laboratory for analysis and graphics in this paper. Two referees provided very helpful reviews of the manuscript. We thank Samar Khatiwala for making the TMM source code available, and for his continuous, kind collaboration and support. 

Edited by: A. Ridgwell

section_name: Appendix A: Stoichiometry of denitrification and anammox

section_content: Consider 1 mole organic matter in phosphorus units (P org ), with the stoichiometric composition C a H b O c N d P. In our global model simulations we assume a composition of organic matter that requires 170 mole oxygen to oxidize 1 mole of organic phosphorus to carbon dioxide, water, phosphate, and nitrate. This value has been derived from geochemical observations (Anderson and Sarmiento, 1994), and has been applied in global biogeochemical models (e. g., Najjar et al., 2007; Moore and Doney, 2007). With d = 16, we define the oxygen demand for oxidation of organic matter to carbon dioxide, water, phosphate, and ammonium via R 0 = 170 -2 d = 138 (see also Paulmier et al., 2009). 

Under suboxic conditions, we first assume complete oxidation of ammonium by nitrate, as suggested by Richards (1965). Table 1 of Paulmier et al. (2009) indicates that for R 0 = 138 oxidation of 1 mole organic phosphorus requires 4/5 R 0 + 3/5 d = 120 moles nitrate, which is reduced completely to dinitrogen, without any surplus of nitrite. 

In an alternative approach we now consider anammox for oxidation of the ammonium released during suboxic degradation of organic matter, and assume that nitrite reduction during denitrification takes place at a rate x of nitrate reduction. We then arrive at the following, bulk stoichiometry for complete remineralization of 1 mole organic phosphorus:

(A1) Case x = 0, i. e., no nitrite reduction via denitrification, results in a considerable surplus of nitrite (2 R 0 -d = 260 mole nitrite for each mole of organic phosphorus remineralized; see Fig. 13, left panel), even though anammox consumes some of it. Even for x = 1 with the given stoichiometry 30 moles, nitrite would be generated per mole of remineralized organic phosphorus. The surplus of nitrite appears because during nitrate reduction the ratio of ammonium released from organic matter to nitrite produced from reduction of nitrate is not 1 : 1, as required for anammox. Only for x = (2 R 0 -d)/(4/3 R 0 + d) = 260/200 = 1. 3 no left-over nitrite would accumulate on the right hand side of equation A1. In this case, oxidation of 1 mole of organic phosphorus requires 4/5 R 0 + 3/5 d = 120 mole nitrate. 

To summarize, with the assumed model stoichiometry, and in the absence of any nitrite or ammonium accumulation, both cases require 120 mole nitrate per 16 moles of ammonium oxidized; i. e., our current model stoichiometry can be regarded to represent either denitrification plus ammonium oxidation by nitrate (with both steps of denitrification proceeding at the same rate) or a combination of denitrification and anammox. The latter case implies that the second step of denitrification happens 1. 3 times faster than the first step, to avoid any nitrite accumulation. 

Unfortunately, very little is known regarding the contribution of the different processes to nitrogen cycling in suboxic waters. Investigations in a Danish fjord rather suggest a dominance of nitrate reduction over nitrite reduction (Jensen et al., 2009) which, in our theoretical framework, would correspond to x < 1. However, given the quite unique hydrographical and biogeochemical conditions of that study (low nitrate and nitrite, sulfide as electron donor) it is not clear whether these findings can be transferred to our study, which focuses on the open ocean. Assuming no nitrite reduction by denitrifiers at all (x = 0) would result in a contribution of anammox of 100%, but would also result in a large surplus of nitrite, which does not seem to agree with observations. In the case x = 1. 3 (no leftover nitrite), the contribution of anammox to total dinitrogen production would only amount to ≈ 24 % (see also Fig. 13 ). If both steps of denitrification proceed at the same rate (x = 1), we would again arrive at a low contribution of anammox of ≈ 26 %.
section_name: Appendix A: Stoichiometry of denitrification and anammox

section_content: The biogeochemical source code of MOPS-1. 0 consists of outer routines (kiel_biogeochem_ini. F, kiel_biogeochem_model. F) that connect to the TMM and translate to the 3-D circulation, and inner routines that contain the actual biogeochemical sources and sinks, and define the biogeochemical parameters (BGC_MODEL. F, BGC_INI. F). They communicate via common blocks in header files BGC_PARAMS. h and BGC_CONTROL. h. Here we present the code that was used to generate plots provided in this paper; a updated version of this code with slightly changed structure, including the coupling to the transport matrix method (TMM) (Khatiwala et al., 2005) is available under https://github. com/samarkhatiwala/tmm. Its changes are presented briefly in Appendix C below, and in more detail in the supplementary material. 

external_forcing_kiel_biogeochem. c connects the outer biogeochemical subroutines to the TMM. It also reads the I/O files and runtime parameters. It calls the following subroutines:

-kiel_biogeochem_ini. F carries out some basic initialization, such as setting the time step length, initializing the tracer fields and vertical model structure, as well as some parts of the carbonate system (option -DCARBON, see below). It calls the following:

-BGC_INI. F, which sets the biogeochemical parameters (e. g., max. growth rate of phytoplankton) and may call CAR_INI. F to define the parameters for the carbon module (option -DCARBON, see below).
section_name: Appendix B: MOPS-1.0 biogeochemical subroutines