Simulated Future Trends in Marine Nitrogen Fixation Are Sensitive to Model Iron Implementation

Biological nitrogen fixation is an important oceanic nitrogen source, potentially stabilizing marine fertility in an increasingly stratified and nutrient‐depleted ocean. Iron limitation of low latitude primary producers has been previously demonstrated to affect simulated regional ecosystem responses to climate warming or nitrogen cycle perturbation. Here we use three biogeochemical models that vary in their representation of the iron cycle to estimate change in the marine nitrogen cycle under a high CO2 emissions future scenario (RCP8.5). The first model neglects explicit iron effects on biology (NoFe), the second utilizes prescribed, seasonally cyclic iron concentrations and associated limitation factors (FeMask), and the third contains a fully dynamic iron cycle (FeDyn). Models were calibrated using observed fields to produce near‐equivalent nutrient and oxygen fits, with productivity ranging from 49 to 75 Pg C yr−1. Global marine nitrogen fixation increases by 71.1% with respect to the preindustrial value by the year 2100 in NoFe, while it remains stable (0.7% decrease in FeMask and 0.3% increase in FeDyn) in explicit iron models. The mitigation of global nitrogen fixation trend in the models that include a representation of iron originates in the Eastern boundary upwelling zones, where the bottom‐up control of iron limitation reduces export production with warming, which shrinks the oxygen deficient volume, and reduces denitrification. Warming‐induced trends in the oxygen deficient volume in the upwelling zones have a cascading effect on the global nitrogen cycle, just as they have previously been shown to affect tropical net primary production.

Improving the constraints on future N 2 fixation trends requires an improved understanding of the processes critical for modern ocean N 2 fixation. Spatial patterns of N 2 fixation are thought to be determined by temperature (Breitbarth et al., 2007;F. Fu et al., 2014), availability of micro and macro nutrients, and grazing pressure. All three drivers are susceptible to warming effects from climate change (Laufkötter et al., 2015). Due to their effects on the macro-nutrient ratio, denitrification and annamox are considered key processes modulating the activity of diazotrophs. Denitrification and anammox that occur in low oxygen environments (such as in tropical upwelling zones) change the relative availability of nitrate (NO 3 ) and phosphate (PO 4 ) in the surface ocean downstream (Deutsch et al., 2007;Weber & Deutsch, 2014). The resulting stoichiometrically excess phosphate (positive P*, defined as PO 4 -NO 3 /16, Deutsch et al., 2007) might provide diazotrophs with a competitive advantage relative to other types of phytoplankton. This competitive advantage is hypothesized to explain why rates of N 2 fixation are enhanced in the NO 3 -depleted low latitude surface ocean (Deutsch et al., 2007;N. S. Garcia et al., 2015;Tyrrell, 1999).
At the same time, iron (Fe) is essential for all phytoplankton as it is needed to form a range of cellular enzymes that are associated with photosynthesis and respiration. However, diazotrophs additionally require iron for producing the N 2 fixation enzyme nitrogenase (Kustka et al., 2003;Raven, 1988;Tagliabue et al., 2017;Zehr & Capone, 2020), causing diazotroph growth to be particularly sensitive to a lack of iron. For diazotrophs, the supply of iron is so essential that its availability helps to define the biogeography of N 2 fixation Knapp et al., 2016; C. M. Moore et al., 2009;Schlosser et al., 2014;Ward et al., 2013). Unfortunately, global ocean models still have a hard time reproducing the observed iron concentrations (Tagliabue et al., 2016) and simulated patterns of N 2 fixation also have large uncertainties (Landolfi et al., 2018).
Iron availability is expected to affect denitrification patters in the ocean (Landolfi et al., 2013). Biological productivity is limited by iron availability over one-third of the surface ocean, including the Eastern Tropical Pacific (ETP) upwelling zones (C. M. Moore et al., 2013). This limitation can impact particulate organic matter export (Martin, 1990). Therefore, a change in iron availability in a highly productive upwelling zone might alter particulate organic matter export and thereby influence the extent and intensity of oxygen deficient zones (ODZs) and of denitrification. In their modeling study, Buchanan et al. (2019) showed that an increase in iron availability in the ETP can lead to a larger ODZ, and a downstream enhancement of diazotroph biomass and N 2 fixation. They inferred that this process could have contributed an additional 7-16 ppm of carbon dioxide (CO 2 ) sequestration in the tropical ocean in the Last Glacial Maximum (Buchanan et al., 2019).
In addition to bottom-up controls via iron availability and denitrification, top-down control of diazotrophy by zooplankton grazing might also influence the global patterns of N 2 fixation (Landolfi et al., 2021;Wang et al., 2019). Landolfi et al. (2021) present a top-down control framework which explains the maintenance of a niche for diazotrophs in nutrient-rich environments by selective grazing pressure on non nitrogen-fixing phytoplankton.
Here we examine what controls the global N 2 fixation response to global warming in three marine biogeochemical models with different representations of the marine iron cycle, all of which perform similarly with respect to reproducing observed NO 3 , PO 4 , and O 2 distribution (Yao et al., 2019). All three models were individually subjected to parameter optimization and their preindustrial steady-states were assessed in Yao et al. (2019). Our objective is not to simply assess the role of iron limitation, as has been done previously, for example, J. K. Moore and Doney (2007), but to examine how differences in nutrient cycling impact transient model behavior. We first describe the models and the new transient global-warming experiments. Then, we compare the different simulated responses of N 2 fixation to global warming. Next, we examine the impact of ODZ development and denitrification on N 2 fixation. We link the development of ODZs to in situ export production trends via the different assumptions about nutrient pathways in the different models, and make some recommendations for further model improvements.

Model Description
In this study, we use the University of Victoria Earth System Climate Model (UVicESCM, version 2.9; Eby et al., 2013;Weaver et al., 2001). The UVicESCM is an earth system model of intermediate complexity. It contains four components: a simple one-layer atmospheric model, a terrestrial model (Meissner et al., 2003), a sea-ice model and a three-dimensional ocean general circulation model. All components have a horizontal grid resolution of 3.6° in meridional and 1.8° in zonal direction. The one-layer atmospheric model is based on an energy-moisture balance scheme (assuming lateral heat and moisture transport by diffusion; Fanning & Weaver, 1996), which calculates heat and water fluxes between the atmosphere and the ocean, land and sea ice dynamically, while applying prescribed (NCAR/NCEP) winds. A thermodynamic sea-ice model (Bitz & Lipscomb, 1999) is used, while constant continental ice sheets are prescribed. The ocean component is the Modular Ocean Model 2 (MOM2). It has 19 vertical layers with increasing thickness from 50 m at the surface to 500 m in the deep ocean. The ocean advection scheme is a non-linear second-order Flux-Corrected Transport scheme (Weaver & Eby, 1997).

Biogeochemical Models and Parameterization
We apply three variants of the marine biogeochemical component in the UVicESCM, which differ in their treatment of the iron cycle. All of them are nutrient-phytoplankton-zooplankton-detritus (NPZD) models with two types of phytoplankton, ordinary phytoplankton and diazotrophs, one zooplankton, and detritus. In the models, the growth of phytoplankton is governed by nutrients, light and temperature. In all three models, diazotrophs grow slower than ordinary phytoplankton but can also fix N 2 to meet their N demand when NO 3 is limiting (equations are given in appendix).
The first model contains no representation of the iron cycle (NoFe, created by disabling iron limitation in the model presented by Keller et al. [2012]). The second model variant is that of Keller et al. (2012) and makes use of a prescribed, seasonally periodic iron concentration mask (FeMask; Keller et al., 2012), by regridding the output of the Biology Light Iron Nutrient and Gases model (BLING; Galbraith et al., 2010). The third variant includes a fully dynamic iron cycle (FeDyn; Nickelsen et al., 2015), which resolves major sources and sinks of iron in the ocean, for example, dust deposition, benthic release, particle scavenging, colloid formation, remineralization as well as hydrothermal release (Yao et al., 2019). In the models that include iron (FeMask and FeDyn), both diazotroph and ordinary phytoplankton growth rates are reduced when iron is limiting, in addition to the usual controlling factors, for example, macro-nutrient(s), light and temperature. In the model without iron (NoFe), the growth rates of diazotrophs and ordinary phytoplankton do not depend on the iron concentration. The use of an externally derived iron mask (BLING) in FeMask, rather than one produced by FeDyn was decided for pragmatic reasons. FeMask has been, and continues to be, widely used in the UVicESCM user community for a variety of research questions (Getzlaff & Dietze, 2013;Mengis et al., 2018;Somes & Oschlies, 2015). Therefore it is useful for us to calibrate it, and to assess how its optimal configuration behaves transiently compared to other model versions. However, application of a FeDyn-derived iron mask to the FeMask calibration would also prove redundant, as the model structures are identical with respect to iron uptake by biology. BLING (used in FeMask) and FeDyn has different iron concentrations in the surface due to a lot of different factors, for example, different dust deposition, sediment release, biological uptake, scavenging parameterization, and different physics. This is why iron half saturation parameters for phytoplankton were included in the model calibration by Yao et al. (2019). Details, including equations, of each model configuration can be found in Nickelsen et al. (2015, FeDyn) and Keller et al. (2012, FeMask and NoFe). Note that the objective of our current study is not a comparison of full, partial or absent consideration of the influence of iron limitation in a preindustrial steady state, which has been dealt with by Yao et al. (2019), but instead a comparison of model behavior under transient forcing. These models have different parameterizations resulting from parameter calibration (see details below) against the same set of observational nutrient and oxygen distributions. These model differences have implications for nitrogen fixation responses to climate warming.

Model Calibration
The models used in this study have different structures and were hand-tuned to different objectives in prior studies (Keller et al., 2012;Nickelsen et al., 2015). Hence, they have different representations of ocean biogeochemical tracer distributions. In order to perform a "fair" comparison between different models, we first needed to calibrate all models using the same objective against the same set of observational data, which was done in an earlier study by Yao et al. (2019).
In Yao et al. (2019), we first chose a subset of biogeochemical parameters for FeDyn (six out of 44; see Table B1) for its calibration. Those parameters are among the most sensitive and/or most uncertain (few or no measurements available) and at the same time provide a large coverage of control in the different biogeochemical domains, for example, iron, nitrogen, oxygen, phytoplankton and zooplankton biomass. For the calibration of NoFe and FeMask, we chose sets of parameters that are as close to the FeDyn set as possible (see Table B1). We then calibrated all three models against the same global set of observational NO 3 , PO 4 , and O 2 concentrations (World Ocean Atlas 2013; H. E. Garcia et al., 2013aGarcia et al., , 2013b with the help of an evolutionary optimization algorithm, which minimizes the model misfits. After calibration, we found that all models have relatively similar performance reproducing the observed NO 3 , PO 4 , and O 2 distributions (less than 5% difference between rms model-data misfits) in an assumed preindustrial steady-state. Residual model misfits for NO 3 , PO 4 , and O 2 are only about 15% of their observed global mean concentration.
We next assessed each model's synthetic global biogeochemical fluxes and indicators, for example, NPP, export production, flux of organic carbon at a depth of 2 km, denitrification, the volume fraction of the ocean with suboxic conditions (oxygen concentrations less than 5 mmol O 2 m −3 ), and the volume fraction of the ocean with hypoxic conditions (oxygen concentrations less than 50 mmol O 2 m −3 ). We found that in steady-state, all models simulate observed global biogeochemical fluxes and indicators (four out of six) to a similar degree. Although NoFe simulates a 50% higher NPP compared to the model with the lowest estimate (FeMask), the range is 49.1-74.8 Pg C yr −1 between the three calibrated models. This spread is smaller than the 32-77 Pg C yr −1 range reported from 24 independent estimates (based either on ocean color observations, on ecosystem models, or on biogeochemical models -a wide range of methodologies) in Carr et al. (2006). Different flux values arise between calibrated models as the model calibrations seek to achieve optimal fit to nutrient distributions using the parameters available in the different models, whose structures vary (presence or absence of iron limitation, namely). For example, in the models which include iron limitation, the calibration exploited this bottom-up control to limit NPP. However, in NoFe the calibration compensated for the absence of iron limitation using top-down control, with a larger resulting grazing parameter value and also greater global NPP rates. It is interesting to note that the wide range of NPP rates achievable with similar model performance with respect to nutrient and oxygen distributions indicates only minor relevance of NPP for general model tuning of these metrics. This is because NPP can be compensated by opposing fluxes in the model. In Yao et al. (2019) the surface total nutrient recycling (sum of the microbial loop recycling, zooplankton excretion and remineralization) compensates the differences of NPP and produces similar export production rates (6.3-7.1 Pg C yr −1 ) out of the euphotic zone (130 m) for all three calibrated models. It also highlights the importance of including additional metrics, as a target for model calibration when attempting to produce a complete picture of ocean biogeochemistry. Another indicator exhibiting variability between the three calibrated models is the volume fraction of the ocean with suboxic conditions. This is not surprising considering the fraction of suboxic water occupies only around 0.2 (World Ocean Atlas 2013; H. E. Garcia et al., 2013b) of the total ocean volume, and our calibrations were focused on the overall NO 3 , PO 4 , and O 2 distributions. More details about our calibration and resulting differences in optimal nutrient pathways for each model can be found in Yao et al. (2019).

Model Experiments
All models are integrated for 10 thousand years with preindustrial (year 1800) climate boundary conditions to reach equilibrium. A prescribed historical and "business-as-usual" RCP8.5 atmospheric CO 2 concentration forcing (van Vuuren et al., 2011) is then applied from year 1800-2100, but with climatological seasonally cycling NCAR/NCEP wind fields. The preindustrial physical circulation and circulation responses are identical in all three model configurations. Atmospheric iron deposition flux (Luo et al., 2008), which is a preindustrial climatological estimate of monthly iron deposition, is not changed during the transient experiment in FeDyn. We do not include atmospheric N deposition in this study, which was found to be of minor importance in an earlier study by Landolfi et al. (2017).

Simulated Trends in Global NPP and Diazotroph Primary Production
Global total primary production in year 1800 is 74.8 (NoFe), 49.4 (FeMask), and 50.4 (FeDyn) Pg C yr −1 . In all simulations, NPP decreases with time ( Figure 1a), in general agreement with the trends of CMIP5 models Laufkötter et al., 2015). Laufkötter et al. (2015) report a range of 15% to +30% (4.3 to +10 Pg C yr 1 ) NPP change by year 2100 from baseline values of 17-83 Pg C yr −1 for a 1998-2007 climatology. In the year 2100, the largest absolute decline is found in NoFe (5.6 Pg C yr −1 ) and the largest percentage decline in FeMask (9.16%). Differences in projected NPP are not due simply to the presence or absence of iron in the models; each model had previously been calibrated against observed biogeochemical tracer distributions in an objective model-tuning exercise to achieve as-similar-as possible model performance against PO 4 , NO 3 , and O 2 . The calibration framework achieved this objective by adjusting biogeochemical parameter values in a way that also produced differences in nutrient pathways and differences in fluxes and biogeochemical rates (Yao et al., 2019). However, despite these remaining differences in fluxes between models, the inter-model NPP spread in our calibrated model suite with identical physics is smaller than that found across CMIP5 models (Yao et al., 2019).

Trends in Simulated N 2 Fixation
Steady-state rates as well as trends in N 2 -fixation also show remarkable differences across the model suite. The simulated preindustrial N 2 fixation rates are between 30.3 and 54.1 Tg N yr −1 (Figure 2a). The calibrated NoFe has the smallest preindustrial N 2 fixation rate, which is also much smaller than earlier, hand-tuned versions of the UVicESCM that applied different parameter values and also did not include iron (UVicESCM version 2.7, 316 Tg N yr −1 , Riche & Christian, 2018; UVicESCM version 2.8 & 2.9, 150 Tg N yr −1 , Keller et al., 2012;Oschlies et al., 2019; UVicESCM version 2.9 with additional structural changes including benthic denitrification, 110-126 Tg N yr −1 , Landolfi et al., 2017). This low N 2 fixation value occurs in NoFe because of the relatively efficient near-surface nutrient recycling in this model, which minimizes denitrification and therefore N 2 fixation in the preindustrial steady-state configuration (Yao et al., 2019). Applying a dynamic iron cycle (FeDyn) in the model results, after calibration, in an optimal parameter set which produces the largest preindustrial N 2 fixation Figure 2. Simulated running 5-year-averaged global N 2 fixation rates (a), their temporal changes (b), temporal changes of denitrification (c) and the changes in ODZ volume (d, change in the percentage of the total ocean volume) between 1800 and 2100. Note: we define ODZ as regions where O 2 concentration is below 5 mmol m −3 . A global map of N 2 fixation is shown in Figure B4. rate (54.1 Tg N yr −1 ) among the three models. However, even the FeDyn N 2 fixation rate is in the lower range of observations and model estimates (between 89 and 213 Tg N yr −1 ; Landolfi et al., 2018). Note, that independent estimates of the N 2 fixation rate were not used as an explicit constraint in the calibration of Yao et al. (2019). Furthermore, none of the models include benthic denitrification, which could increase global total denitrification and enhance global N 2 fixation (Weber & Deutsch, 2014).
The simulated evolution of N 2 fixation in response to global warming differs considerably between simulations (Figures 2a and 2b). While N 2 fixation increases by 21.6 Tg N yr −1 (71%) in NoFe between 1800 and 2100, fixation rates remain almost unaffected by global warming in FeMask and FeDyn. Landolfi et al. (2017) reported an increase by 0.3% in N 2 fixation due to the increasing PO 4 limitation for the same simulated period in a version of UVicESCM with prescribed iron concentrations (which includes a representation of nitrous oxide and has a different set of parameter values comparing to FeMask). Another similar result is found in a model comparison study by Wrightson and Tagliabue (2020), where the only model that had a strong increase in N 2 fixation (58 Tg N yr −1 ) between years 1800 and 2100 was a model with neither iron nor phosphate limitation on (diagnosed) N 2 fixation. However, two models with neither iron nor explicit diazotrophs in their study showed large decreases in N 2 fixation (38.8-50.1 Tg N yr −1 ) over this time period, due to the tight spatial coupling in the upwelling zones between increasing phosphate limitation, reducing denitrification rates, and the subsequently diagnosed nitrogen fixation (where nitrate is restored to the Redfield ratio). As in our study, models compared in Wrightson and Tagliabue (2020), which include a representation of iron limitation project no substantial change in marine nitrogen fixation (a slight increase by up to 0.4 Tg N yr −1 ) to decreasing N 2 fixation (up to 23.2 Tg N yr −1 ). For iron-limited models, they attribute discrepancies across simulated global N 2 fixation trends to differences in the balance between decreasing fixation rates in the Atlantic and Indian Ocean sectors (for which their model ensemble shows a high degree of agreement) and fixation trends in the tropical Pacific sector (for which their model ensemble disagrees on both magnitude and sign). J. K. Moore and Doney (2007) demonstrated that iron stress in diazotrophs have a mitigating effect on global nitrogen cycle responses to perturbation by limiting the ability of diazotrophs to fix nitrogen when upstream denitrification is enhanced, with stabilizing feedbacks between N 2 fixation and denitrification operating least efficiently in the tropical Pacific.

Simulated Eastern Tropical Pacific Ecosystem Response to Warming
We next look at differences in the model responses to global warming in the ETP, due to its large impact on global N 2 fixation trends J. K. Moore & Doney, 2007;Wrightson & Tagliabue, 2020 and discussed below). This region has also been previously described as displaying large inter-model differences in nitrogen fixation response to warming (Wrightson & Tagliabue, 2020). Furthermore, Tagliabue et al. (2020) demonstrate biological iron limitation within this region to be of critical importance in determining productivity trends across the equatorial Pacific, which we also suspect to apply to nitrogen fixation trends, based on the global differences between iron-limited models and NoFe, shown above. However, it is also important to note that the different model behavior we report for the ETP also applies to the eastern tropical Atlantic upwelling region. Thus differences in global model trends are not only resulting from differences within the ETP.
As the temperature increases, NPP, remineralization, microbial loop recycling, and other temperature-sensitive fluxes accelerate accordingly (Figure 3). In all models, the surface ocean nutrient recycling pathways through enhanced and shallower detrital remineralization and the enhanced microbial loop gain importance under warming. Also, the increase in phytoplankton potential growth rate is greater than the increase in zooplankton potential grazing rate (which is capped above 20°C; Keller et al. (2012), in all models). Hence, total plankton biomass (ordinary phytoplankton, diazotrophs and zooplankton) increases in all models in the ETP.
In the NoFe model, light limitation due to phytoplankton self-shading is the dominant control on phytoplankton growth in the ETP ( Figure B5). Grazing control on the biomass is also stronger compared to the other models (NoFe has the highest optimized maximum grazing rate at 0°C, see Table B1). With ocean warming, the increase in phytoplankton potential growth rate is larger than the increase in zooplankton potential grazing rate in the ETP. This creates a 23% higher total plankton biomass in the region, which leads to a 35% higher net detritus production (sloppy feeding of zooplankton plus plankton mortality and lysis, minus zooplankton grazing on detritus, Figure 3). Warming also enhances detrital remineralization. However, a 41% increase in detrital remineralization 10.1029/2020GB006851 8 of 25 ( Figure 3) is not enough to keep up with the 401.0 Tg C yr −1 increase in net detritus production. Therefore, export production in the ETP increases by 15% (Figure 3), contributing to an increase in the simulated ODZ volume.
Phytoplankton are limited in the ETP by Fe and NO 3 ( Figure B5) in both models with iron, which makes their production more sensitive to decreasing nutrient supply under warming-induced intensification of stratification. In FeDyn and FeMask, euphotic zone detrital remineralization accounts for around 23% of the NPP compared to 17.2% in NoFe (Figure 3). Warming increases production in the ETP in FeDyn and FeMask, but also particle remineralization and microbial-loop nutrient recycling. The net effect is a decrease in export production in both models (5.1% for FeDyn and 0.8% for FeMask, Figure 3). The optimized microbial loop constant is larger in FeDyn compared to FeMask (Table B1), and this causes a larger decline in export production in the ETP in this model (Figure 3). In both iron models there is a reduction in ODZ volume in the ETP (Figures B7h and B7i). This different and decreasing trend in ODZ volume in the models with iron has cascading consequences for downstream nitrogen cycling across the tropical Pacific. Reduction of the ODZ volume leads to a decline in total denitrification, which causes water with a relatively greater nitrate-to-phosphate ratio to advect westward. This excess nitrate in the tropical central and western Pacific tends to suppress N 2 fixation by diazotrophs. However, because surface warming also enhances the N 2 fixation rate (see Equation B3), it contributes to the maintenance of N 2 fixation fluxes in FeMask and FeDyn, despite the declining of total denitrification. This decoupling also occurs because according to our model formulation, diazotrophs will always fix some nitrogen (see . Surface (0-130 m) Surface N cycle (number before + or − sign) and its trend (number after + or − sign) in the Eastern Tropical Pacific (Longitude: 115°W:60°W, Latitude: 8°S:20°N) between years 1800 and 2100. The unit for inventory is Tg C and for fluxes is Tg C yr −1 . The thickness of the arrow is approximately proportional to the respective flux. Colors denote different models: NoFe (red), FeMask (green), and FeDyn (blue). Note that all temperature dependent parameters are kept constant across all models so that the different model responses toward warming are the result of different nutrient pathways after model calibration. The parameter list can be found in Table B1.
Equations B1-B3). Therefore, increasing diazotroph primary production (the drivers of which are discussed in the following section) enhances the N 2 fixation rate in FeMask and FeDyn despite the increasingly unfavorable N:P ratio.
Strong parallels can be found between our results in the ETP and another recent model study: Tagliabue et al. (2020) found that sustained iron limitation in the ETP mitigates future declines in tropical Pacific NPP under global warming, whereas a switch to nitrate limitation could result in large declines in productivity. The driver they identified is different than what we find in our model; they describe a positive feedback in which stratification reduces nitrate resupply, which reduces biomass at the surface and allows deeper light penetration. Greater light levels at depth push the depth of the chlorophyll maximum deeper into the subsurface Equatorial Undercurrent, thereby removing additional nitrate from upwelling water and further reducing surface productivity . However, if iron remains a stronger control than nitrate on phytoplankton growth (in their experimental setup, achieved by prescribing different biological Fe uptake parameter values), this positive feedback is mitigated by the sustained addition of iron off the South American continent . As already described, our experimental setup introduces additional differences in nutrient pathways (by different parameter values) beyond sensitivity to biological iron uptake. Still it is interesting to note that even the global NPP trend in our models demonstrates increasing sensitivity to climate change with decreasing ETP iron limitation (comparing Figure 1 and Figure B5), which suggests that controls on phytoplankton growth in this region act as a "gatekeeper" for surface nutrient availability to an extent which can dominate globally averaged responses to warming. This may be a point of concern, as both iron-limited models (as well as others, such as the hand-tuned Nickelsen et al. [2015] and Kvale et al. [2020] variants, as well as CMIP5 models; W. Fu et al., 2016) broadly predict increasing retention of iron in the upper ocean with warming. Regionally, changes in primary production limitation factors into the future in the ETP can produce a cascade of divergent consequences across the tropical Pacific. In Tagliabue et al. (2020), the consequences for NPP and higher trophic levels are described. Here, we explore the consequences for nitrogen cycling. Also similar to our study, Tagliabue et al. (2020) found increasing discrepancies in transient model behavior, and drew their conclusions using a suite of models which performed similarly with respect to NPP over the historical period.

Decoupling of N 2 -Fixation and Diazotroph Primary Production in Models With Iron
In the NoFe model, we find that increasing N 2 fixation correlates with an increase in diazotroph primary production, which is intuitive as diazotrophs fix N 2 to support their growth when NO 3 is in short supply. In the tropics (30S:30N), global warming increases the temperature-dependent phytoplankton (both diazotrophs and ordinary phytoplankton) potential growth rate and the detritus remineralization rate by around 40% from preindustrial to year 2100, which accelerates surface nutrient cycling and promotes phytoplankton growth in oligotrophic waters. Simulated growth rates in the tropical upwelling zones increase faster under global warming than simulated grazing rates (that are capped at 20°C, see Equation 28 in Keller et al. [2012]), which leads to higher detritus production and then an increase in particle export out of the nutrient-replete upwelling region despite the enhanced remineralization (Figures 3 and 4a), which decreases water column oxygen and increases denitrification (Figures 5a-5c), producing increasingly elevated P* (Figures 5d-5i) conditions downstream that favor diazotrophs ( Figure 6a) and promote N 2 fixation ( Figure B4). However, in FeMask and FeDyn, diazotroph primary production increases over time but N 2 fixation does not. This is because diazotroph primary production and N 2 -fixation can be decoupled (Knapp et al., 2012). In our models, diazotrophs can utilize NO 3 when it is available. Similar to NoFe, surface nutrient cycling is enhanced by warming in the oligotrophic regions for FeMask and FeDyn (including diazotroph primary production). However, persistent iron limitation in the tropical upwelling zones ( Figure B5) mitigates the particle export response to regional warming ( Figure 4; the net export in the three ODZ regions declines by 66.5 Tg C yr −1 for FeMask and 69.0 Tg C yr −1 for FeDyn by year 2100). In line with these trends in particle export, a net contraction of the ODZ volume occurs in both models, which reduces total denitrification (Figures 2c, 2d, 5b, and 5c). This decline in denitrification over time results in reduced P* (Figures 5d, 5e, 5g, and 5i) for the downstream surface ocean, which tends to suppress N 2 fixation associated with the temperature-induced increase in the growth rate of diazotrophs in the stratified, oligotrophic regions. The strength of this effect is subject to model parameter differences. The sensitivity of the modeled denitrification response to changes in water column oxygen depends on the prescribed O 2 :N ratio, which was determined in the parameter calibrations (Yao et al., 2019) and varies between 9.54 and 10.30 among the models (Table B1). Despite the loss of a competitive advantage over regular phytoplankton due to a declining P*, diazotrophs are able to maintain a competitive advantage in the oligotrophic regions due to a relative aversion of zooplankton for diazotrophs as a food source (a top-down control; Landolfi et al., 2021). Figures 6a-6f show that the increased diazotroph primary production is usually associated with a higher effective grazing rate on ordinary phytoplankton (more than 0.2 d −1 higher compare to the rate on diazotrophs) in the region in the year 2100, for example, the Indian Ocean, the Western Tropical Pacific, and the Tropical Atlantic. This top-down control on ordinary phytoplankton and the temperature-driven increase in growth rates allow diazotrophs to enter regions that are relatively abundant with NO 3 , where they do not have to fix all the N they need to grow (Figures 6g-6i). This also contributes to the decoupling of the trends of diazotroph NPP and N 2 fixation.

FeDyn Versus FeMask-Is FeMask a Pragmatic Choice for Climate Projections?
We find no significant differences in the global N 2 fixation trends between the FeDyn and FeMask models over the time period simulated. This similarity in model results may be due to relatively similar parameter values in the two calibrated model versions, or it may be a result of compensating effects keeping surface iron concentrations relatively unchanged under global warming in the dynamic iron model (e.g., iron scavenging and remineralization). While warming enhances remineralization and recycling of iron from particles, iron scavenging and colloid formation can quickly remove the recycled iron from the dissolved phase, which keeps the simulated change in iron concentration under global warming small in the ETP (a 0.4% increase). This small change is not able to qualitatively change the modeled export production in the upwelling zone, hence there is little impact on denitrification and N 2 fixation trends. However, we do see regional differences in other quantities, such as export production trends in the Southern Ocean (Figure 4), which suggest that model-model differences may be regionally significant. In addition, FeMask has prescribed surface iron concentrations, but surface iron concentrations are expected to increase with climate warming (W. Fu et al., 2016), including UVicESCM variants like FeDyn. These dynamic changes are not represented in FeMask, hence the FeDyn model is recommended for climate projections where regional trends are of interest. However if computer resources are limiting and if the focus is on the tropical and subtropical oceans, we can recommend a transient mask of surface dissolved iron concentrations, derived from a single integration of FeDyn, could be applied to FeMask for a physically consistent, but computationally cheaper, representation of iron limitation.
Despite all calibration efforts, it is worth noting that FeDyn also has deficiencies, the most prominent being its poor skill in reproducing observed iron concentrations (Yao et al., 2019). Considering the globally averaged 0.57year turnover time of iron in FeDyn (compare to other models ranging from a few years to more than 500 years;  Tagliabue et al., 2016), the dissolved iron concentration depends strongly on local iron sources, for example, dust deposition, sediment release, and hydro-thermal activity. Sediment release is the largest iron source in our model and a better parameterization for the sediment release might be beneficial for reproducing observed iron concentrations. A recent laboratory study showed that certain types of diatoms can reduce their iron requirement under higher temperatures and maintain their growth rate (Jabre & Bertrand, 2020). This finding suggests a more flexible, even a temperature dependent iron demand. Appropriate parameterizations should be tested in future studies. Likewise, the Fe dust deposition is projected to increase in the future (Hamilton et al., 2020), which we did not account for in this study. The more straightforward option to include more realistic parameterizations of climate-sensitive iron dynamics as well as time-dependent external Fe sources, might shift the balance between the two iron models in favor of FeDyn. In the top three panels, the contour interval is 4 gN m −2 yr −1 , and in the panels (d-i) the contour interval is 0.1 mmol m −3 . While PO 4 and NO 3 were used as data constraints in the prior model calibration study (Yao et al., 2019), P* itself was not.
There are still other factors beyond iron parameterization which might further improve the representation of nitrogen cycling in our model. Wrightson and Tagliabue (2020) compiled a list of non-iron-model related potentially significant drivers of N 2 fixation, including CO 2 fertilization of diazotrophs, and a possible upper thermal threshold for optimal growth for phytoplankton, which are not yet implemented in our model and should be explored in the future. Figure 6. The change in diazotroph NPP (top, gC m −2 yr −1 ), the preferential grazing pressure on ordinary phytoplankton (middle, d −1 ), and the percentage of diazotroph growth depending on nitrogen fixation (bottom, %) in three models. The contour lines represent the change in diazotroph NPP in all panels, and the solid thick contour line denotes 0, and the solid (dashed) contour denotes positive (negative). The contour interval is 8 gC m −2 yr −1 . The preferential grazing pressure on ordinary phytoplankton is calculated as the effective grazing rate on ordinary phytoplankton minus the effective grazing rate on diazotrophs. The percentage of diazotroph growth depending on nitrogen fixation is a number between 0 and 100, where 0 means that diazotrophs only take up nitrogen from surrounding water for their growth and 100 means that diazotrophs use solely newly fixed nitrogen for its growth. Note that 0 percent dependence on N 2 fixation cannot be reached due to the formulation of diazotroph nitrogen uptake (see Equation B1).

Conclusions
The results presented here demonstrate divergent trends in simulated future global N 2 fixation and denitrification that are determined by the representation of bottom-up and top-down phytoplankton growth limitation in tropical upwelling regions (schematic in Figure 7). When iron limitation is neglected, warming induces positive feedbacks of enhanced regional export production and denitrification that triggers enhanced N 2 fixation across the central and western basin. When iron limitation is accounted for, the response of export production to warming is muted and even produces a net decline in denitrification in the ETP, leading to essentially no apparent change in global marine N 2 fixation.
While the simulated responses of N 2 fixation diverge between the three models tested, all models simulate an increase (albeit, of variable magnitude) in diazotroph primary production. That the N 2 fixation trends hinge upon the treatment of iron limitation in the tropical upwelling zones, but simulated primary production of diazotrophs does not highlights the importance of model structure for simulated future trends. This finding goes beyond J. K. Moore and Doney (2007) in that our simulations do not explore global responses to the presence or absence of diazotroph iron limitation, but instead consider a suite of ecosystem pathway differences in models that perform similarly in steady-state. By using the same physical base model, we are able to isolate and compare ecological pathway responses (which is not possible using multi-model ensembles, e.g., Wrightson & Tagliabue, 2020). Our analysis suggests flux and rate data may be just as (or even more) relevant for calibrating models as are data of biogeochemical tracer concentrations. It also implies that a different choice in calibration objective (say, application of a mix of flux and nutrient distribution data), or weighting of surface versus deep ocean attributes, or a different choice of parameter sets, could produce alternative optimal parameter sets with unique transient behaviors. However, we expect that even with different ecosystem structures, the same basic relationship between ETP iron limitation and global N 2 fixation response will apply, as Tagliabue   The upwelling ecosystem is controlled by a combination of top-down grazing and light limitation when iron does not limit primary production (left side). As the temperature increases, upwelling is reduced ("−" in the center of the figure). In the scenario not considering iron (left part of the figure) primary production increases faster than grazing pressure, which leads to increased particle export (+) and higher remineralization (+). Higher temperatures and increased remineralization intensify water column ODZs (+) and associated denitrification (+). Elevated denitrification leads to higher N 2 fixation, and an increase in diazotroph primary production downstream (+). However, including iron limitation (right side) mitigates the export production response to warming in the upwelling ecosystem. Export production declines (−), and the ODZ volume also shrinks (−). More oxygen in the water column leads to less denitrification (−) and reduced P*, which suppresses N 2 fixation downstream (−/+, no apparent change) despite temperature-driven increases in diazotroph primary production. effect in their model ensemble for tropical NPP. An interesting additional test, which we leave to the future, would be to explore the role of model parameter choice in long-term recovery from the decoupling between diazotroph primary production and nitrogen fixation.
Our finding of a strong connection between global N 2 fixation rates and iron limitation in upwelling zones sheds new light on two controlling factors of N 2 fixation, namely iron limitation and denitrification. From a global view, spatial patterns of N 2 fixation are not only controlled by iron availability in oligotrophic regions (as proposed by Ward et al. [2013]), but also by the iron limitation upstream, which regulates the denitrification (and hence, P*) pattern. Lastly, our study demonstrates a role of temperature in determining N 2 fixation perturbation patterns. In our earlier study we demonstrated a correlation between denitrification and N 2 fixation as well as globally balanced rates in all model steady states (Yao et al., 2019). However, this linkage is reduced in transient conditions, such as ongoing climate warming. As temperature increases, temperature-sensitive rates regulating bottom-up and top-down controls on export production respond differently. The net result is a decreasing sensitivity of N 2 fixation to denitrification (seen particularly in the models with iron) that can to lead to at least a temporary imbalance in global rates. Furthermore, there is growing evidence that heterotrophs may play a role in N 2 fixation (Moisander et al., 2017), and that N 2 fixation in key cyanobacterial diazotrophs may not be stimulated by the lack of NO 3 (Mills et al., 2020). Both of these factors might contribute further to the spatial decoupling between N 2 fixation, and denitrification, and the temporal and spatial decoupling of N 2 fixation and diazotroph primary production which might operate on decadal and longer timescales (J. K. Moore & Doney, 2007). The long-term implications of potential N 2 fixation and denitrification imbalance would be an interesting topic worthy of further study.
This study has shown that even models that are calibrated essentially equally well to global compilations of biogeochemical tracer distributions in an assumed seasonally cycling climatological state, can lead to quantitatively and qualitatively different projections of key biogeochemical processes like nitrogen fixation in transient global warming simulations. Until there are sufficient observations to constrain the sensitivity of transient model simulations to environmental change, it will remain difficult to estimate uncertainties of model projections. Further developing, testing and applying mechanistically sound representations of relevant model features, such as the marine iron cycle and marine nitrogen fixation can be a promising way forward to reduce uncertainties. This may work best when the calibration process can be extended to include large-amplitude historical variability that may resolve part of the sensitivity of marine biogeochemistry to changes in the environment.

Appendix A: Nutrient and Light Limitation on Phytoplankton Growth
The nutrient limitation factor is defined in the model as a function of nutrient concentration and phytoplankton nutrient uptake half-saturation rate: X is NO 3 , PO 4 or Fe, C X is the concentration of X and k X (O or D) , is the phytoplankton half-saturation of the uptake rate for X. Light limitation is derived from Keller et al. (2012, for NoFe and FeMask) and Nickelsen et al. (2015, for FeDyn) as: where α (O or D) is the initial slope of the photosynthesis-irradiance curve (in FeDyn, α (O or D) is dependent on Fe Lim(O or D) , details see Nickelsen et al., 2015), I is the in situ irradiance and J max (O or D) is the in situ maximum potential growth rate.
a is the phytoplankton growth rate constant at 0°C. F temp (O or D) is a temperature factor for biological process.
T is the situ temperature and T b is the e-folding temperature of biological rates, C d is a handicap for diazotroph growth and the 2.61 temperature scaling is to inhibit diazotrophs growth below 15°C.
The iron is not incorporated in the nutrient minimum function but as a factor that directly governs the growth rate. The in situ growth rate of ordinary phytoplankton (J O ) and diazotrophs (J D ) are defined as: More details of model equations can be found in Keller et al. (2012, FeMask and NoFe) and Nickelsen et al. (2015, FeDyn).

Appendix B: Nitrogen Fixation by Diazotrophs
As described in Schmittner et al. (2008), the diazotrophs can take up nitrogen from the surrounding water when it is available, and it can fix nitrogen when the preformed nitrogen does not satisfy the demand for growth. The growth rate for diazotrophs is described in Equation A7, and does not depend on preformed nitrogen concentration in the water. The nitrogen uptake rate from the surrounding water is described as in Schmittner et al. (2008) as (Table B1; Figures B1, B2, B3, B4, B5, B6, and B7):  Note. For computational and conceptual reasons, the number of tuneable parameters had to be constrained to a relatively small number. We selected those parameters that showed largest impact on the model results as shown in Yao et al. (2019). The microbial loop constant describes the rate of bacteria (not explicit modeled) recycling the phytoplankton biomass back to the resolved nutrient pool. The diazotrophs' handicap means that the diazotrophs maximum growth rate is around 40% of ordinary phytoplankton (a). For more details of the parameter of the model, please see Keller et al. (2012) and Nickelsen et al. (2015). a Parameter value optimized for the model. where C NO3 is the nitrate concentration in the water, and k N is the half saturation constant for N uptake. Hence the N 2 fixation rate is: More details of the model equation for N 2 fixation can be found in Schmittner et al. (2008).

Data Availability Statement
Model code and model data used in writing this manuscript are available on OceanRep-GEOMAR Repository (https://data.geomar.de/downloads/20.500.12085/673e7de0-20ab-4dd3-afe9-c4bfb00b1faf/). Figure B7. Annually averaged total NPP, export, and oxygen concentration (depth 350 m) in the Eastern Tropic Pacific. Color shading represents the difference between year 2100 and 1800 and contours show simulated values in year 2100. The unit for NPP and export production is g C m −2 yr −1 and the unit for O 2 concentration is mmol m −3 . The interval between neighboring lines are 50 g C m −2 yr −1 for NPP, 2 g C m −2 yr −1 for export, and 5 mmol m −3 for O 2 .