This manuscript assesses simulations with five different global ocean biogeochemical models to assess uncertainties about the role that biological dinitrogen (N2) fixation will play in changes in ocean productivity and biogeochemistry in a warming climate. Generally, the analysis is sound, but the presentation is weak in places. I offer some suggestions below for ways to revise the paper to make the overall conclusions more compelling.

(1a) I think the comparison between IPSL-CM5A-LR and the offline simulations is presented in a misleading way. This is not an apples-to-apples comparison. An offline simulation will never reproduce the parent model exactly. Therefore the suite should really include PISCES-v1 run in the offline mode. Possibly the differences from the ESM would be small. But it would be nice if the reader could verify that. Consider cumulative CO2 uptake (Table 3). It differs by only 16 PgC between the two ESMs despite the massive increase in N2 fixation in CM6A (and different atmosphere and ocean models, and different emission scenarios). But all 3 offline simulations differ from CM5A by 50-90 PgC, despite the 'identical' (328) circulation. I suspect this has more to do with the inline/offline configuration than with the bgc model structure. If the authors want to test this, I think they could achieve this with an offline PISCES-v1 experiment.

(1b) I find some of the text explaining the results shown in Figure 5 vague or misleading. Some of the results are clearly robust to the differences in the physical ocean environment. For example, when we compare offline PISCES-v2 to CM6A the change in total N2 fixation is very similar (Figure 5b). But the net change in NPP differs by almost a factor of 2 (Figure 5a), and I think this glossed over in the text (316-321). The explanations of the differences among the offline models are also vague in places (e.g., 306-309).

(2a) I don't think the analysis of mechanisms underlying the differences among models is as well presented as it could be. On 341-343 the text seems to be saying that N2 fixation declines due to warming in IPSL-CM5A-LR, which should not be the case, and is unaffected by DIN-inhibition. But the latter seems to be contradicted in the very next paragraph (348) and by the data in Figure 6. In the ESM L_N only goes down to ~0.8. But in the other models it doesn't change at all. So if 348 isn't referring to the ESM, what is it referring to? Similarly, in Section 3.5 I think there are several assertions that are questionable and not really supported with data. It isn't obvious to me why we would have a >2X larger decline in subsurface O2 in IPSL-CM6A-LR than IPSL-CM5A-LR, when export production is about the same and shows a very small (and negative) trend in IPSL-CM6A-LR (Table 3). How does an increase in N2 fixation result in subsurface O2 depletion without affecting export? Possibly via DOM, but that is not substantiated or even discussed. And why would we assume that it is due to remineralization and not to circulation given the different climate models used? As in 3.3, I find the presentation here a bit careless.

I also do not think that the Western Pacific box in Figure 6a is as representative as is implied. If we take the North Pacific subtropical gyre as whole, the southwestern corner is probably the warmest part. At HOT, for example, annual mean SST is around 25C. So with the Breitbarth temperature function N2 fixation would initially increase, although it would decline if net warming exceeds ~2C. In the Western Pacific box SST exceeds 26C right from the beginning, so the decline is monotonic (Figure 6cd).

I also don't think including N* in Figure 6 is useful. N* is never mentioned up to this point (e.g., the Methods do not mention any of the models as having a dependence of N2 fixation on N*). Anyway, why would N* be negative at the surface in a region where a lot of N2 fixation occurs? Possibly this could be useful for evaluating the realism of the models wrt N/P stoichiometry, but at present it adds little to the analysis.

(2b) Similarly, I'm not sure including N* in Figure 7 and Table 2 is a good idea. When I look at the observational distribution I think I understand why negative values occur in nutrient-depleted surface waters, but this is probably an artefact of the kind of data used. In the subtropical Atlantic, for example, Gruber and Sarmiento (10.1029/97GB00077) calculated that N* is positive in the subsurface waters and therefore that there was probably net N2 fixation in the overlying surface waters. But in this plot we see a broad, contiguous area of negative N*. I think this is an artefact for two reasons. (1) In nutrient depleted surface waters almost all the nutrients are recycled. Therefore much of the DIN may be NH4, which is rarely measured. But P has no such redox chemistry: recycled PO4 and 'new' PO4 are the same, analytically. So using the gridded NO3+NO2 data product underestimates DIN and creates an artificially low N/P. (2) Concentrations are often lower than the analytical detection limit (ADL) for standard methods. Because the ADL for PO4 is not 16X lower than that for NO3, this again creates an artificially low DIN/DIP ratio. Why this also happens in the model is interesting. Do the model N* estimates include NH4?  (Also, why are the columns ordered differently in Figure 7 and Table 2?)

(3) I general I do not think that the cluster analysis is adequately explained. Figure 4 is quite a lot of information to digest, and the presentation could be improved. Firstly, I would suggest that Figure 4(c) be moved to the top, as it contains the definitions of the colours. I think the black "n.c." segment on the colour bar should be removed (the only black on the map is over land), and the threshold for no change should be stated. The text states that  cluster 4 – pink indicates that "growth rates increase without any significant change in N-fixation" (264-265) but does not state what is the criterion for a change to be considered significant. For panel (a) I would consider (1) using common ranges for the x and y axes, (2) stating in the caption that the greyscale applies (equally?) to all colours, and (3) making the 'no data' squares white instead of black. If using a common x and y range would interfere with the visual presentation, consider including it as an additional Supplementary version. I would change the caption to Figure 4(c) to something like "Global maps of the distribution of the five clusters". (Note that in my copy there is abcd in the caption, but no actual labels on the subplots. Also the labels in the figure all say \mu but the caption says m; this might just be due to PDF rendition.)

Conceptually, the description of the cluster analysis methodology is not very clear or specific. The only literature reference is to a 771 page textbook from 1991. So presumably this is a well-established methodology, but it will be unfamiliar to many readers, and I do not think the explanation given in the text is very illuminating. Nor can the reader easily trace it back to its cited source. Does P represent a probability? Probability of what? That X or Y will fall within a specific increment within (-1,1)? That it will be positive, negative, or NSDZ? What do X and Y represent? Normalized anomalies at the individual model grid points? Remapped to a regular grid? If it is point-by-point on whatever grid (Figure 4c+d), what defines a probability? The anomaly at each point has a unique value; sampling over some range of inputs is required to generate a probability distribution. Sorry I'm just not following exactly what was done here. A few sentences of explanation can go a long way in helping the reader to understand what is being presented.

(4) I find Section 3.4 a bit disjointed. It almost feels like it is two separate sections spliced together and might better be split in two. 375-379 is like a wrapping-up of one topic and then a whole new one is introduced. It might also be a good idea to combine the paragraphs on 394-398 and 405-410 into one, so that the general background on emergent constraints leads directly to the application that is directly relevant here.

Methods

In 2.1 I would add a few sentences about the setup of the physical ocean models. On 126 it is stated that NEMO consists of ocean dynamics, sea ice and biogeochemistry. But there are many options for various physical process parameterizations and some configuration of these was used in the 'frozen' versions used for CMIP5/6.

It would be good idea to state something in the Methods about how advection and mixing are done in the offline simulations.

The description of the modified N/P stoichiometry is a bit confusing. On 165-168 it is stated that using a high N/P for organic matter derived from diazotrophy is a change in PISCES-quota from the base case. But on 156 is says that PISCES-v2 has the same N/P=46 for organic matter derived from diazotrophy. It appears that the only difference is that in the base case the organic matter is implicit.

Some specifics:

26 change "model" to "models"

56 CESM2s is an official model name?

65 delete "atmospheric" (see also 94)

97 "unconstrained" misspelled (see also 359, 502)

160 change "advanced" to "modified"

206 "phytoplankton-realised growth" do not hyphenate

253 change "ocean physics" to "physical ocean"

267, 269 change "independent from" to  "independent of"

270 '' the reinforced relationship" Odd choice of words. How about: "confirms the intensified effect of diazotrophy on phytoplankton growth ..."?

280 change "oceanic" to "ocean" (actually I think all 4 occurrences of "oceanic" could be changed to "ocean", but this one in particular)

297 "The increase in N-fixation is dampened everywhere as compared to PISCES-v2, with even regions where N-fixation decreases in PISCES-v2fix" How about: "The increase in N-fixation is small everywhere compared to PISCES-v2; there are even regions where N-fixation decreases in PISCES-v2fix"?

303 change "only slightly increases" to " increases only slightly"

362 "better performance scores" I would not use this term unless the 'score' is defined somewhere.

467 change "under a similar high-emission scenario" to "under similar high-emission scenarios"

468 "additional" misspelled

471 "of the of the"

479 change "participated to" to "participated in"

481/496 Riche and Christian: 2017 or 2018? (see also 98)

485-494 It might be a good idea to mention Pahlow et al and Inomura et al in this paragraph, as they were described on 423 as presenting "more mechanistic models of N fixation".

Figure 7 is too low resolution: when I expand it to be readable, the fonts are blurry

This study by Bopp et al. investigates how changes to N2 fixation affects NPP in a suite of climate model simulations (IPSL) with PISCES biogeochemistry under historical and future warming scenarios. They simulate 5 different model versions, 4 of which apply different formulations for biogeochemistry, most notably N2 fixation, and 1 simulation with different ocean physical circulation and resolution. The simulation with higher resolution (IPSL-CM6A) did not significantly affect N2 fixation compared to the simulation with similar biogeochemistry and lower resolution (PISCES-v2) so they mainly focused this study on the differences in N2 fixation parameterizations.

They implemented a variety of different N2 fixation parameterizations including phosphorus limitation, temperature-dependent growth, elemental stoichiometry, and underlying biogeochemistry (the latter not well described in the paper). Their model simulation with a strong increase in N2 fixation (PISCES-v2) caused a strong increase in NPP, whereas the other simulations with a slight increase or decrease in N2 fixation projected a decrease in NPP throughout the 21^st century. They focus on a region in the western tropical North Pacific to better understand the mechanisms responsible for the difference in these simulations. The one clear driver for driving the high N2 fixation rates in PISCES-v2 was the exponential temperature-dependence growth rate, whereas the PISCES-quota simulations using the Breitbarth et al. 2007 bell-shape reduced diazotroph growth rate at high >26°C temperatures.

Overall I think this is a fine manuscript that highlights the importance of N2 fixation on NPP, which is often neglected or not analyzed in model simulations.  Another aspect I liked was how the authors distinguished between NPP and phytoplankton biomass instead of only focusing on NPP like most studies. My one major issue with the manuscript is the insufficient description of which processes are causing N vs. P vs. Fe limitation and thus driving the N2 fixation response in the different biogeochemistry model simulations, but this can be addressed in revisions. This information is included in some of the references, but given the emphasis on N2 fixation in the paper I think some additional details should be included and discussed in this manuscript as well.

Major Comment

The one critical issue I do not understand in this study is why N2 fixation in IPSL-CM5A-LR does not increase more similarly to IPSL-CM6A-LR. Since they have the same temperature-dependent growth rate, the authors state this is “due to the Ln term (limitation of excess inorganic nitrogen)” (line 339). However, all N2 fixation parameterization contain this Ln term (Figure 2). Why does Ln only decrease in IPSL-CM5A-LR and not any of the other simulations?

My guess is that it is caused by some changes in the underlying biogeochemistry (e.g. iron limitation or denitrification), which is only mentioned in one brief sentence in the description section 2.2 (lines 143-145) and not further discussed. The processes that contribute to the decreasing Ln term in IPSL-CM5A-LR but not the other simulations should be described and discussed to understand these results presented here.

Minor Comments:

Lines 144-145: “… external sources of nutrients and the treatment of water-sediment interactions”

Following my comment above, these processes can significantly impact N2 fixation rates. While I understand they are described in previous papers, the important processes controlling N2 fixation in these simulations presented here should be described and discussed.

Line 155: “… phosphorus is added to the phosphate pool”

I thought diazotrophs are responsible for consuming PO4 and DOP to very low levels e.g (Mather et al., 2008), not providing a source. These fluxes are generally small among the entire circulation-biogeochemical system so it is likely not significant, but I am concerned that this PO4 source could be supporting additional NPP in an unrealistic way.

Line 158: “… annual restoring of global mean PO4 concentration …”

I am surprised to see a PO4 restoring term in a prognostic climate model. I assume this is applied only during the spin-up and not historical/future projections? Does this occur at all depths?

Line 306: “The regional changes changes of NPP resemble those simulated in IPSL-CM5A”

Is this also true for N2 fixation?

Lines 307-309: “… differences between PISCES-quota and two other offline simulations originate from the major developments in PISCES-quota such as variable C:N:P stoichiometry and the inclusion of a third phytoplankton …”

This comes back to my major comment. It is interesting that newest (PISCES-quota) and the oldest formulation presented (PISCES-v1 in IPSL-CM5A-LR) give very similar projections of changes to NPP and N2 fixation (Figure 5a,b), despite their large differences in biogeochemical model structure/parameterizations, whereas the in-between version (PISCES-v2) shows a significant difference. I would be interested to know which changes mainly contributed to the large increase in N2 fixation and NPP and PISCES-v2 compared to v1.

This statement specifically mentions stoichiometry and the additional phytoplankton type, but I thought it was pretty clear that the change in temperature-dependent N2 fixation growth rate was responsible for these decreases in N2 fixation and NPP in PISCES-quota compared to v2 as described in the following section so I am confused by this statement. Please describe how variable C:N:P stoichiometry and a third phytoplankton class causes higher N2 fixation because this is not clear. I think this emphasizes the need to better understand the key processes causing the N2 fixation differences in the different simulations.

Mather, R. L., Reynolds, S. E., Wolff, G. A., Williams, R. G., Torres-Valdes, S., Woodward, E. M. S., Landolfi, A., Pan, X., Sanders, R., and Achterberg, E. P.: Phosphorus cycling in the North and South Atlantic Ocean subtropical gyres, Nature Geoscience, 1, 439-443, 10.1038/ngeo232, 2008.