Contrasting juxtaposition of two paradigms for diazotrophy in an Earth System Model of intermediate complexity

Nitrogen fixers, or diazotrophs, play a key role in the carbon and nitrogen cycle of the world oceans, but the controlling mechanisms are not comprehensively understood yet. The present study compares two paradigms on the ecological niche of diazotrophs in an Earth System Model (ESM). In our standard model configuration, which is representative for most of the state-of-the-art pelagic ecosystem models, diazotrophs take advantage of zooplankton featuring a lower food preference for diazotrophs than for ordinary phytoplankton. We compare this paradigm with the idea that diazotrophs are more competitive 5 under oligotrophic conditions, characterized by low (dissolved, particulate, organic and inorganic) phosphorous availability. Both paradigms are supported by observational evidence and lead to a similar good agreement to the most recent and advanced observation-based nitrogen fixation estimate in our ESM framework. Further, we illustrate that the similarity between the two paradigms breaks in a RCP 8.5 anthropogenic emission scenario. We conclude that a more advanced understanding of the ecological niche of diazotrophs is mandatory for assessing the cycling of essential nutrients, especially under changing envi10 ronmental conditions. Our results call for more in-situ measurements of cyanobacteria biomass if major controls of nitrogen fixation in the oceans are to be dissected.


Introduction
Nitrogen is an essential element in the metabolism of plants. Even though it is much more abundant in air than for example carbon dioxide, the assimilation of carbon by phytoplankton in vast regions of the ocean is considered to be limited by the 15 availability of nitrogen. Among the reasons for this is that most phytoplankton cannot use the molecular nitrogen, that is so abundant in air, because it cannot break the exceptionally strong chemical bond between the two nitrogen atoms constituting molecular nitrogen (N 2 ). An exception to this rule are so-called nitrogen fixing microorganisms (bacteria and archaea) that are capable of breaking this bond and thus can utilize the molecular nitrogen of the air. Their total input of bioavailable nitrogen to the ocean is substantial: current estimates range from 70 to 200 T gN yr −1 (c.f., Großkopf et al., 2012;Gruber 20 and Sarmiento, 1997; Tang et al., 2019b;Wang et al., 2019). Paleo-records suggest that for the last several thousand years, this input has been balanced by denitrification and anaerobic ammonium oxidation (anammox) which converts bioavailable forms of N (e.g. NO − 3 , NO − 2 ) under low-oxygen conditions back to N 2 (Altabet et al., 2012). The apparent balance of sources and sinks of bioavailable N suggests a coupling mechanism between these sources and sinks. Such coupling mechanisms and related controls on the overall oceanic N-budget are discussed controversially (e.g., Gruber, 2016;Wang et al., 2019) -also 25 because the controls on nitrogen fixation are not comprehensively understood.
There is, however, consensus, that anthropogenic forcing, such as global warming and input of bioavailable nitrogen to the ocean, modulates the turnover of bioavailable nitrogen in the ocean which might have far-reaching consequences (Duce et al., 2008;Krishnamurthy et al., 2007;Krishnamurthy. et al., 2009;McMahon et al., 2015). Some studies suspect that N 2fixation might increase regionally in the future, adding extra nutrients to already-overfertilized systems (e.g., 30 Wasmund et al., 2001), or even result in vicious cycles (e.g., Vahtera et al., 2007). This comprehension has triggered efforts to capture the respective dynamics in numerical models which explicitly resolve oceanic transport mechanisms (such as currents and mixing) along with the turnover of biogeochemical species (nutrients, carbon and oxygen) -with the aim to forecast the effects of changing environmental conditions and thereby facilitating management and mitigation measures.
A generic approach to capture the major aspects of the marine N-cycle in numerical models is to include an explicit rep-35 resentation of diazotrophs as an extra functional group, that is distinctly different to the representation of ordinary non-fixing phytoplankton (e.g., Fennel et al., 2001;Moore et al., 2001;Paulsen et al., 2017). The underlying model assumptions and respective mathematical formulations vary from one model to another, but they typically agree on (1) diazotrophs are capable of utilizing N 2 , while ordinary phytoplankton is not and (2) diazotrophs grow slower (or at least not faster) than ordinary phytoplankton, because of the metabolic cost involved in providing means to break the strong bond between the two nitrogen atoms 40 constituting N 2 . When taken alone, a direct deduction of these two assumptions is that diazotrophs only have an advantage over ordinary phytoplankton in regions where bioavailable nitrogen supply is below the demands of ordinary phytoplankton, relative to the phosphorus supply (Deutsch et al., 2007;Schmittner et al., 2008). In the remainder of the ocean, ordinary phytoplankton would suppress the relatively slow growing diazotrophs by competing more successfully for resources, such as bioavailable phosphorus (P), iron (for which diazotrophs have a higher demand) and light (e.g., Falcón et al., 2005;LaRoche and Breitbarth, with the UVic 2.9 Earth System Model (Weaver et al., 2001). The major aim is to illustrate the far reaching implications of both paradigms in a generic model that is currently used to project into our warming future and to assess geoengineering options (see, e.g., Keller et al., 2014;Kemena et al., 2019;Mengis et al., 2016;Reith et al., 2016). Our experiments are designed to constrain the envelope of model responses as set by each of the paradigms. We find (and illustrate in this study) that the development of a reliable model must be preceded by additional in-situ observations.

75
We use the UVic 2.9 Earth System Model, documented and extensively assessed in Keller et al. (2012). We refer to this standard model version since it is extensively used (e.g., Keller et al., 2014;Kemena et al., 2019;Reith et al., 2016), well documented (Keller et al., 2012Schmittner et al., 2008) and was used earlier to study diazotrophy (e.g., Landolfi et al. (2017); Somes and Oschlies (2015)). The model is of intermediate complexity and thus allows for running a multitude of simulations to quasi-equilibrium. The horizontal resolution is 1.8 • in latitude and 3.6 • in longitude in all submodules (ocean, including 80 biogeochemistry, dynamic-thermodynamic sea ice module, atmospheric energy-moisture balance model, land module, with an active terrestrial vegetation component, and simple land ice). The ocean submodule is based on a three-dimensional primitiveequation model (Pacanowski, 1995). The vertical discretization of the ocean comprises 19 levels and increases gradually from 50 m at the surface to 500 m at depth. The vertical background mixing parameter, κ h , is 0.15 cm 2 s −1 throughout the water column. In the Southern Ocean (south of 40 • S) this background value, κ h , is increased by 1.0 cm 2 s −1 . An anisotropic viscosity 85 scheme, suggested by Large et al. (2001), is implemented to improve the equatorial circulation (Somes et al., 2010). A marine pelagic biogeochemical model is coupled to the ocean circulation component. Its prognostic variables are ordinary (i.e. nondiazotrophic) phytoplankton (P o ), diazotrophic phytoplankton (P d ), zooplankton (Z), detritus (D), nitrate (NO 3 ), phosphate (PO 4 ), dissolved oxygen (O 2 ), dissolved inorganic carbon (DIN) and alkalinity (ALK). The model is described in detail in Keller et al. (2012) and Schmittner et al. (2008). Note that the dispersive numerics of the oceanic advection occasionally produces small negative values in phytoplankton and zooplankton concentrations. These values are set to zero in the model evaluation.
In the following, we present a choice of details relevant for the simulated dynamics of diazotrophs (relevant model parameters are listed in Tab.1): phytoplankton growth is generally controlled by the availability of light and nutrients (here, nitrate, phosphate and iron, where the latter is parameterized rather than explicitly resolved). Phytoplankton blooms are terminated 95 by zooplankton grazing, once the essential nutrients are depleted and the related reduction in phytoplankton growth does not longer keep up with the grazing pressure, that has built up during the bloom. For ordinary phytoplankton the maximum potential growth rate is where k F e =0.1 mmol F e/m 3 is the half-saturation constant for iron (Fe)-limitation and a = 0.6 d −1 determines the maximum potential growth rate at an ambient temperature T = 0 • C.

100
For diazotrophs the formulation of the maximum potential growth is similar, but diazotrophs are disadvantaged in nitratereplete waters by C d = 0.4. In addition, growth is stalled in waters colder than 15 • C: The actual growth rate of non-diazotrophic phytoplankton, J O , is, in case of low irradiance (I) and/or nutrient depleted conditions, the maximum potential rate J max o reduced by the following implementation of Liebig's law of the minimum: where k N and k P are the half saturation constants for N O 3 and P O 4 , respectively.

105
The actual growth rate of diazotrophs, J D , is similar but it is not affected by nitrate deficiency: The light-limited growth J IO and J ID are both determined by: The initial slope of the P-I curve (α) is set to 0.1 (W/m 2 ) −1 d −1 . Diazotrophs take up available nitrogen following: where P D denotes the biomass of diazotrophs (in units molN/m 3 ). When no or not enough bioavailable N is available, nitrogen fixation tops the respective N pool up to a Redfield-ratio of N:P=16. Note that the above formulation differs slightly from the Our minor change ensures a more realistic behavior where no nitrogen is fixed under nitrate-replete conditions. When applied to the reference model version, the difference to the original Keller et al. (2012) simulation turned out to be negligible.
with a maximum growth rate µ max = 0.4 d −1 . Grazing on ordinary phytoplankton, P O , is calculated by setting θ = θ o = 0.3.
Grazing on diazotrophs, P D , is calculated with a lower grazing preference θ = θ d = 0.1.

Implementation of the two Paradigms
We integrate two sets of numerical model configurations. The first set explores the paradigm selective grazing, represented by the model setup GRAZ and assumes that diazotrophs are grazed less than ordinary phytoplankton. The second set explores 120 the paradigm low P-demands, represented by the model setup OLIGO. We consider both paradigms individually to constrain the envelope of model responses as set by each of the paradigms. To ensure comparability between the paradigms we identify parameters for both paradigms which represent available observations in an "optimal" way (see below). Finally we compare the respective "optimal" representatives of the two paradigms, OLIGO and GRAZ, in Section 3.
The selective grazing paradigm is implemented by adjusting the grazing preferences for diazotrophs and ordinary phyto- More specifically we calculate the sum of absolute deviations between the locally estimated nitrogen fixation of Wang et al.

130
(2019) and our simulated nitrogen fixation (mean absolute error): N m denotes the number of estimated values for nitrogen fixation (on the model grid), a j stands for an 'observation' at location j andâ j for the corresponding model result. We consider climatological annual means. Note that this choice of misfit metrics inevitably adds a subjective element and the respective choice will impact the results of the optimization procedure (see Dietze, 2015, 2017). A choice has to be made, however, and here we follow a standard approach (c.f., Pontius et al., 2008;135 Sauerland et al., 2018;Willmott and Matsuura, 2005). Note that solutions with very low biomass of diazotrophs (below an integrated value of 10 T gC) are not considered.
The second paradigm, represented by OLIGO, assumes that diazotrophs can cope better with oligotrophic conditions than ordinary phytoplankton. This paradigm is implemented by (1)  In addition, we will show results from simulations of GRAZ and OLIGO, covering the period 1800-2150. They start with the historical state and are continued with the CO 2 emission scenario RCP 8.5 (Riahi et al., 2011).

155
In situ observations of diazotrophic biomass and nitrogen fixation rates, such as compiled by Luo et al. (2012), are very useful but still too sparse to allow for robust comparisons between model and observations. This has been a serious drawback hindering the identification of the major controls of diazotrophy. This study exploits the most recent and most comprehensive product from Wang et al. (2019), that fills data gaps by combining observations with inverse biogeochemical and prognostic ocean modelling -an approach pioneered by Deutsch et al. (2007).

Results
Our tuning exercises for the two paradigms selective grazing and low P-demands identified two respective "best" simulations out of two sets of 24 model integrations, which, according to Equation 9, feature the highest similarity to the most recent estimate of global oceanic nitrogen fixation (setups GRAZ and OLIGO; c.f., Sect.2.2). In the following, we present the respective model results with respect to diazotrophy under historic climate conditions. These results are put into perspective to the 165 reference model version REF. We continue with the results from RCP 8.5 scenario integrations.
Both "best" representatives of the two considered paradigms, GRAZ and OLIGO, feature a very similar fit to the observational estimate for nitrogen fixation, clearly superior to the reference simulation REF (the respective costs are 18.2 and 17.6, versus 27.7 mmol N/m 2 /yr). A comparison of our simulations with the, unfortunately very sparse, but independent (in the sense that this data has not been used to identify our the "best" simulations with Equation 9) data collection of diazotrophic 170 biomass from Luo et al. (2012), indicates that the simulated distribution of diazotrophs based on OLIGO is more realistic than in GRAZ. The biomass of diazotrophs in GRAZ is apparently too low. When taking the upper 10% percentiles as a measure for bloom intensities (derived from the respective histograms of positive local annual mean values), OLIGO lies with 716 mg C/m 2 much closer to the observations (598 mg C/m 2 ) than GRAZ (282 mg C/m 2 ) and REF (320 mg C/m 2 ). Please note, however, that the model data are climatological while the observations are rather anecdotal such that the comparison 175 relies on the implicit assumption that, e.g., seasonal variations are small (also assumed by, e.g., Landolfi et al., 2015). Even though OLIGO and GRAZ are very similar, some smaller differences among the two are evident: while GRAZ features generally rather low fixation rates in the Southern Hemisphere, this model version seems most realistic in the North Atlantic.

185
With respect to diazotroph biomass, the differences between the setups are more pronounced. While the simulated pattern in nitrogen fixation roughly matches the distribution of diazotrophs in REF (Fig.1a, compared to Fig.2a) and GRAZ (Fig.1b, compared to Fig.2b), this relation is not as evident for OLIGO: e.g., the amount of fixed nitrogen in the western tropical Pacific is rather low compared to abundance of diazotrophs (Fig.1c, compared to Fig.2c). The major reason why the overall pattern in dizotroph abundance and fixation do not fully match in OLIGO is that diazotrohs find an ecological niche in oligotrophic 190 regions (as they require less P than ordinary phytoplankton) without the need to fix atmospheric nitrogen. As desert dwellers they occupy more area but, even so, they can not necessarily sustain high production rates since access to essential bioavailable phosphorous is limited. Correspondingly, the quota of annual mean fixation rates versus the biomass of diazotrophs is in OLIGO more than 3 smaller than in REF and factor 2 smaller than in GRAZ. Also, in contrast to REF and GRAZ, minimum values even go down to zero for OLIGO. Note that these quotas were based on locations with a diazotrophic biomass of more 195 than 5 mgC/m 2 only (to avoid division by very small numbers).
In a next step, we determine the parameter sensitivities for the two paradigms, selective grazing and low P-demands (which optimal solutions are represented by GRAZ and OLIGO, resp.). Specifically, we explore the sensitivity of changes in the half saturation constant, k d P , versus changes in the food preference, θ d (based on the respective ensembles of simulations described in Sect.2.2, while C d is kept constant on the respective optimal level). Changing both parameters, θ d and k d P , by 20% of their 200 covered range (0.02 and 0.0088, resp.) leads on average to a much larger change in the amount of fixed nitrogen for the selective grazing paradigm, compared to low P-demands (16.8 and 0.8 T g N yr −1 , respectively). For the integrated biomass of diazotrophs however, the respective changes are similar for both paradigms (4.9 T g C for selective grazing v.s. 5.8 T g C for low P-demands). Note, in this context that for both paradigms, the biomass of diazotrophs reacts strongly nonlinear to parametric changes. Interestingly, decreasing the food preference for diazotrophs can even lead to a tipping point where the 205 spatial distribution of diazotrophs is shifted entirely to the upwelling regions.
We conclude that the paradigm low P-demands produces a model behavior that is much more robust towards parameter changes than the selective grazing paradigm, in the sense that small changes in (rather unconstrained) model parameters effect moderate changes in the quality of the fit to observations. In contrast, small changes to the selective grazing paradigm might 210 even induce a regime shift.
A question which arises in this context, is whether the relatively large parameter sensitivity of the paradigm selective grazing is also reflected in the model sensitivities towards environmental changes. To explore the effects of changing climate conditions, we projected the model versions OLIGO and GRAZ into a warming future, corresponding to the RCP 8.5 emission scenario.
The response of both, OLIGO and GRAZ, are closely tied to an increasing vertical stratification effected by an ocean, that 215 is warmed from above. The water expands at the surface and this increases the vertical density gradient, such that the energy requirements for vertical mixing are increased, because mixing is now associated with pushing lighter, more buoyant water downwards. The increased energy demands for mixing result in a dampening of mixing events and, overall, in less nutrients mixed upwards to the sun-lit surface. Among the processes setting in are (1)

225
Because of their underlying paradigm, OLIGO and GRAZ respond very differently to this increase in oligotrophic regionseven though they share a very similar behavior under historical conditions: in the simulation OLIGO, the diazotrophs can take advantage of the increasing vertical stratification because the diazotrophs are expertly exploiting oligotrophic conditions. The black line in Fig.3a shows that the biomass of diazotrophs increases globally along with an increase in size of the oligotrophic regions. In addition, Fig.3b illustrates that, even though the size of oligotrophic regions and the diazotrophic biomass increases, 230 the nitrogen fixation decreases in OLIGO. The reason is a reduced total supply of phosphorous to the surface in (expanded) oligotrophic regions which, in our model, puts an upper limit on nitrogen fixation. Additionally, nitrogen fixers have in the extended oligothrophic regions a competitive advantage over ordinary phytoplankton, because they require less P for their subsistence. Subsistence is, however, not necessarily correlated with high production rates which, ultimately, are limited by the availability of phosphorous -whose limitation provided the niche for diazotrophs in OLIGO.

235
In contrast, the development of nitrogen fixation and abundance of diazotrophs goes much more hand-in-hand in GRAZ and the projected evolution of both variables differs strongly from OLIGO: an initial increase in both variables is followed by a subsequent decay. Such switching behavior at arbitrary tipping points (that are set by the model parameters) are typical for the rather generic zooplankton formulation used here and have been described, e.g., in Löptien and Dietze (2017), their Sect.3.2.
Note, however, that for both model versions the Bay of Bengal features a sudden pronounced onset of nitrogen fixation in the 240 mid 2000s, which is visible even in the global mean. We anticipate that this abrupt change might be to a particularly pronounced sea surface temperature increase in this region.
We have been exploring two paradigms that are proposed in the literature to essentially define the environmental conditions or niches where oceanic diazotrophs may thrive -even though they typically grow slower than ordinary phytoplankton, with 245 which they compete for essential resources such as phosphate. The first paradigm, a de-facto standard in global earth system modelling, builds on selective grazing, where the diazotrophs are exposed to a grazing pressure that is reduced relative to other, ordinary phytoplankton.  . Our experiments also show that the similarity of the model versions based on the differing paradigms breaks in projections into a warming future. This suggests that current observational data are not sufficient to constrain the sensitivity of the dynamics of bioavailable nitrogen under changing environmental conditions. We find that for assessing the fidelity of respective paradigm additional observations of diazotrophic biomass are needed most. To 255 this end, the quota of nitrogen fixation versus the respective biomass of diazotrophs may prove especially helpful.
Data availability. The model output is preliminary archived at https://nextcloud.ifg.uni-kiel.de/index.php/s/iYx7CWaGN8rx8KE and will be provided in the Institute repository after review. The observational data collected by Luo et al. (2012) are available via https://doi.pangaea.de/10.1594/PANGAEA.774851 and cited in the references.
Author contributions. Both authors were involved in the design of the work, in data analysis, in data interpretation and in drafting the article.