General comments:

Stappard and co-authors detailed the NutGEnIE biogeochemistry module extension to an existing Earth system model. Their motivation was to create a modelling system capable of extended simulations in order to investigate controls on ocean primary productivity over long timescales. While the model would represent a substantial contribution to the field, I found several issues with the manuscript in its current form, which I’m detailing below, that warrant addressing by Stappard et al. before continued evaluation. In particular, the methods section should be heavily revised to improve clarity, and the discussion section should be greatly expanded to include a more thoughtful analysis of the model biases as compared to observational products.

Specific comments:

(1). I’m still uncertain of the total number of explicit biogeochemical tracers the model carries. Typically, model development manuscripts include a full list of explicit tracer equations (often in the supplemental material). I recommend this exercise here to improve clarity, especially as it pertains to the iron cycle in the model.    

(2). Upon first read, I was confused how POM was treated in the model, since it was discussed heavily in Sections 2.3.1 to 2.3.3. It was not clear until page 8 of the manuscript (Section 2.3.4) that it is implicitly represented. I recommend a reorganization of section 2.3 to (1) first describe the equations for nutrient uptake, (2) its portioning into DOM/POM, and (3) the remineralization scheme of DOM/POM.

(3). It is unclear why NutGEnIE necessitates a burial fraction parameter (k_BF). Shouldn’t some portion of the POM flux make it to the bottom of the deepest grid cell, and shouldn’t that represent the portion that is buried? Explaining the rationale for such an enhanced burial flux is needed.

(4). The authors briefly introduce top-down control of autotrophs by grazers in Section 1.2, yet do not discuss any model caveats by ignoring this process in NutGEnIE. In particular, this can be an important control on surface primary production in HNLC regions (i.e., some suggest it is ‘bottom-up’ via iron limitation, while others suggest top-down controls lead to higher observed surface nutrients compared to other regions). This caveat should be outlined in greater detail within the discussion section, especially as it relates to biases in model performance (e.g., surface nitrate and phosphate comparisons in HNLC regions).

(5). I recommend elaborating more, or at least clarifying, the iron cycle in the model. Including the equations for each explicit iron tracer would especially helpful in clarifying the underlying dynamics of iron within the model.

(6). L. 184 - 185: Based on the wording here, it is unclear if anammox is represented in NutGEnIE at all. In Figure 2., it does not seem represented. If the authors choose to omit anammox, this should be discussed as a caveat to the model, since anammox is responsible for a considerable portion (~28% based on remineralization stoichiometry, see Babbin et al. 2014 - Science) of fixed-N removal in oxygen-minimum-zones.

(7). L. 224: The authors briefly mention ‘particle concentration’ (Cp) here, but it wasn’t clear to me what they meant by this. POM is only implicitly represented here, so how do they quantify particle concentration?

(8). Eq (12 - 14): I’m confused by these equations. For example, in Eq 12, this is just a fraction (Michaelis-Menten function ranging from 0 - 1), not the actual rate of aerobic remineralization. It would be helpful to first define the total depth-dependent remineralization rate as the divergence of POM flux in each grid cell (i.e., R_remin) . Then, for example, equation 12 would be better represented as [O2]/(K_O2 + [O2]) * R_remin. Also, it is unclear how the sum of the rates in Equations 12 - 14 equals the total remineralization rate of POM.

(9). Eq (15 - 17): Is particulate organic iron (POFe) not also remineralized in the model? Also, it is unclear where and how the Gibbs free energy yield values are used (I also recommend adding subscripts to the distinct Gibbs values, e.g. DeltaG^{o}_{O2}). Are they folded into the calculation of the inhibition constants? If so, it would be important to include their formulas, either here or in the supplementary material. 

(10). Figure 4: I believe the authors are missing the nitrification source term for the nitrate panel. Also, here they specify nitrogen fixation as an input to the nitrate budget, but it is not represented in the tracer equation (Equation (4)).

(11). In general, text below equations are often missing units for parameters. It would helpful to include the units, and reference tables when mentioning parameters throughout the methods section.

(12). In Section 2.3.7, only Fe is represented in equations (9) and (10), yet the authors mention that complexed iron (FeL) is also available for biological uptake. Did they mean to write FeL in equations (9) and (10)? If so, wherever the authors mention Fe in the text, should they instead write FeL for clarity?

(13). In Section 2.3.7, it is unclear if ligands are an explicit tracer. If they are, then Figure 3 suggests particulate iron can be created in depth cells below the euphotic zone, yet equation (11) states POM flux is only set by euphotic zone values. Please clarify.

(14). L. 459: There are missing details regarding the initial conditions of the model runs.

(15). L. 477: “ocean physics is not the focus here, so the properties are not discussed in detail”. While I generally agree with this statement, there are some notable biases (compared to WOA, treated here as reality) in both temperature and salinity that could cause stratification and other errors in NutGEnIE. For example, in the Pacific, there is anomalously low temperature in between the subtropics (Figure 5c). Would that not have an influence on the delivery of deep nutrients to the surface there? For instance, there appears to be a very similar bias distribution for surface phosphate and nitrate (Figures 7c and 8c). Similarly, the transect biases in temperature match the patterns in the biases of nutrients. For example, the model appears to be too warm in the deep Atlantic, and is too low in PO4 and NO3. In the deep Pacific, the model is too cold, yet also high in PO4 and, to a lesser extent, NO3. The authors could potentially strengthen their validation exercises by relating some of the nutrient biases to stratification differences.

(16). Similarly, in Section 3.3.4, I believe some of these biases in surface oxygen are related to co-located biases in surface temperature. For example, Figure 5c suggests there is a cold bias in the surface of the subtropical Pacific, where there is also anomalously high oxygen values (Figure 10C).

(17). Figures 7c and 8c: These surface biases are quite large, especially considering the model is designed to study limiting surface nutrients. For example, in Table 6, the authors report a surface mean nitrate value from WOA of 6.0 umol/kg, whereas NutGEnIE reports nearly half that value (2.9 umol/kg). I would have liked to see a more thoughtful discussion on why these biases don’t impact the authors’ confidence in the model’s performance.

(18). L. 673: I’m not sure the model results support this conclusion, especially considering the model greatly underestimates both N-fixation (L. 645) and denitrification (L. 654) when compared to other studies (Section 3.1.1.). Why do the authors think that is the case in their model runs? If they believe this does not impact model performance, this should be elaborated on in the discussion.

(19). I would have liked to have seen some additional validation figures. Comparing model AOU (apparent oxygen utilization) to WOA estimates could help improve confidence in the representation of remineralization in the model.  Similarly, N* can be extracted from both the model and WOA to assess model performance in generating spatially-varying N-fixation and denitrification signatures. Finally, oxygen-minimimum-zone (OMZ) thickness comparisons (e.g., thickness of waters within each cell-profile that are less than 60 umol/kg O2) would improve overall confidence in denitrification within the model, since OMZs are crucial regions for balancing the global N-budget.

Technical comments and corrections:



L.30: “Net primary production (NPP) represents the total rate...” (suggested edit)



L.32: “phytoplankton produce biomass” (typo?)



L.37: “Nutrient supply to the euphotic zone acts as a fundamental control on ocean PP levels”. (typo, I also recommend adding a semicolon before “this supply and subsequent growth limitation”)

L.41: “grazing reduces the total amount of photosynthesis” (this could be rewritten for clarity).

L.48: “Again, studies have proposed methods of modeling temperature limitation” (I suggest rewording or merging with previous sentence for improved flow of manuscript).

L.50: “Elements (C, H, N, P, O and S)” (I suggest defining these explicitly or omit )

L.51: “such as proteins and nucleic acids.” (I suggest removing ‘etc.’)

L.55: “Addition of the proximate limiting nutrient (PLN) stimulates immediate growth” (I suggest removing the comma after (PLN))

L. 72: “Their work suggests the stratified subtropical gyres” (suggested text addition)

L.85: ‘The modelling’ should be ‘the modeling’ (no capitalization after semicolon). Or, could be rewritten as “Deutsch et al. (2007) only conducted model simulations over short timescales to a modern ocean steady-state”

L.90: Recommend removing “but at the same time”

L.95 - 99: I’m not sure this paragraph is necessary.

L. 102: Consider rewriting the sentence starting with “Such investigations” to place the information and citations outside of the parentheses.

L. 105: Remove comma after “(cGEnIE)”

L. 110: Fix bad reference to section.

Figure 1: Consider adding longitude and latitude tick labels to this and other map figures.

L. 138: “biogenically-induced chemical fluxes (ref) and is capable of” (suggest adding hyphen and ‘is’)

L. 142 - 151: Please be consistent when using chemical abbreviations here and in the rest of the manuscript. Typically, it is best to define chemical abbreviations before using abbreviations throughout the rest of the manuscript (i.e., “that include dinitrogen (N2), ammonia (NH4), calcium (Ca), and sulphate (SO4)”). Also, N2 is italicized here, when other forms in this paragraph are not. Please also include any negative or positive charges on NO3, PO4, NH4, SO4, etc, wherever they appear in the document.

L.149: “By default, results are output as annual averages for each grid cell” (suggest removing ‘figures’, since models only output numerical data, not visuals).

L.150: “giving the possibility of results output relating to shorter timeframes” (suggest rewording)

L.153 - 156: Please be consistent when using parentheses to split sentences. For example, here the authors use ‘a)’ ‘b)’ when elsewhere they use ‘(a)’ and ‘(b)’.

L. 156: “we detail the most pertinent features” (suggested edit)

L. 157: “nutrients taken by phytoplankton are instantly converted to POM and DOM in the surface ocean” (suggested edit for clarity)

L. 166: ”The nutrient uptake terms (…) only have a value in the surface layer”. Is this by design? Does the model restrict uptake to the top grid cell, or is this a result of the model? Please clarify.

L. 183: “Like denitrification, anammox converts…” (add comma, remove capitalization of Anammox)

L. 188: Please rewrite this sentence for clarity (there are some typos), and please define RP, RN, DDFe, and BBFe in the text before using their associated abbreviations. The authors could also reference the associated supplemental figures in the caption of Figure 2.

Figure 2: Please define POP, PON, and POFe in the caption, since they are labeled on the Figure panels. I also recommend redesigning these figure panels such that only explicit tracers, and the fluxes that couple them, are represented. For example, here the authors use green circles to define other phytoplankton and diazotrophs, which may confuse readers into thinking these are explicit tracers.

L. 204: Please consider adding units when describing the components of equations (i.e., mmol / m3 d).

L. 208: “bit” should be “but”

Eqs (2) and (4): How are ’S_{PO4}’ and ’S_{NO3}’ different from ‘RP’ and ‘RN’? If they are identical, consider just using one term. Also, the authors mention ’S_{PO4} being configurable by a parameter on line 210? In Figure S1, surface nutrient inputs are supplied via a forcing field rather than a parameter. Please clarify what is meant here. Also, I recommend adding supplementary figures (similar to Figure S3) that detail the magnitude of the surface forcing for both nitrate and phosphate.

Eq (4): I recommend using consistent terminology for reaction rates. For example, why use the delta symbol for nitrification and the ‘R’ symbol for denitrification, when they are both reactions in the model? Perhaps it is easier to use ‘R_nit’ and ‘R_den’ for clarity?

Eq (6); ‘B_{Fe}’ is labeled as ‘BFe’ in Figure 2. Please be consistent between figure labels/captions and equations. Also, is ‘free dissolved iron III’ another explicit tracer in the model?

L. 217: Please reference Figure S3 after mentioning ‘DDFe’, either here for after describing the re-gridding of Mahowald et al.

L. 232: It would be helpful to define units for these terms in the equation. For example, the units of V^{OPhy}_{max} are unclear since the model does not explicitly represent biomass. The authors mention the parameter values on line 255, but please move them to earlier in the text when they are first defined.

L. 234: I’m confused by the light limitation term. In Figure S37 the text mentions that “lower values of ‘K_light’ indicate that light is more limiting to nutrient uptake”. Where does ‘K_light’ fit into the equation on line 234? Also, this relates to my previous comment about light limitation. Is it restricted to the surface grid cell?

L. 239: ‘Therefore, DIN and Fe uptake are scaled’? (Did the authors mean to include ‘uptake’ here?)

L. 242: Please remove the comma after ‘atmospheric transfer’, and the authors have already used ‘N2’ earlier in the text but are just now defining it as dinitrogen.

L. 245: I recommend moving the discussion of N_thresh (Section 2.3.5) here to improve clarity.

L. 256 - 257: I was confused by this last sentence at first. I recommend reorganizing the text to specifically say that the V^{*}_{max} terms are temperature-dependent (i.e., can reach values both higher and lower than what are reported here). Then, temperature-dependent uptake is scaled by the combined limitation terms. Also, it is unclear why the authors state that “the maximum percentage of the grid cell nutrient concentration taken up by other phytoplankton each time step is 80%”. That would assume that the other limitation terms are equal to 1, but the authors state that the other maximum values are 0.7 and 0.5. Please reword or consider omitting this sentence from the text to improve clarity.

Eq (11): Please include units for the flux. I also think the authors can remove the (z = z_o) from the equation and just use z_o. Also, from this equation, it does seem like the model restricts uptake to the top cell of the model. Please detail this earlier in the manuscript so it is more clear.

L. 266: This is the first time the authors have mentioned that sulfate is an explicit tracer in the model. It would help improve clarity to have mentioned this much earlier in the methods section, and to provide a tracer equation for sulfate (and all other tracers) in the supplemental material.

L. 277: Is the model configurable to represent non-Redfield stoichiometry?

L. 286: Since bacteria are not an explicit tracer, it is not necessary to say ‘by bacteria’ here.

Eq (18): It might improve clarity to separate the two Michaelis-Menten functions rather than showing the product. Also, please include the units for maximum rate of nitrification, and consider using a different symbol (i.e., R_{nit}) to match the style of other reactions in the manuscript.

L. 301: Please reword this sentence.

Figure 3: Please be consistent with terminology elsewhere in the manuscript. Does ‘PartFe’ represent ‘POFe’ in Figure 2?

L. 310: It was confusing where the caption ends and the next sentence begins (minor comment).

L. 312: Please move this text on the equilibrium between Fe, ligands, and complexed iron to line 304 for clarity.

Tables 1 - 5: References to these tables should be placed in the appropriate locations in the text when first mentioning these parameters.

L. 340 - 344: “Parameters were adjusted to result in a combination that showed best agreement with observed nutrient distributions...” (suggested edit)

L. 404: I recommend converting +- 26% into Tg N yr to match the other estimates.

L. 405: “Wang et al. (2019) provide location of fixed nitrogen due to…”. Do the authors mean fixed nitrogen loss?

L. 420: Here the authors italicized PO4, NO3, and O2, whereas in other areas they are not italicized. Please be consistent.

Figures 5 - 8: Please provide labels on the maps and transects so that it is easier to identify which panels represent model results vs. which panels represent validation products. In all Figures, it would also help to extend the colorbar limits slightly to better represent values beyond their current ranges, since the values often reach the maximum/minimum limits. This is most notable in Figures 7 and 8, since the representation of these nutrients are a central point to the manuscript. Finally, please include latitude ticks (higher priority) and longitude ticks on map figures.

Figure 9: Can the authors please convert longitude labels from 0 - 360 format to -180 to 180 format (with E and W labels)? The inset map is also quite small and could be resized for clarity. 

L. 610: “For all variables…” (typo)

Table 6: Please use the same numerical precision between surface and interior reported values (i.e., surface value of 0.58 vs 2 for PO4).

L. 651: “Denitrification can occur throughout the water column”. Ideally, this won’t be the case. Instead, NutGEnIE should restrict denitrification to only very low oxygen regions. Perhaps reword this to be clear. Also, I suggest removing ‘by bacteria’ since the model does not resolve bacterial biomass or their metabolisms.

Figure 15: Consider rewording the caption. 

Review of "NutGEnIE 1.0: nutrient cycle extensions to the cGEnIE Earth system model to examine the long-term influence of nutrients on oceanic primary production", by Stappard, Wilson, Yool and Tyrrell, submitted to Geoscientific Model Development in 2025

Scope of the manuscript, major comments, and recommendation

-----------------------------------------------------------

Stappard et al. discuss a multi-nutrient-cycling extension to the established cGENIE ocean biogeochemical model. The extended model is simple in the sense that the biological cycling of the three nutrients phosphorous, fixed nitrogen, and iron  is strictly implicit, without explicit description of living or dead biomass. Instead, export production by non-nitrogen-fixing organisms for example is described as being proportional to the most limiting nutrient, assuming a fixed Redfield-like stoichiometry, with the proportionality factor taking into account temperature, light and nutrient effects on the maximum growth rate of phytoplankton in a way as it is done in many other models that do describe biomasses explicitly.

This simplification, and the comparatively coarse resolution, allows the model to be integrated with large timesteps, making it possible to be integrated over tens of thousand of years, longer than the residence time of phosphorus in the ocean, and hence allowing it to treat all major nutrient cycles in the ocean as an open system, prescribing just the external inputs of nutrients, e.g. from riverine input, and letting the system decide which average nutrient levels are ultimately reached in the ocean. Most other models at least treat the phosphorous cycle, some also the nitrogen cycle as closed systems, neglecting inputs to and losses from the ocean and setting a fixed average nutrient concentration.

Indeed, the main motivation that the authors give for developing this model is to investigate what is the 'ultimate limiting nutrient' in the ocean, in the sense that its inventory sets the overall strength of the biological carbon cycling. But clearly, the usage of that model needs not be limited to that rather specific geochemical question: The model could also be used how nutrient inventories in the ocean change over time when e.g. external climate and nutrient influxes change, such as happened over glacial-interglacial timescales.

The model described in this manuscript, with its simplicity and consequent speed, fills a niche at one end of the different approaches for modelling the ocean carbon cycle and is therefore a quite useful addition to the literature. Model results for present-day climate are compared to nutrient and oxygen climatologies; I especially liked that modelled iron distributions (for which no climatology exists yet) are assessed against an important subset of the GEOTRACES intermediate data product. Overall, the manuscript is well written. It fits thescope of the journal, and in the end I think it should be published in

Geoscientific Model Development.  I have three main criticisms, however, and the manuscript should be revised accordingly before being published.

The first and most important criticism is that there are several unclear points in the model description of the modelling of different nutrient cycles involved. These need to be clarified in a manner that a reader can understand the critical details of the model without having to dig into the model code. I will detail the points where I found something unclear in the minor comments below.

The second criticism is that the discussion of the limitations of the model is still a bit weak. One major point here is that it is not really discussed how much assumptions in the model parameterizations, such as the choice of stoichiometric N:P:Fe ratios influence the results, especially with respect to the question of the ultimate limiting nutrient. This will probably be done in detail in subsequent papers that use the model for that purpose, but some discussion here

would be in order. 

And thirdly, I think the manuscript should describe a bit the differences and similarities of their model to the carbon cycle component of the CLIMBER-X model (Willeit et al, 2023, doi: 10.5194/gmd-16-3501-2023), because it fills a rather similar niche in the ecology of carbon cycle models.

Minor comments

--------------

It is unclear to me whether the removal of a fixed proportion of nutrient uptake, mentioned in line 168ff is done uniformly over the ocean or locally. 

In line 229 it is mentioned that the nutrient framework is configurable; I understood this as that there is an option to add further possibly limiting nutrients. Is that so? It would be good to give a few more details on this.

In equation (2) and (3) why is the somewhat uncommon small delta is used instead of the more common partial derivative sign?

In equation (4) a nitrification term is mentioned; but the model, as far as I can see does not include ammonia as an explicit variable, only nitrate. On the other hand, if I look to equation (18) I see an explicit dependence on ammonia concentration. Can you explain? Does the model contain a prognostic ammonia variable that is also advected by the currents? I also find it slightly confusing that in Figure 2, middle panel, the heading says nitrate, but the arrows are annotated as DIN, which to me is the sum of nitrate/nitrite and ammonia.

In equation (7) a term proportional to the biological degradation of dissolved organic Fe is given. It is is unclear to me how to think about this organic Fe. Is it different from the ligand-bound Fe? How would one distinguish the two?

In equation (8) scavenging is made proportional to particle concentration, but is not explained whether this particle concentration varies with space and time, and what determines it.

In line 239 it is mentioned that phytoplankton nutrient uptake of Fe and P happens in a fixed proportion. As approximation this is ok, although the Fe:P ratio is much more variable in phytoplankton than the N:P ratio. But given the variability of Fe:P it may be a good idea to discuss a bit how this assumption may influence model results, especially concerning the question of the ultimate limiting nutrient. There are several recent papers discussing the effects of a variable C:Fe ratio, e.g. Wiseman et al. (2023, doi:10.1029/2022GB007491).

In line 259 ff it could be mentioned that an approximation of the vertical dependence with two exponentials is actually fairly similar to what one would obtain from the classical Martin (1987) curve, but is maybe somewhat more mechanistic. Just a suggestion.

In line 268 it is said that the consumption rate of electron acceptors in the remineralization is given in equations (12) to (14), but in fact these equations only give the inhibitory factors (dimensionless numbers between zero and one), not the rates themselves.

Also, it is mentioned that sulfate reduction is included here as a degradation process, but it is not made clear whether there is an explicit equation for sulfate concentration. Is sulfate simply made proportional to salinity or is it a prognostic state variable?

The section on ligand scheme enhancement (lines 300 ff) is somewhat confusingly written in several aspects. First, the classical one-ligand scheme with constant total ligand concentration is explained in some detail, without clearly marking that this is the old state of how iron is treated. Then in lines 318 ff the new scheme is explained, without clearly delineating it from the older approach. And then, it is described that the ligand stability constant is made dependent on depth. So is this the way that two ligands are described, by having one ligand, but with vertically varying stability constant? And if so, how is the vertical variation described, is it a step function, high at the surface, low below a certain depth, or something different? And what is the concentration of the ligand(s), is it/are they constant? Looking at the table 5 and carefully re-reading the text, my conjecture is that there is a constant ligand concentration

(but unclear what it's value is), and that the stability constant is approximately doubled in the uppermost surface layer. Is that it?

In line 376 ff it is explained that the modeled Fe distribution is compared to profiles obtained from GEOTRACES data, which is a feature that I really liked. But 'GEOTRACES data' is a bit vague; it should be noted which data set exactly is used (i.e. which intermediate data product) and the data should also be cited.

In line 383 ff it is stated that the model NPP is compared to a composite of several satellite-based estimates of NPP. But is it unclear to me how this comparison works: The model calculates net nutrient uptake in the uppermost model layer, and that can of course be converted to a net carbon uptake. But this, in my view at least, is more a calculation of the export production at the lower depth of the first layer at 88m, approximating the net vertical flux of carbon out of the euphotic zone, and not net primary production, which is significantly higher and also includes the carbon that is heterotrophically respired within the upper ocean. Or did I misunderstand something here?

Line 502: I first stumbled across the description; maybe just add that you are talking primarily about the deep phosphate concentrations here, not the surface.

Concerning the discrepancies to WOA surface fields mentioned for phosphate and nitrate (e.g. lines 515 ff): I think that some of the patterns look as if iron limitation is somewhat to weak in the major Fe-limited areas, like the Southern Ocean, the Equatorial and subpolar North Pacific. Maybe you can check this briefly? I would not be too surprised, especially in the Southern Ocean, and given the vertical resolution of the model. Overall the comparison looks quite good to me for a relatively simple model.

line 575: Maybe it should be mentioned that the fact that WOA does not contain iron data has to do with the limited amount of data and is not just an oversight by the makers of WOA data.

Figure 13: It is interesting that the nutGENIE model does a quite reasonable job in reproducing the patterns of N2 fixation, although it does not use a stronger temperature dependence for nitrogen fixation than for other phytoplankton, as many other models do. This is encouraging.

There are a few smaller typos in the references, e.g. in Ballantyne et al., and a few missing subscripts and capitalizations. Please check this once more.