Modelling the biogeochemical effects of heterotrophic and autotrophic N2 fixation in the Gulf of Aqaba (Israel), Red Sea

Recent studies demonstrate that marine N2 fixation can be carried out without light by heterotrophic N2-fixers (diazotrophs). However, direct measurements of N2 fixation in aphotic environments are relatively scarce. Heterotrophic, as well as unicellular and colonial photoautotrophic diazotrophs, are present in the oligotrophic Gulf of Aqaba (northern Red Sea). This study evaluates the relative importance of these different diazotrophs by combining biogeochemical models with time series measurements at a 700m-deep monitoring station in the Gulf of Aqaba. At this location, an excess of nitrate, relative 15 to phosphate, is present throughout most of the water column and especially in deep waters during stratified conditions. A relative excess of phosphate occurs only at the water surface during nutrient-starved conditions in summer. We show that a model without N2 fixation can replicate the observed surface chlorophyll but fails to accurately simulate inorganic nutrient concentrations throughout the water column. Models with N2 fixation improve simulated deep nitrate by enriching sinking organic matter in nitrogen, suggesting that N2 fixation is necessary to explain the observations. The observed vertical structure 20 of nutrient ratios and oxygen is reproduced best with a model that includes heterotrophic, colonial and unicellular autotrophic diazotrophs. These results suggest that heterotrophic N2 fixation contributes to the observed excess nitrogen in deep water at this location. If heterotrophic diazotrophs are generally present in oligotrophic ocean regions, their consideration would increase current estimates of global N2 fixation and may require explicit representation in large-scale models.


Introduction 25
Biological nitrogen fixation refers to the conversion of dinitrogen gas (N2) via reduction to ammonium (NH4) into bioavailable forms of nitrogen by a specialized group of microbes containing the nitrogenase enzyme complex. On geological timescales, the size of the oceanic reservoir of bioavailable nitrogen, and thus the ocean's capacity for exporting carbon to depth, is controlled by the balance between removal of fixed nitrogen by denitrification and input by N2 fixation (Falkowski, 1997;Haug et al., 1998;Deutsch et al., 2007;Gruber and Galloway, 2008;Fennel et al., 2005). The amount of organic matter 30 exported from the surface to the deep ocean (i.e., export production) depends on allochthonous inputs of nitrogen (i.e., "new nitrogen") into the euphotic zone (Eppley and Peterson, 1979). These new nitrogen inputs determine the amount of "new production", which is directly related to the exported fraction. Locally the supply of new nitrogen can occur through several mechanisms, including microbially mediated N2 fixation, diapycnal mixing injecting deep nitrate (NO3) into the surface, lateral transport, atmospheric sources and riverine input. While the injection of deep NO3 is often regarded as the dominant source of new nitrogen that drives the seasonal cycle of marine primary production; there is significant interest in quantifying the 5 contribution of N2 fixation to primary production, particularly in oligotrophic areas (Karl, 2002;Zehr and Ward, 2002;Capone et al., 2005;Luo et al., 2012).
Diazotrophs contain nitrogenase, the catalyst enzyme for N2 reduction, which is encoded by nif genes. Trichodesmium spp., a group of non-heterocystous filamentous cyanobacteria that form large colonies, were traditionally considered the main contributors to N2 fixation in the surface subtropical and tropical ocean (Carpenter and McCarthy, 1975;Capone et al., 2005). 10 Increased sampling efforts and method improvements subsequently led to the discovery of diverse diazotroph groups including heterocystous endosymbiotic cyanobacteria (Zehr et al., 1998;Carpenter et al., 1999), free-living unicellular cyanobacteria (Zehr et al., 2001;Montoya, 2004;Moisander et al., 2010), and other cyanobacterial symbionts (Zehr et al., 2000). Most recently, genetic techniques have allowed the detection of nif genes in a number of anaerobic and heterotrophic phylotypes (Zehr et al., 2008;Zehr, 2011;Rahav et al., 2013Rahav et al., , 2015. The abundance of nif genes does not necessarily imply that these 15 organisms are actively fixing N2 (Zehr et al., 2000;Moisander et al., 2017); however, the correlation between bacterial productivity and N2 fixation rates suggests that significant aphotic N2 fixation may occur in the Red Sea (Rahav et al., 2013(Rahav et al., , 2015. The differences in size and physiology of these diverse diazotrophs also suggest that they occupy distinct niches, and thus may affect primary productivity and export production differently (Bonnet et al., 2016;Moisander et al., 2010).
Most biogeochemical models treat N2 fixation as a purely light-dependent, autotrophic process. These models either use 20 standard formulations of light limitation for the diazotrophic groups (e.g., Fennel et al., 2002;Moore, et al., 2004;Gregg, 2008;Dutkiewicz et al., 2012), or include theoretical considerations to introduce an empirical N2 fixation flux in the model (e.g., Bisset et al., 1999). Some approaches neglect light limitation on diazotrophy and instead infer global N2 fixation patterns from the distribution of dissolved inorganic nitrogen and phosphorus, and estimates of ocean circulation (e.g., Deutsch et al., 2007).
In trait-based models, diazotrophs are usually accounted for by a single functional group with biological rate parameters 25 intended to represent either Trichodesmium spp., unicellular cyanobacteria, or a generic autotrophic diazotroph. More complex models of diazotroph growth have explored cellular function and the effects of a variable cellular stoichiometry, however still limited by light (Kreuss et al., 2015;Fernández-Castro et al., 2016). Only a few modelling studies have evaluated multiple autotrophic diazotroph groups simultaneously, by considering separate groups for Trichodesmium spp., unicellular cyanobacteria and diatoms-cyanobacterial associations (e.g., Monteiro et al., 2010;Dutkiewicz et al., 2012Dutkiewicz et al., , 2015. To our 30 knowledge, heterotrophic N2 fixation has not yet been considered explicitly in biogeochemical models. Understanding the ecological dynamics of different types of diazotrophs should significantly improve predictive capabilities in biogeochemical models, and lead to more accurate estimates of global N2 fixation rates. It has been suggested that N2 fixation rates are underestimated globally due to limited knowledge about the distribution and characteristics of N2 fixing organisms (Montoya, 2004;Zehr, 2011). It is also assumed that marine N2 fixation may increase globally, as a result of ocean warming and higher concentrations of dissolved CO2 in sea water (Hutchins et al., 2007;Levitan et al., 2007;Dutkiewicz et al., 2015). 5 To our knowledge, these laboratory experiments have only explored the response of Trichodesmium. Less information is available about the effects of climate trends on other diazotrophic organisms.
In this study we explore the biogeochemical signatures that result from different assumptions about the ecological niches occupied by diazotrophs. Our study area is the Gulf of Aqaba, a northern extension of the Red Sea. Aside from the reported presence of diverse diazotroph types, the morphology of the Gulf of Aqaba limits horizontal transport of deep waters, thus 10 allowing us to simplify the physical model's complexity and focus on the biological component. We aim to answer the following two questions: i) How important is N2 fixation as a source of new nitrogen in the Gulf of Aqaba? ii) How important is heterotrophic, light-independent N2 fixation? To address these questions, we implemented a one-dimensional model at a monitoring station for which monthly quality-controlled measurements of physical and biogeochemical variables are available from 2004 onward. We then systematically tested different model assumptions about diazotrophy and calibrated selected model 15 parameters to facilitate an objective comparison between the different biogeochemical model versions. The different assumptions about diazotrophy consider the characteristics of organisms identified in the Gulf of Aqaba, including heterotrophic fixing α and γ proteobacteria (Rahav et al., 2013), unicellular cyanobacteria, and Trichodesmium spp. (Post et al., 2002;Foster et al., 2009;Rubin et al., 2012). Our most important conclusion is that aphotic N2 fixation is necessary to reproduce the observed excess nitrogen in deep waters of the Gulf, while maintaining reasonable surface N2 fixation rates. Our 20 best performing biogeochemical model of the Gulf of Aqaba estimates annual N2 fixation rates in overall agreement with local and large-scale estimates in the literature.

Study Area: The Gulf of Aqaba
The Gulf of Aqaba is a quasi-rectangular, 200-km long, 20-km wide, semi-enclosed basin in the northeast region of the Red Sea ( Figure 1). The average depth of the Gulf of Aqaba is 800 m, its deepest point approximately 1800 m, and it is surrounded 25 by arid mountains that steer the dominantly northerly winds (Berman et al., 2003). Two shallow sills, the Bab el Mandeb (~140 m) and the Strait of Tiran (~240 m), inhibit the entrance of cold and dense deep waters from the Indian Ocean. Since inflow is restricted to warm surface waters, the Gulf's deep-water masses (>300 m) are locally formed (Wolf-Vecht et al., 1992;Biton et al., 2008) and have negligible horizontal transport toward the exterior (Klinker et al., 1976;Manasrah et al., 2006).
The annual hydrographic cycle exhibits a well-defined seasonality where vertical temperature and salinity distributions are 30 dominantly affected by surface heat fluxes and modified by surface advective fluxes (Carlson et al., 2014). During winter (December to March), convective vertical mixing usually extends to depths >300 m (Labiosa et al., 2003), and even reaches the bottom (~700-800 m bottom depth) in some extreme years ( Figure 2). From April to September, the water column is thermally stratified, and inflowing warm surface waters from outside the Gulf occupy the layer above the thermocline (Genin and Paldor, 1998;Berman et al., 2000;Biton and Gildor 2011). During fall (October to December), surface cooling and high evaporation rates erode the seasonal stratification and re-establish a well-mixed water column (Berman et al., 2003;Monismith 5 and Genin, 2004). The Gulf experiences net evaporation of approximately 1.6 m yr -1 (Ben-Sasson et al., 2009) due to negligible precipitation and run-off (Wolf-Vecht et al., 1992).
The Gulf is oligotrophic, with surface nitrate (NO3) and phosphate (PO4) concentrations usually close to their detection limits during summer stratification (Fuller et al., 2005;Mackey et al., 2009;Meeder et al., 2012). Deep winter mixing supplies inorganic nutrients to the surface, and NO3 and PO4 reach ~0.1 µM and ~2 µM, respectively (Figure 2;Lindell and Post 1995;10 Lazar et al. 2008). Dust from the desert provides a sufficient atmospheric source of soluble iron (Fe) for microbial growth in the Gulf (Chase et al., 2006;Chen et al., 2007). Phytoplankton spring blooms reach a maximum chlorophyll concentration of ~2 mg Chl-a m -3 and their initiation is strongly correlated with the termination of winter cooling of the water column (Zarubin et al., 2017). Interannual variability in the depth of winter convective mixing results in periods of nutrient accumulation in deep waters (Figure 2;Wolf-Vecht et al., 1992;Lazar et al., 2008;Carlson et al., 2012), which is re-set during extreme winter 15 mixing events approximately every four years . The periodicity of these extreme mixing events has been associated with regional weather patterns that modify the Red Sea water temperatures .

Methods
We analysed the role of autotrophic and heterotrophic N2-fixing organisms in determining biogeochemical patterns at an open pelagic site (Station A), located in the northern Gulf of Aqaba, by testing four alternative ecosystem model versions. The 20 ecosystem models are evaluated in terms of their ability to replicate observations of oxygen (O2), NO3, PO4, and chlorophyll.
In this section, we first describe the available observations, then the models, and finally the systematic model calibration method.

Observations
Meteorological and oceanographic observations are available from the Inter-University Institute (IUI) for Marine Sciences in 25 Eilat, Israel (http://www.iui-eilat.ac.il/Research/ NMPMeteoData.aspx). Meteorological observations are used to calculate surface heat and momentum fluxes for the physical model and incoming light for the biological models. Observed meteorological variables include wind speed, air temperature, air humidity, air pressure, irradiance and cloud cover. This data has been collected continuously and automatically at 10 min intervals by the meteorological instrumental array at the end of the IUI pier since 2006. 30 Monthly CTD and bio-chemical profiles at Station A (29.5° N, 34.9° E) were collected during monthly surveys of the National Monitoring Program (NMP) from 2004 to 2014 (http://www.iui-eilat.ac.il/Research/NMPMeteoData.aspx). CTD profiles are used to nudge temperature and salinity in the physical model (see section 3.2). Bio-chemical profiles, including NO3, nitrite (NO2), ammonium (NH4), PO4, O2 and chlorophyll-a (Chl-a), are used for biogeochemical model calibration (years 2006 to 2010) and model validation against unassimilated data (years 2011 to 2014). Nutrients were measured using spectrophotometry 5 (QuickChem 8000 flow injection), O2 was determined by Winkler titrations, and Chl-a concentrations are estimated using fluorometry (Turner Designs 10-AU).

Model Descriptions
The ecosystem models are implemented within the General Ocean Turbulence Model (GOTM), a one-dimensional physical model that computes solutions to differential equations for the vertical transport of momentum, salt and heat (Burchard et al., 10 1999  phytoplankton, zooplankton and detritus are initially set to a homogeneous small value of 0.1 mmol N m -3 , and also readjust rapidly because the adjustment timescales for these variables are short (Fennel et al., 2006). Diazotrophs initial values are set to lower densities than non-fixing organisms, with a homogenous total value of 0.03 mmol N m -3 (i.e., models with multiple diazotrophs maintain the same amount of initial diazotrophic biomass by dividing initial values equally among organisms). Four main ecosystem model versions of increasing complexity (referred to as H0, H1, H2 and H3) are treated as alternative 25 hypotheses of how biological processes, especially diazotrophy, control the vertical distribution and temporal variability of dissolved inorganic nutrients and oxygen. H0 is the base model without diazotrophics and follows the model equations described in Fennel et al., (2006Fennel et al., ( , 2013. In general, the model follows Monod kinetics using a fixed N:P ratio ( ":$ %& = 16). Sensitivity to the constant N:P ratio is explored in section 4.2.3. We test the H0 model with and without a sediment denitrification flux (model versions H0 and H0', respectively). H0 includes denitrification but no N2 fixation, as no diazotrophs 30 are considered. H0' does include neither denitrification, nor diazotrophic organisms, and thus the underlying assumption of this model version is that there is a balance between inputs from N2 fixation and losses of fixed nitrogen due to denitrification.
When present, the denitrification flux follows Fennel et al. (2013) with a loss fraction of 6 mol N2 per mol of organic matter remineralized at the sediment-water interface. This generates an average sediment denitrification flux of 0.25 ± 0.46 mmol N m -2 d -1 , with a maximum value of 3.01 mmol N m -2 d -1 .
H1, H2 and H3 are modified versions of H0, in which different groups of diazotrophic organisms are added sequentially.
Diazotrophic growth formulations are similar to those of non-fixing phytoplankton, except that they are not limited by nitrogen 5 and have a higher N:P ratio ( ":$ & = 45). H1 introduces a generic autotrophic diazotroph; H2 replaces H1's generic diazotroph with two autotrophs representing unicellular and colonial (e.g., Trichodesmium spp.) cyanobacteria. The unicellular group overall follows the same formulation as the generic diazotrophs, except that no aggregation term is included. We simplify unicellular diazotroph behaviour in this manner because this group represents free-living picoplanktonic cells that typically do not form large colonies although they can aggregate (Bonnet et al. 2016). Instead, we assume this group is grazed by 10 zooplankton at similar rates as non-fixing phytoplankton. This difference between colonial and unicellular groups is consistent with studies suggesting that colonies represent an evolutionary adaptation to decrease grazing pressure (Nielsen, 2006). Aside from their size, Trichodesmium spp. colonies may be less palatable and harder to digest due to toxins (Kerbrat et al. 2010(Kerbrat et al. , 2011. Grazing is not a major fate of this group (O'Neil and Roman, 1994).
The last model version, H3, adds a heterotrophic group to the model structure of H2. This functional group is not limited by 15 light and grows by consuming both dissolved inorganic and organic forms of nutrients. An intermediate version H3' is used as a control, where the heterotrophic organisms do not fix nitrogen and are limited by the availability of nitrogen in inorganic forms and from small detritus. Model H3 eliminates the nitrogen limitation and the heterotrophic group becomes a heterotrophic diazotroph group. In a subsequent set of four experiments (H3a, H3b, H3c, and H3d), we remove complexity from H3. The heterotrophic group remains, but we sequentially remove the autotrophic groups one-at-a-time: first the colonial 20 cyanobacteria (H3a), then the unicellular cyanobacteria (H3b), then the generic autotrophic diazotroph (H3c), and finally we remove all autotrophic diazotrophs (H3d). A summary of all model versions is given in Table 1 and a description of state variables and full model equations in the Supplement.
We acknowledge that some of these model assumptions are still simplifications of diazotroph behaviour. For example, Trichodesmium spp., as well as other autotrophic and diazotrophic organisms, may also take up dissolved organic matter to 25 support growth (Benavides et al., 2017). Nevertheless, model assumptions and formulations used here are in line with the most commonly accepted understanding of the dominant controls on microbial and diazotrophic growth.

Parameter Optimization Method
Parameter optimization refers to the minimization of misfit between model and observations by adjusting model parameters. 30 We applied the method first to systematically calibrate the most sensitive parameters of H0 (see supplement), and then to independently re-calibrate parameters in H1 to H3 after the introduction of diazotrophs. We used an evolutionary algorithm, where changes in the parameter values follow a set of rules inspired by the process of natural selection (Houck et al., 1995;Kuhn et al., 2015). The algorithm starts with a randomly generated "population" of 30 parameter sets ( ⃗), which are iteratively modified over a number of generations. During each generation of the population, the cost ( ⃗) of the model with parameter set ⃗ is calculated as: 5 ∑ 4 6 7,9 − 7,9 ; where 6 represents a model estimate and the corresponding observation. N is the number of observations included for each variable . Here the number of variables V is 5 (nitrate + nitrite, ammonium, phosphate, chlorophyll-a, and oxygen measured as profiles at Station A between 2006 and 2010). Model-data misfits are weighted by the factor 7 = 1/ 7 , i.e, the inverse standard deviation of each variable. Half of the parameter sets with the lowest J value "survive" to the next generation. The 10 other half of the population is regenerated from new parameter sets obtained by recombination of two random "parent" sets drawn from the better performing half (i.e., the "survivors" of the previous generation). Parameters also "mutate", i.e. random noise is added, for additional variability in the parameter space. An allowable range of values is set for each parameter based on the literature (Table 2).

Optimized Parameters 15
The parameter optimization method has limitations. Most importantly, the optimization cannot estimate with confidence parameters that are unconstrained by the observations (Fennel et al. 2001;Schartau and Oschlies, 2003;Ward et al., 2010). To avoid this, a subset of H0's most sensitive parameters was selected for optimization through a preliminary sensitivity analysis.
Optimized parameters for H0 are identified in Table 2 along with the optimal values. The optimization was replicated 10 times over 100 generations using the algorithm described in section 3.3.1. Non-optimized parameters are fixed at their a priori 20 estimates based on Fennel et al. (2006Fennel et al. ( , 2013. For each model version with diazotrophs (H1, H2 and H3), some of the parameters already optimized for H0 required recalibration to properly accommodate the changes in system dynamics. Re-calibrated parameters for each model version are presented in Table 3. No re-calibration was performed for model versions (H0' and H3a-d), as they are aimed to test the relative importance of individual model components. 25

Diazotroph Parameters
Since none of the parameters directly related to the diazotroph groups are constrained by the available observations, they were predefined for H1, H2 and H3, based on the observational and modelling literature (Table 3). Previous modelling studies have used maximum growth rates of generic N2 fixers ranging from 0.4 d -1 (Moore, et al., 2004) to 1.25 d -1 (Ward et al., 2013).
cultured under various combinations of Fe and light availability exhibit maximum rates around 0.3 ± 0.05 d -1 (Capone et al., 1997;Berman-Frank et al., 2001;Hutchins et al., 2007). Growth rates can be higher (up to 0.5 d -1 ) at high CO2 and high light availability (Kranz et al., 2010;Hong et al., 2017). We chose a common reference maximum growth rate of 0.25 d -1 for all photosynthetic diazotrophs, such that differences between the model versions result only from the different assumptions about 5 the losses of each group (e.g., predation of unicellular cyanobacteria vs. sinking of large aggregates). Based on growth rates measured for cultured heterotrophic bacteria, we chose a value of 0.2 d -1 for the heterotrophic diazotrophs (Pomeroy and Wiebe, 2001). Observational and modelling studies were also considered to set the photosynthetic initial slope of photosynthetic diazotrophs (Geider et al., 1997;Moore et al., 2004;Hutchins et al., 2007). Other parameters are based on Fennel et al. (2002).

Sensitivity to Physical Nudging
Model runs, with and without temperature and salinity nudging towards observations, demonstrate that nudging has a 25 negligible effect below 200 m, indicating that horizontal advection does not modify the lower part of the water column in a significant way (Supplement). Above 200 m, nudging corrected model errors in the representation of vertical mixing and surface forcing. The average magnitude of the differences due to nudging in the top 200 m is 0.20±0.45 ºC and 0.5±0.16 kg m -3 . Since these effects are small and limited to the surface, we conclude that neither nudging nor the neglect of horizontal advection affects our conclusions significantly. Figure 3 shows simulated NO3 and PO4 concentrations from models H0', H0, H1, H2, H3' and H3, along with the corresponding measurements. Observed NO3 and PO4 concentrations exhibit a marked increase in deep water after the strong 5 winter mixing of 2008. Weaker winter mixing after 2008 results in deep-nutrient accumulation, which is more pronounced for NO3 than PO4. Model H0' reproduces some deep-nitrate accumulation, but underestimates NO3 concentrations in comparison to observations. Model H0 strongly underestimates inorganic nitrogen below the nutricline. Model H1, where N2 fixation was introduced via a generic autotroph, generates only small changes in the vertical distribution of nutrients. In model H2 the representation of NO3 below the nutricline is slightly improved; however, underestimation of mid-water NO3 is still noticeable.

Sensitivity to Planktonic Stoichiometry
To investigate whether the vertical distribution of dissolved inorganic nutrients was affected by our assumption of fixed N:P 25 phytoplankton and diazotroph ratios, we explored the sensitivity of NO3 and PO4 to changes in the non-fixing phytoplankton N:P ratio ( ":$ %& = 16) and the diazotrophs N:P ratio ( ":$ & = 45). In this analysis we used model version H1, which includes a single non-fixing group and a single N2 fixing group and varied the ratios one-at-a-time. The range of values for R D:E FG varied from 10 to 28 and the range for ":$ & from 19 to 51. Figure 5 shows examples of the results obtained by increasing and decreasing each ratio.
Changes in the N:P ratios had negligible effects on the vertical distribution of NO3, but strongly affected PO4 distribution. In general, lower than Redfield ":$ & of diazotrophs increases PO4 below 300 m, most strongly at depth. This occurs possibly because more phosphorus returns to the dissolved pool per unit nitrogen trough excretion and remineralization. In contrast, a 5 decrease of more than half the ":$ %& of non-diazotrophs produces only a minor decline in deep PO4, while increases in the ratio did not have a significant effect. The decline in deep PO4 occurs because, in the absence of nitrogen limitation, diazotrophs can utilize additional phosphorus.

Validation against Independent Observations
Observations from 2010 to 2014 (outside the optimization period) are used to independently validate the models. The rootmean-square errors (RMSEs) in Table 4 show that, in terms of chlorophyll, PO4 and surface O2, all models behave similarly and achieve similar agreement for assimilated and independent observations. As demonstrated in the previous sections, the model versions mainly diverge in their behaviour with respect to NO3, with some differences in O2 concentrations. Between 0 30 and 100 m, H3 has the largest RMSEs for NO3, but below 100 m it has the lowest values, particularly against unassimilated NO3. H3 also has the lowest RMSEs for surface and deep oxygen (Table 4). Compared against the other model versions, H3 increasingly overestimates surface NO3 over time. However, the deep NO3 inventory is best represented by H3. By the end of the observed time series, between 2013 and 2014, H3 starts to also 5 overestimate deep NO3.

Primary Production and N2 Fixation Rates
We now compare the simulated rates of primary production with those reported for the Gulf of Aqaba by Rahav et al. (2015) and Iluz et al. (2009) (Figure 8a) and the simulated rates of N2 fixation with those measured by Rahav et al. (2015) andFoster et al., (2009) (Figure 8b). Following Rahav et al. (2015), we show the rates at the DCM and their averages above and below 10 the DCM. The depth-resolved in situ primary production rates reported by Iluz et al. (2009) were also averaged in the same way for comparison. Where necessary, observations of primary production in carbon units were converted to the model's nitrogen units using the Redfield ratio.
Simulated primary production above the DCM ranges from 0.02 to 0.85 mmol N m -3 d -1 and exhibits an annual cycle with peaks of productivity in October and April. A prolonged period of low primary production extends from April to September 15 in most model versions. Model versions H3b and H3d maintain rates twice as large as the rest of the models during the summer/fall period. Aside from H3b and H3d, differences between models are small and simulated rates agree with those measured by Iluz et al. (2009) and Rahav et al. (2015).
Above the DCM, models H1, H2, H3 and H3a show a well-defined N2 fixation peak during summer months (i.e., after the peak in primary production). Maximum rates in these models range from 0.001 to 0.1 mmol N m -3 d -1 , which agrees with the 20 observed rates by Foster et al. (2009) and Rahav et al. (2015). In general, simulated N2 fixation rates are low during winter and spring. Similar temporal patterns and differences between model versions occur at the DCM and below. Peaks in N2 fixation at these depth levels occur after the surface peak, and have a shorter duration and smaller amplitude. Deep N2 fixation rates estimated by models without heterotrophic diazotrophs do not match the observed rates by Rahav et al. (2015).

Is N2 Fixation Relevant in the Gulf of Aqaba?
In this study we implemented and optimized a series of models with different assumptions about N2 fixation in the Gulf of Aqaba. The models range from one neglecting N2 fixation to another assuming that, in addition to two autotrophic diazotroph groups, N2 fixation can occur in the entire water column (i.e., independent of light availability). While the models are very similar in their abilities to replicate chlorophyll and PO4, model H3 performed the best in reproducing the observed pattern of 30 deep-NO3 accumulation and O2. Overall, all models that consider N2 fixation accumulate nitrogen at different rates, as they enrich the nitrogen content of detritus, which is then remineralized at depth over time.
The best model performance was obtained with two groups of autotrophic organisms and a group of heterotrophic organisms (H3). A model without explicit N2 fixation, but in the absence of sediment denitrification, also increases the accumulation of deep NO3 in a similar fashion as version H2 but not sufficiently high to match the observations. This suggests that N2 fixation in the area must exceed denitrification rates. In the models with denitrification, the average sediment denitrification flux is 5 0.25 ± 0.46 mmol N m -2 d -1 , with a maximum value of 3.01 mmol N m -2 d -1 . These values at are the lower end of a global compilation of sediment denitrification rates by Fennel et al. (2009), which have a mean of 2.2 mmol N m -2 d -1 and maximum values exceeding 10 mmol N m -2 d -1 .
Observations from the Gulf of Aqaba exhibit an excess of nitrogen that contrasts to exterior waters from the Arabian Sea and Indian Ocean, which are considered net nitrogen sink regions (Gruber and Sarmiento, 1997). There are too few reported values 10 of dissolved inorganic nitrogen-to-phosphorus ratios for the Red Sea region from Bab-el-Mandeb to the Strait of Tiran to provide a complete idea of its spatial distribution; however, the limited available information supports our conclusions. Naqvi et al. (1986) found a significant difference in N:P ratios between surface incoming and sub-surface outflowing waters at Babel-Mandeb, concluding that N2 fixation was a process required to explain these anomalies in the nitrogen budget. Higher nitrogen concentrations have been observed in the Red Sea in comparison to the Arabian Sea and Indian Ocean, where a strong 15 deficit of nitrogen develops as losses due to denitrification exceed the input of newly fixed nitrogen (Burkill et al., 1993;Naqvi, 1994;Gruber and Sarmiento, 1997;Morrison et al., 1998Morrison et al., , 1999. Close to the entrance of the Persian Gulf, average excess phosphate has been estimated to be around 5 mmol m -3 at all depths and seasons reported, with maximum excess values on the order of 8 mmol m -3 (Gruber and Sarmiento, 1997). Thus, it has been hypothesized that limited deep-water exchange at Babel-Mandeb allows waters of the Red Sea outside of the Gulf of Aqaba to acquire different characteristics from inflowing 20 Arabian Sea waters (Naqvi et al., 1986). Our model results support this hypothesis and suggest that N2 fixation is key for the formation of the distinct biochemical characteristics in the Gulf of Aqaba. Considering the regional context, our models suggest that, despite low rates, N2 fixation is necessary to explain the nitrogen vertical distribution in the Gulf of Aqaba, and the interannual accumulation of deep nitrate during years with weak convection.

How does N2 fixation Contribute to Primary Production?
In this section we discuss the contribution of N2 fixation to primary production in the Gulf of Aqaba, and our quantitative estimates of N2 fixation with respect to global rates (Figures 8-9). Our estimates of surface primary productivity agree with those reported by Iluz et al. (2009) for March-April of 2008. However, our models overestimate surface primary productivity values in 2010 when compared to those reported by Rahav et al. (2015). On average, our best-performing model version yields 30 an annual primary production rate of 304±56.9 g C m -2 yr -1 (H3). This rates is higher than previously published annual averages, which range from 80 g C m -2 y -1 (Levanon-Spanier et al., 1979) to 170 g C m -2 y -1 (Iluz, 1991), whereas more recent unpublished primary production estimates at IUI range between 141 and 197 g C m -2 y -1 (pers. comm. Y. Shaked).
The ratio of new to total primary production (f-ratio) in our experiments ranges from 15% to 80%. Maximum f-ratios are estimated in January and February due to significant contributions from deep NO3, whereas f-ratios are at their minimum during stratified conditions (June -August). Our best-performing model version, H3, estimates a summer minimum f-ratio 5 0.22. The average f-ratio for all scenarios is 0.47. This agrees with published estimates for the Gulf of Aqaba of 0.5 during the stratified period as determined from a nitrate-diffusion model (Badran et al., 2005).
Total annual N2 fixation rates from our best-performing model versions (H3 and H3a) are similar to high estimates reported for other regions (Capone and Carpenter, 1982;Michaels et al., 1996;Lee et al., 2002), while those obtained in the other model experiments are within the range of values reported for the Gulf of Aqaba. The intensity of winter mixing has a minor effect 10 on N2 fixation rates; the largest effect occurred in H3a where N2 fixation increased by 15% after deep winter mixing. Based on our best-performing model version (H3), we estimate that 10% to 14% of the total primary production is supported by N2 fixation.

Are Heterotrophic N2 Fixers Important?
In contrast to previous models (e.g., Hood et al., 2001;Fennel et al., 2002;Monteiro et al., 2010;Moore et al., 2004), our 15 model version H3 relaxes the assumption of light dependence for diazotrophy, through the inclusion of heterotrophic diazotrophs in addition to two groups of autotrophic diazotrophs. This model improves the representation of NO3 and O2 at depth (Figures 3, 6). Changes in deep NO3 can be explained through the enrichment of detritus, while changes in O2 may reflect the additional sink of O2 at depth due to the heterotrophic group. All model versions with heterotrophic organisms also match observed estimates of N2 fixation in deep waters of the Gulf of Aqaba, and without them N2 fixation rates below the 20 DCM are underestimated (Figure 8). Heterotrophic N2 fixation also impacts total N2 fixation ( Figure 9).
For instance, N2 fixation rates in mesopelagic and abyssopelagic waters down to 2000 m (Fernandez et al., 2011;Bonnet et al., 2013;Loescher et al., 2014) have been attributed to non-cyanobacterial organisms, including proteobacteria (Turk-Kubo et al., 2014). These nifH-expressing heterotrophic phylotypes can be as abundant as unicellular cyanobacterial groups and dominate 25 the deep and dark zones of the water column (Church et al. 2005;Langlois et al., 2005;Riemann et al., 2010). Genetic evidence and rate estimates from the Gulf of Aqaba suggest that nifH expressing heterotrophic proteobacteria α and γ may explain the correlation of bacterial productivity rates with N2 fixation rates (Rahav et al., 2013;. Aside from the Gulf of Aqaba, aphotic N2 fixation and nifH gene expression have also been reported in the Baltic Sea (Farnelid et al., 2013), Arabian Sea (Jayakumar et al., 2012) and Mediterranean Sea (Rahav et al., 2013). 30 Despite the importance of heterotrophic diazotrophs in our model, the simulated colonial diazotroph blooms are responsible for the highest N2 fixation rates, so they are a necessary model aspect to achieve resemblance with the observations. This is in line with evidence of extensive blooms of Trichodesmium spp. being responsible for the high N2 fixation rates observed in the Arabian Sea and Red Sea (Capone et al., 1998;Post et al., 2002;Foster et al., 2009). In the northern Gulf of Aqaba, colonies and free trichomes of Trichodesmium spp. are found throughout the year down to 100 m depth (Post et al., 2002). Ephemeral blooms of T. erythraeum and T. thiebautii have been documented near the coast of Eilat (Post et al., 2002, Gordon et al., 1994Kimor and Golandsky, 1977). However, massive blooms are rare in the Gulf of Aqaba (Foster et al., 2009;Mackey et al., 5 2007, pers. Communication Berman-Frank) and the model probably overestimates the contribution of Trichodesmium spp.'s annual blooming to total N2 fixation rates, as seen in the much larger surface N2 fixation rates generated by H2 and H3. As new observational information is collected, further model refinements may be necessary to better reflect the actual contribution of different diazotrophic groups in the Gulf of Aqaba.

Limitations and Uncertainties 10
The one-dimensional nature of our physical setting, which neglects the contribution of horizontal advection to the vertical structure of simulated tracers, can be considered a limitation of this study. This simplification is, however, necessary to perform model calibration and test multiple model structures at a manageable computational expense. We applied temperature and salinity nudging to ensure accurate representation of the vertical density structure. Comparison of the simulated vertical structure with and without nudging shows that this correction has negligible effects on deep waters, where the effects of N2 15 fixation are the most relevant. This is consistent with the existing literature about circulation of the Gulf of Aqaba, which describes how geomorphology and bathymetry limit water exchange between the Gulf of Aqaba and the Red Sea to the upper 300 m (Wolf-Vecht et al., 1992;Biton and Gildor, 2011). It is, therefore, unlikely that horizontal transport could explain the observed accumulation of deep NO3. Nevertheless, transport of nitrogen-enriched sub-surface waters from the Gulf of Aqaba towards the exterior may dampen the long-term accumulation of nitrogen (Figure 7). 20 There are other sources of nitrogen that were not explored in the present study that we discuss here briefly. For instance, we did not include contributions to N2 fixation by diatom-diazotroph associations, which are significant in other regions. While diatom-diazotroph associations have been detected in the Gulf of Aqaba, they are not as abundant as unicellular diazotrophs, Trichodesmium and proteobacteria (Kimor et al. 1992, Foster et al., 2009. In general, due to the oligotrophic characteristics of the region, small phytoplankton species (<8 µm) contribute more than 90% of the chlorophyll-25 a standing stock (Lindell & Post 1995, Yahel et al. 1998. Dinoflagellates and diatoms together correspond to less than 5% of the phytoplankton biomass, except during ephemeral diatom blooms during spring when they can account for nearly 50% of the total biomass (Al-Najjar et al., 2006).
Another source of nitrogen which could significantly affect this region is atmospheric deposition, as the Gulf receives considerable dust input from the surrounding deserts. Recently, it has been shown that atmospheric dust input does not correlate 30 with chlorophyll variability in surface waters of the Gulf of Aqaba (Torfstein and Kienast, 2018). A previous study suggested that atmospheric deposition of nitrogen could support over 10% of surface primary production in the region, based on measurements of local aerosol composition and a dust deposition model (Chen et al., 2007). However, this estimate had a relatively large uncertainty due to errors associated with the deposition flux calculation and the temporal variability in dust flux (Chen et al., 2007). Moreover, very low nitrogen concentrations and N:P ratios lower than Redfield from the surface down to 80 m were observed during the same time period (Foster et al., 2009). Therefore, the role of atmospheric nitrogen inputs remains uncertain. 5 In contrast to the improvement in NO3 distributions with the addition of heterotrophic diazotrophs, all other model versions exhibit similar underestimation of deep total PO4. This suggests that their structure lacks a process affecting this nutrient.
Several processes can independently affect PO4. Prior knowledge suggests that seasonal advection of nitrogen depleted surface waters may contribute to changes in PO4, while having minimum impact on NO3 concentrations. Other relevant missing processes include variable plankton stoichiometry and PO4 remineralization. Our sensitivity analysis highlighted that varying 10 the planktonic N:P ratios has an evident effect on PO4 but does not affect NO3 distributions significantly. A related effect was noticed when comparing fixed constant vs. variable plankton nutrient uptake stoichiometry in a model of the central Baltic Sea: seasonal particulate organic C:N ratios were very similar, but C:P diverged (Kreuss et al., 2015). Moreover, preferential PO4 remineralization may more directly increase deep-PO4 concentrations as organic matter is decomposed in the bottom layers. The inclusion of this process in a model of the North Atlantic Ocean improved the representations of biogeochemical 15 characteristics of the area and increased the N2 fixation rates obtained by the model . Variable nitrogen allocation within diazotrophs also improved the performance of a model representing annual cycles in the North Atlantic sub-tropical gyre (Fernández-Castro et al., 2016). Further investigation of this process is recommended.
Finally, an intrinsic caveat that should accompany all biological models is the uncertainty associated with parameter values (Denman, 2003). We reduced this uncertainty by using systematic parameter optimization for a subset of parameters, yet others 20 were kept fixed at a priori assumed values. In the case of N:P ratios, a simple sensitivity analysis suggests that they do not affect conclusions with respect to inorganic dissolved nitrogen but modify phosphorus concentrations. In the case of parameters related to diazotrophs, these are largely unconstrained by the observations we used, which are typically available in long-term data sets (i.e., Chl-a, NO3, PO4, O2). We emphasize that results in this study are experimental, testing the effects of assumptions about diazotroph behaviour rather than modifying such behaviour by subjectively tuning parameter values. 25

Conclusions
We implemented and optimized biogeochemical models that represent a range of different assumptions about diazotrophy in a 700 m-deep pelagic station from the northern Gulf of Aqaba. Our model results demonstrate the importance of N2 fixation in replicating the observed water-column-integrated nitrogen and oxygen inventories. The model without N2 fixation is unable to replicate the observed vertical structure of inorganic nitrogen. The models that include diazotrophs significantly modify this 30 variable by increasing the fraction of remineralized nitrogen from organic matter decomposition. The effect of N2 fixation on O2 distributions depends on the type of nitrogen fixer. N2 fixation by autotrophic-photosynthetic organisms increases oxygen concentrations, while heterotrophic organisms decrease deep-oxygen due to increased respiration and organic matter decomposition. The observed vertical structure of NO3 and oxygen is reproduced best with a model that includes heterotrophic, and colonial and unicellular autotropic diazotrophs suggesting that heterotrophic N2 fixation is necessary to explain the observed excess nitrogen at this location. N2 fixation assumptions do not affect PO4 concentrations significantly, but they are affected by assumptions about N:P ratios of organic matter. The N2 fixation rates simulated by this model are similar to the 5 highest observational estimates from the Gulf of Aqaba. Aphotic N2 fixation is simulated to occur at lower rates than maximum autotrophic N2 fixation yet occurs continuously over a large portion of the water column. This suggests that heterotrophic diazotrophs set a background rate of N2 fixation in the ocean that should be considered further in global estimates and biogeochemical models.         Table 1for key to different simulations). A summary of previous estimates of N2 fixation rates in observational and model studies is included in (b). Redfield C:N ratio is used to transform the model results nitrogen units to carbon units. Control versions H0' and H3' are not shown.   Fahnenstiel et al. (1995) c. Veldhuis et al. (2005) d. Fennel et al. (2006) e. (Lima and Doney (2004) f. Ward et al. (2013) g. Moore, et al.(2002) h. Geider et al. (1997) i. Gifford et al. (1995) j. Landry et al. (1984) k. Tande and Slagstad (1985) l. Amon and Benner (1996) m. Enríquez et al. (1993) n. Smayda and Bienfang (1983)