Response and list of changes 

We appreciate the feedback and the many concrete suggestions for changes from the 3 reviewers. We have addressed most of them to the best of our ability. The major changes to the manuscript are: 

• Rewritten large parts of the introduction 

• Updated the “down-regulation” algorithm to make it simpler and ensure conservation of C, N, and Si under all situations. This required a re-calibration of the model, however, the results are largely similar 

• Fixed the issue with missing description of the transport matrices 

• Strengthened the justification for parameters for diatoms (new figure B1), validation (new figure 11), and sensitivity analysis (new figure F3). Further, the lower size limit was increased to become more realistic. 

• Added a new figure F8 to demonstrate that the bottom boundary condition does not control the results 

• Detailed responses are provided below. 

Specific responses to Reviewer #1 

The authors describe a library for their Nutrient-Unicellular-Multicellular (NUM) model, which can be used from Matlab.  The ms describes the model and compares predictions with observations and provides a few hints how it can be used.  The NUM model is advertised as a tool for ecological modellers without a biogeochemical or physical background. 

The model is built on principles and assumptions similar to previous plankton ecosystem models but does not appear to achieve a similar performance as similarly or less complex global ocean models representing the current state of the art.  While this is OK given the focus on simplicity, I see several problems both with the ms and the model itself.  To begin with, the authors vastly overstate the novelty and scientific rigor of the model or understate the current state of the art.  For example, this is neither the first nor the only model built on size-based parametrisation of multiple phytoplankton and zooplankton groups, the same holds for mixotrophy, and many of the underlying assumptions are not based on first principles, contrary to what is written in the ms.  Related to this is a lack of consideration of and comparison with the pertinent literature.  Further, essential information about the model is either missing or difficult to access, e.g., how the transport matrix was extracted and from which model, or how the model was calibrated.  Given the many deficiencies of the model (basically any current model) I find it inappropriate to advertise the use by modellers without biogeochemical or physical background, as this background is essential for a meaningful application.  Finally, the model is strongly nudged towards bottom concentrations of N and Si.  As it stands, this not only makes the (essential) analysis of mass balance extremely difficult, which is further complicated by fluxes to unresolved compartments (higher trophic levels), but also put into the question the interpretation of and conclusions drawn from the model results.  These issues, in particular the ability to check the mass balance, must resolved to make this ms publishable. 

We appreciate the rigorous feedback and the detection of important issues that has helped us improve the transparency of the model documentation, its presentation, framing and discussion.  

It is not our intention to claim that the model itself is novel. It is part of a class of size-based models, and we find that these models are well referenced in the manuscript. The novelty is the presentation of a self-contained implementation of the model. 

The issues about the missing transport matrix were due to an unfortunate error where a paragraph was omitted from the submitted manuscript. That is fixed now 

We disagree about the importance of the nudging. In this regard, we have added a new figure (F8) to demonstrate that results are indeed very insensitive to the value of the bottom boundary condition.  

We do agree that there were unclear issues about mass conservation in figure 4 and have specified that the equations indeed conserve mass of all nutrients and carbon. 

Another, though comparatively minor, problem is the choice of Matlab as the high-level interface.  This restricts the use to those who have access to a Matlab license.  More importantly, however, Matlab has repeatedly introduced backwards-incompatible changes, making older programs difficult or impossible to use with recent Matlab versions, e.g., the ROMS-AGRIF toolbox.  This is difficult to deal with also because Matlab is not open source.  Thus, for future developments, I would suggest to switch to open-source software, e.g., Python or R, both of which have advanced numerical integration packages and can be used with Fortran. 

We agree that the choice of matlab is not ideal. We would indeed like to switch to python.



Specific evaluation 

In the Abstract, l4: "The distinguishing feature is that body size is used for structuring predator-prey interactions and to scale all parameters".  This makes it sound as if this was new or specific to the presented model, but the same concept has been employed previously by Ward and Follows (2016) and Ward et al. (2012: L&O 57), and is was also outlined in Andersen et al. (2015: JPR 37). 

We have changed “distinguishing” to “central” to avoid the misinterpretation that basing models on size is novel. In the introduction where we refer to size, we indeed refer to Ward and Follows and several other models (lines 31). 

Several statements in the Introduction regarding the current state of the art are misleading at the very least: 

l20: "focus on phytoplankton" This is not quite correct, ESMs generally focus on export and most of them place similar focus on phytoplankton and zooplankton, e.g., PISCESv2 (Aumont et al, 2015, GMD 8) has 2 phyto and 2 zoo groups and COBALTv2 (Dunne et al., 2020, JAMES 12) has 3 each.  In fact, I am not aware of a current Earth system model merely treating zooplankton as a closure term (l22).  Some global modelling studies even focus specifically on zooplankton, for example Prowe et al., (2012: Progr. Oceanogr. 101; 2019: L&O 64) and Ward et al., (2012).  In particular, Ward et al., (2012) and Prowe et al., (2019) appear also highly relevant for the present manuscript and should be discussed/contrasted with the present model, also Ward and Follows (2016).  There are likely more.  Without relating the present to previous models is will be very difficult to assess the usefulness of the present study. 

l21: "fast turnover rates mean that they have the highest loss rates to the biological carbon pump": This is not true, the fast turnover rates characterise quick recycling in the microbial loop. Only export enters the biological pump and this is largely channelled through zooplankton. 

l32: "most eukaryotic plankton are neither phytoplankton or zooplankton but mixotrophic jacks-of-all-trades": This is not correct. While many species are indeed mixotrophs, most members of the groups of major importance to the biological pump and higher trophic levels, i.e., diatoms, coccolithophorids, copepods, are indeed predominantly autotrophs or heterotrophs. Mixotrophy is best established for dinoflagellates and photosynthetic bacteria and even there one trophic mode is usually dominant and the other restricted to some auxiliary processes, meaning that, with few exceptions, the classification as functional autotrophs and heterotrophs is still valid and useful.  Also, the ability to switch between autotrophy and heterotrophy seems to be a rare exception.  This is a problem in models representing mixotrophy, including the one under consideration here, since they assume that this switching is a regular feature of mixotrophy. 

l36: "sufficient representation of the functions of marine planktonic food webs ... requires a representation of the mixotrophic nature of eukaryotic plankton":  This statement, while I have seen it quite often already in several forms, has never been substantiated.  Just that a phenomenon exists and is not represented in models does not mean that it is important, even though it may have strong effects in the model (and several modelling studies have shown that this is the case).  Either remove this statement or at least tone in down. 

Response to the four comments above. We agree that this part of the introduction does not work well. It tried to frame the NUM model as being a solution to problems in other models. And some of the alleged problems were, as the reviewer correctly points out, not actually problems. We have rewritten this part of the introduction completely to explain the NUM model’s place in the landscape of plankton models. We believe it is more balanced now. 

The importance of mixotrophy for the ecological efficiency was a key result of Ward and Follows (2016), which we had repeatedly referenced. Nevertheless, this sentence is removed in the new introduction.  

The NUM library and model description (Sections 2 and 3) is very difficult to follow, at least in part because many model equations are collected, without equation numbers, in tables in Appendix B.  Several statements are either not correct, unclear, or misleading. 

We find that model description would be harder to read if the equations are included in the text. We have therefore maintained the current structure with only the important equations in the main text and the full set of equations in tables in the appendices. 

l78: "All the setups are derived from a global transport matrix": Where does this transport matrix (TM) come from?  How was it extracted and from which model? TM extraction is not a trivial task and only very few researchers can do this currently.  Thus, some details are required here (or at least somewhere in the model description).  Also missing is the temporal and spatial resolution of the TM. 

A key part of the section about the transport matrix was unfortunately omitted from the submitted manuscript. It is now included (section 3.7) and here we describe the transport matrix, including its spatial and temporal resolution. 

l82: "DOC, nitrogen, and silicate": Unclear why these have been chosen. What about P and Fe? 

We have added a justification: “The 3 nutrient variables are the minimal set needed to represent growth limitation for phytoplankton (incl. diatoms which are limited by silicate) and the carbon source for osmoheterotrophs (bacteria).” 

l109: "The vacuole requires a silica shell":  No, it only requires a cell wall.  The vacuole is not unique to diatoms but fairly common in other groups as well, e.g., green algae. 

Rephrased as “Further, diatoms have a silicate shell, which makes silicate another limiting nutrient”  

l107: "This trait allows them to increase their volume per carbon without additional nutrient demands, which in turn leads to lowered predation pressure": This statement, while frequently encountered in the literature, lacks observational support and should be toned down.  In fact, diatoms seem to be the preferred food for many copepod species, and their success appears more related to their high growth rates rather than predation avoidance. 

Is it the increase in volume or the lowered predation pressure due to size or the silicate shell which needs observational support? Regarding the first, we have written that: “in the model, this trait allows then to increase their volume”. Regarding the second, we have erased the mention of the lowered predation pressure due to the increased size. We have instead added a more recent reference with direct support for reduced copepod grazing due to the silicate shell: “The hard silica shell also reduces the predation pressure on the diatoms \citep{hamm2003architecture, pancic_silicified_2019}”. 

l119: "spring blooms are controlled by a time lag between the maximum consumption by adult copepods and the rise of juveniles (Vidal and Smith, 1986)": This is not what Vidal and Smith wrote.   

What they did write is that the development of one herbivores (copepods) community was affected by this time lag, not that of the (phytoplankton) spring bloom.   

We agree; this is an incorrect citation. Vidal and Smith demonstrated how copepod growth was determined by the spring blooms and not the other way around. We have removed the sentence. 

Interestingly for the present work, they also wrote that "growth rates of planktonic crustaceans must not be scaled just on the basis of body size", which is what current size-based models do, including the present.  I think discussing this kind of caveats of allometric scaling is essential for the advancement of size-based models.  Sweeping these problems under the carpet is only going to aggravate problems in the future. 

In our interpretation, Vidal and Smith refer to development (moulting) rates when they talk about “growth”. First of all, growth (increase in mass) in the NUM model depends on the food available. Body size only enters through the maximum growth rates. Growth itself is therefore not determined alone on body size but crucially about food conditions and temperatrure. How maximum growth rates scale with mass across species is well demonstrated in the literature (see e.g. Serra-Pompei et al 2020 fig. 2a). We disagree that we sweep anything under the carpet. In fact, we resolve some of the additional non-size related diversity with the feeding-type trait of the copepods. 

l125: "representative populations": Representative of what? 

We have removed “representative”, such that the sentence reads: “Modelling the diversity of copepods with these two main traits makes it possible to represent the entire copepod community in the global ocean with just a few populations that spans the trait space of adult size and feeding modes.” 

l129: "we have used a single POM group": Unclear: Does this mean only one size class? 

Yes. Clarified now: “POM may be represented by a size spectrum to enable the calculation of sinking velocities for different POM sizes. However, in the simulations presented here we have just a single size class of POM” 

l160: "comply": does not make sense, use a better word, e.g., allow for 

We have changed to “respect the increased prey vulnerability”. 

l167: "they [diatoms] have a lowered predation risk": They are the preferred prey of many copepods (see above). 

See response above about the reduced predation risk due to the silicate shell. 

l196: "To avoid unnecessary respiration, the cell regulates the actual uptakes jX from the encountered uptakes jenc.X to maximize cell division rate under the constraint of a fixed C:N:Si stoichiometry": According to Table B1 this seems to be a simple law of the minimum, dictated by the fixed-stoichiometry assumption and should be described as such.  The formulation "to maximize cell division rate" is misleading. 

It is not just a simple law of the minimum. The cells may have different sources of carbon (light, DOC, food) and different sources of N (N, food). They choose the combination of sources than minimizes the metabolic costs of uptake. This is essentially a fitness optimization. We have explained this in the text now (l226-232) and included a reference to the specific equations in table B2. 

l247 and Eq. (19): What is the vertical resolution?  This equation will exert a very strong relaxation on the whole ocean grid, including the surface layer, which is usually much thinner than the bottom layer, and there the gradient with respect to the bottom concentration will be strongest also. 

An entire paragraph leading up to this equation was missing from the submitted version. We apologize for this. We have now included the missing text (l286-297). 

Regarding the bottom boundary condition, we do not fully see the basis of the reviewer’s concern. Surface nutrient concentrations in our model are set by plankton uptake and by the mixing rates defined in the transport matrix. Given this structure, we are unsure which alternative boundary condition the reviewer considers preferable. One alternative could be a no-gradient boundary condition, however, that would leave the bottom boundary concentrations completely free. Another alternative would be to remineralize all nutrients from POM at the bottom. While this would ensure conservation, it would also require excessively long run times to equilibrate the nutrient concentrations. Representing a complete nutrient cycle  is not the aim of our model and it  is not required for the more ecological questions that the model is designed to address.  

To address the reviewer’s point explicitly, we performed a sensitivity analysis of the bottom boundary concentration. The model shows very low sensitivity to this parameter. We have added a new figure (F7) that demonstrates that varying the bottom boundary concentration produces negligible changes in the model results. 

l282: "All parameters ... are based on first principles or cross-species comparisons": This is simply not true.  For example, representing temperature dependence by one or more Q10s is not based on first principles and also often does not hold among species, and the same applies to prey preference, viral lysis, etc.  In fact, very little of the theory (assumptions) behind the present model has a sound empirical foundation, which should be obvious from the above and below comments. 

We do not share the reviewer’s assessment that the model lacks a solid empirical or theoretical foundation. To fortify our description, we have edited the section regarding temperature to specify the theory our formulation is based on. We have also rewritten the relevant paragraph on parametrization (lines 67–73) to describe more clearly the theoretical and empirical bases for the remaining parameters and assumptions, and to briefly elaborate on our approach based on first principles. 

The temperature scaling implemented in our model applies temperature corrections directly on the processes. This is the topic of the paper Serra-Pompei et al (2019) that is referenced in the section about temperature dependence and further analyzed in Andersen and Visser (2023). Temperature dependence of growth is therefore an emergent property, and it does not follow Q10 relations. We have now emphasized this is the section dealing with temperature corrections (l282-284). Indeed, the temperature correction with Q10s does not hold among species, yet our model does not represent individual species, and  the correction is applied on the metabolic processes of the organisms.   

l289: "Lower higher": Unclear which one, lower or higher? 

Rephrased to “Decreasing the mortality by higher trophic levels increases the copepod biomass…” 

l291: "We calibrate by minimizing the error between average modeled and observed values": How was this done? 

We have rewritten this part, by presenting the steps in the calibration algorithm in a way that is more clear and easier to follow.  

Did you account for possible correlations between the three calibrated parameters or among the variables used in the cost function? 

No. We have added a reference to the figure that shows how the calibration is affected by each parameter (Figure F6). 

Eq. (21): What are Q and P? 

We have specified that in the text now. 

l296: "three sites that represent markedly different environments": Which sites are these?  They are never identified. 

The locations are specified below in lines 360-362. 

l310: "innocuous": makes no sense.  Did you mean simple? 

We have changed to “appears simple” 

caption Fig. 4: "all fluxes into POM are directly mineralized into dissolved silicate": The arrow through (1 – γ2) in the right panel points to POM.  Should it not point to Si? 

Correct, this is a mistake (and a few other) in the figure which has been corrected. 

  How are the fluxes to higher trophic levels accounted for?  For example, what happens to the N going into HTL? 

It is written in the equations, but not on the figure, which was a mistake. We have now added the missing arrows. 

This, together with the nudging (Eq. 19), raises questions regarding mass balance.  Checking mass balance is essential for asserting the internal consistency and integrity of numerical models.  Models employing nudging or fluxes into unresolved compartments thus need to provide special tools for accounting for the fluxes implied by the nudging and outside of the model domain.  Without these, the NUM library will be completely useless. 

We have added: “We have checked conservation of all nutrients in Eq. 12-17 and they are conserved within a factor close to the numerical noise.”. In addition, the model code includes a function named “checkConservation” that performs the calculation to assert the internal consistency and integrity. As written above, we have added a figure to demonstrate the robustness of the bottom boundary condition (Fig. F7).  

caption Fig. 5: "Sheldon size distributions": It should be explained briefly what this is. 

We have added a brief explanation in the caption noting that the Sheldon size distribution refers to “Biomass normalized by the width of each size bin” and included a reference. 

Results: 

l325: "with a weak effect on the biomass distributions": This is difficult to believe, owing to the very strong effect on surface nutrients according to Eq. (19). 

We included figure F7 to demonstrate that the effect is indeed weak. 

l330: "(3.8)": What does this number refer to? 

We have specified “(Section 3.8).” 

l375: "diatoms can also take up dissolved DOC and therefore a fraction of their biomass will be categorized as bacteria": This is very confusing.  On l365 it is stated and Fig. 4 (right) indicates that diatoms are obligate autotrophs, i.e., do not use DOC.  Table B1 also does not specify which parameters apply to which groups, so this should be clarified and/or the statements and Fig. 4 corrected regarding diatoms. 

This phrase is indeed confusing. We have specified that “[diatoms] can therefore not perform phagotrophy and they obtain nutrients only from diffusive uptakes of dissolved matter”. We have also specified that when we say “obligate phytoplankton” DOC uptake is also included. 

We have updated the figure 4 caption to state the “U” in the figure represents all unicellular plankton: generalists + diatoms. 

The equations in Table B1 applies to both generalists and diatoms, as stated in the caption. There is one exception which is flagged with a footnote.  

l379: "The copepod community is dominated by active copepods in higher latitudes (Fig. 13b) and by passive (ambushing) copepods in oligotrophic gyres":  This is the opposite of what has been reported by Prowe et al. (2019).  What differences between the two models could be responsible for this difference? 

First, we have removed the part of the text here that was more suited to the discussion.  

In the discussion paragraph on copepods, we have now added Prowe et al (2019) to the list of data analyses of copepod feeding traits. There is a clear difference in the observed patterns between Prowe et al and Benedetti et al. We believe that this is due to fewer data in Prowe et al and have added: “The pattern of active copepods in high latitudes and passive copepods in low latitudes is different that observed and modelled by \citet{prowe_biogeography_2019}.  We believe that the pattern in passive/active copepods were due to a very limited set of observations, in particular at low latitudes.” (l564-566). We have no explanation for why the two models give different results. 

l382: "In the Equatorial Pacific though, passive copepods dominate, possibly because there small primary producers are more common ..., whose size is preferred by these copepods": This does not appear to make sense.  Ambush feeding is effective only for motile prey, which is why ambush-feeding copepods are usually relatively large, carnivorous species with a preference for large prey. 

We have deleted this section of the text. The observations were also redundant with the section above which discusses the high active-copepod-ratio – when this ratio is high the passive ratio is low. 

l407: "In less seasonal environments, the absence of advection appears to result in too low production": Vertical mixing in oligotrophic regions is often driven by internal waves.  Are these included in the 1D simulations? 

The water-column simulations are using the same transport matrices that are used to drive the global model. The only difference is the absence of lateral advection. We have now specified that it is “lateral advection” which was missing. 

l421: "generalists combine phototrophy (yellow lines) for carbon uptake with predation (red lines) for nitrogen": Why do they not just use the carbon contained in their prey? 

We have specified: “… generalists combine phototrophy (yellow lines) for carbon uptake for respiration with predation (red lines) for nitrogen and carbon for biosynthesis (and possibly respiration)”. 

Fig. 14: Please add labels (global, water column, chemostat) to the panels. 

Done. 

Fig. 15: axis label missing for the x-axes 

Done. 

Discussion: 

l434: "Overall the NUM model is an intermediate-complexity model laboratory for ecological simulations that is accessible to marine ecologists without biogeochemical or physical oceanography background": Given the many deficiencies of the present model and the general lack of consensus about the formulation of the plankton-ecosystem components of ocean models, I find it inappropriate to advertise the use by people without the scientific training required for a meaningful employment of this kind of tool. 

We agree that ecological models require informed use, yet we maintain that the NUM model is appropriately described as accessible to marine ecologists. Its design focuses on process-based ecological interactions structured by body size, with formulations for growth, feeding, and trophic transfer that do not require specialised training in physical oceanography or full biogeochemical modelling. In revising the manuscript we have clarified several parts of the documentation, including parameter choices and mass conservation (see comments on Fig. 4 and Eq. 19), which improves transparency for users. Regarding accessibility, we acknowledge that Matlab is not ideal and we plan to transition the code to Python. The model captures key large-scale biomass, production and nutrient patterns shown in the validation section, which supports its usefulness for ecological questions. Its intermediate complexity, transparent equations and limited parameter set allow ecologists to explore hypotheses in a controlled framework while remaining aware of the model’s assumptions and limits.We also note that Referee #2 explicitly highlights the model’s potential usefulness for plankton ecology. Having addressed the main concerns raised, we have retained the original sentence 

l449: "a deeper production maximum": This is most likely a model artefact (these have never been observed), resulting from the assumption of a fixed stoichiometry and Chl:C ratio.  Models with variable Chl:C ratio generate deep chlorophyll maxima but these are not associated with production maxima. 

This statement does not refer to observations, but to what occurs in the three versions of the model. We have adjusted the text to specify “deeper plankton maximum”. Nevertheless, we disagree that deep Chl maxima have never been observed in seasonal environments. We have added a reference (Yasunaka, Sayaka, et al. "Global distribution and variability of subsurface chlorophyll concentration." Ocean Science Discussions 2021 (2021): 1-22). 

l449: "the global simulation reproduces the deep maximum less accurately": In order to assess the accuracy of a model feature, it must be compared with observations, not other model configurations, but this is not done here. 

Rephrased to “However, the global simulations reproduction of the deep maximum is coarser than in the water column” 

l455: "This finding points to the importance of a better understanding of plankton’s ability to maintain a position in the photic zone despite vertical mixing": I don't think so.  Mixing losses apply identically to other unicellular species as well.  Much more likely it is a misrepresentation of the microbial loop. 

This is a statement about comparing water-column and global simulations with the chemostat. We have rephrased to clarify this. 

l463: "The smallest diatoms are predominantly osmoheterotrophs (bacteria)": No, diatoms (of any size) are not bacteria. 

Indeed, they are not bacteria. We have increased the lower size limit of diatoms in the updated model to eliminate the unrealistically small diatoms, so this comment has been removed. 

l469: "We conjecture that the absence of large diatoms is due to mixing losses": See above comment for l455 

See our response to comment for l455. 

l477: "the model produces global patterns of biomass, production, and surface nutrients reasonably": This statement is in stark contrast to what is shown in Figs. 8, 10, 11. 

We have borrowed the statement from reviewer #2 to rephrase that sentence to: “Despite this simplification, the model reproduces important global-scale patterns of plankton biomass, nutrient fields, and net primary production with some accuracy.”. 

l514: "what mechanism provides diatoms a faster cell division rate than other phototrophic plankton": This has been linked to the energetic efficiency of a siliceous cell wall (Martin-Jézéquel et al., 2000: J. Phycol. 36) 

Thank you for this reference. We have removed the question. 

Appendix B: 

Please add equation numbers to the equations in the tables. 

Done 

Appendix F: 

Fig. F3: Please add labels in the panels and units to the colour bars.  It could also be helpful to use a 2-column layout.  Also, the numbers on the colour bars are nonsensical, indicating negative concentrations.  There is only one white area (in panel b) and it is unclear why it is white and not green (as 0.01 should be according to the colour bar). 

We have put headings on each panel, fixed the white area, and changed the color bar labeling to a logarithmic scale.

We appreciate the feedback and the many concrete suggestions for changes from the reviewers. We have addressed most of them to the best of our ability.

Specific responses

The submitted manuscript presents a significant development of the Nutrient-Unicellular-Multicellular (NUM) plankton modelling framework. This includes the extension of the physical framework to global, water column and modelling environments, and the addition of a diatom functional group. A unique feature of this global model is the emergence of different trophic strategies as a function of organism traits and environmental conditions. 

The model is calibrated against multiple datasets, but there is some ambiguity regarding the spatial resoloution of the assimilated data (see comments). Nonetheless, the model reproduces important global-scale patterns of plankton biomass, nutrient fields and net primary production with some accuracy. 

I found the paper to be well written, subject to some corrections, as suggested in an annotated copy of the submitted pdf. Of these suggestions, I have highlighted a few of the more important points below. 

• Please could the authors say a little more about what they mean by “first principles"? I am not convinced that all the phytoplankton parameters are derived purely from physical constraints, but it is possible the authors mean something different. Regardless off this, the authors later step back from this ambition slightly, stating on Line 282 that "All parameters of the plankton groups are based on first principles or cross-species comparisons". The latter is more realistic and this qualifier should be included throughout.  

We have now elaborated on this: “First principles are those relations that can be related directly down to physical laws (diffusion, fluid mechanics, etc.), energetics of chemical reactions, geometry, or mass conservation.  This ambition is realized for many of the parameters for the unicellular compartment \citep{andersen_cell_2023}.  It more difficult for the multicellular compartment, and here most parameters are determined by cross-species analyses.  The description of the parameter values are done elsewhere \citep{serra-pompei_general_2020, cadier_competition_2020, hansen_seasonal_2019, andersen_cell_2023}“ 

• Section 3.4 suggests that there is no non-advective export of Si included in the model (i.e. particulate Si does not sink. This seems very unrealistic, and potentially an important omission. I did not see any indication that living diatoms had a sinking speed. Please could the authors clarify and justify this choice. 

It is certainly correct that the handling of silicate POM is not ideal. We have now considered that silicate which should have been in POM is quickly sinking and therefore lost. We have added more text in the discussion about the diatoms to highlight this omission: “Further, silicate is not incorporated into the POM, which means that it is lost from the model.  In light of this very simple description of diatoms, and the important omission of non-sinking silicate POM,…”  and “Finally, silicate should be included into the POM group. This would require another state variable to account for the variable stoichiometry of POM.” 

• When discussing or plotting size spectra, could the authors please use the length dimension instead of (or as well as) mass. It is much, much easier for readers to understand how big the organisms are in terms of diameter. I have no idea how big a 0.01 µg C plankton is. 

We have chosen mass as the dominant axis for plots because it reflects the most important attribute of the organisms. Plotting with a length-axis in not just a change of scale, because the relation between length and mass differs between generalists, diatoms, and copepods. To assist the readers who prefer lengths, we have now added three length scales on the first figure with size spectra (Fig. 5). We further mention length in parentheses when we mention masses in the text. 

• After reading the ms I was unsure if these first three data types shown in equation 21 are spatially resolved or a single global mean. Much of the text suggests the latter, but if this is not clear from the description in section 3.9. This requires clarification. 

We have rewritten the description of this equation to clarify that the comparison is made at each observation time and space. 

• I expect that the low N and Si concentrations in the ocean surface are better explained by unrealistically efficient nutrient uptake by phytoplankton, rather than underestimated external supply. For surface Si, the inclusion of unrealistically small diatoms could certainly contribute to this. 

We have now excluded the unrealistically small diatoms by changing the minimum size of diatoms. We also agree that the low concentrations are due to the efficient uptake of Si (and N) and we have rewritten this paragraph to center on this explanation. 

• Please use micromolar units for nutrient concentrations (and possibly biomass), which are much more standard in ocean biogeochemistry. 

We have now done it for nutrients. 

We appreciate the feedback and the many concrete suggestions for changes from the reviewers. We have addressed most of them to the best of our ability. The major changes to the manuscript are: 

• Rewritten large parts of the introduction 

• Updated the “down-regulation” algorithm to make it simpler and ensure conservation of C, N, and Si under all situations. This required a re-calibration of the model, however, the results are largely similar 

• Fixed the issue with missing description of the transport matrices 

• Strengthened the justification for parameters for diatoms (new figure B1), validation (new figure 11), and sensitivity analysis (new figure F3). Further, the lower size limit was increased to become more realistic. 

• Added a new figure F8 to demonstrate that the bottom boundary condition does not control the results 

• Detailed responses are provided below in. 

Specific responses to Referee#2

The manuscript, “Computational library for the Nutrient-Unicellular-Multicellular plankton modeling framework v. 1.0”, by Amalia Papapostolou et al. presents a computational modeling framework for the recently published NUM (Nutrient-Unicellular-Multicellular) plankton model and example results from its application to lower-trophic ocean ecosystems at sites typical of many relevant oceanic observations of plankton ecosystems and biogeochemistry. This manuscript seems very well suited to the scope of Geoscientific Model Development (GMD) and promises to be of interest to a substantial portion of this journal’s readers. 

As that authors state, this manuscript does present what appears to be a useful computational framework, which can make the NUM model accessible to ecologists who are not model developers / computational oceanographers. As such, it has the potential foster new studies of the plankton ecology and biogeochemical cycles together with the ocean circulation and mixing that drive important aspects of their dynamics. However, some of the writing is vague and incomplete (i.e., leaving out important facts very relevant to the ecological dynamics, such as the limitation of diatom growth/cell division by the availability of dissolved silica). Also, some of the figures present a great deal of information, which is challenging to interpret with the minimal labeling provided. More labels for the panels (or rows/columns thereof) would make it easier to understand what is presented. 

Overall, the manuscript is reasonably well written. Please see my specific concerns below for some suggestions to improve the clarity and accuracy of the writing and figures. My assessment is that after moderate revisions, this manuscript should be suitable for publication in GMD.   

Specific Concerns: 

Line 59: This is really a technical, not scientific challenge. 

Correct; fixed. 

Line 74: The model equations are ordinary differential equations (ODEs) only for the chemostat (0D or ‘box’) model configuration. When solved as spatially explicit equations (1D water column model or 3D ocean model), they are Partial Differential Equations (PDEs), which include spatial as well as temporal derivatives. 

This is not the case. Since we use transport matrices to resolve space, solving the plankton model is really “just” about solving ODEs. The transport matrices represent a discretized version of the PDEs of the global flow fields. We have now clarified that in the text: “All the setups are derived from a global transport matrix that represent a discretized version of the underlying geophysical partial differential transport equations.” 

Section numbering seems to be off for Section 2:  Should it not be Section 2.1 rather than Section 2.01, etc? 

Correct; fixed. 

Section 2.02 Diatoms: It should be made clear here that the model does account for the limitation of diatoms’ growth rate the availability of dissolved Si (silicic acid) and that the reproduction (cell division of diatoms is well known to be limited by the availability of dissovled Si (silicic acid). 

Done. We have written: “The silicate shell also means that diatom growth is limited by dissolved Si, besides light and nutrients.” 

Line 67: delete the first “and”. 

Fixed. 

Also, should this not be Section 2.2, rather than 2.02? 

Correct; fixed. 

Section 2.0.3 Copepods: It would be good to mention the substantial inter-specific differences in size for copepods (e.g., very large in nutrient-rich subarctic waters versus very small in nutrient-poor subtropical waters), and how the NUMI model could capture this. 

We have added a note and a reference: “In this way, the model represents the large differences in adult size observed in the oceans \citep{brun_trait_2016}.”. 

Section 3.2 Cell model: Mention imitation of growth by dissolved Si (silicic acid) for diatoms, which is an ecologically important difference between diatoms and other unicellular plankton. Therefore, the vacuole is not the only difference between diatoms and other unicellular plankton in this model. 

We have added in the section 2.2 Diatom “Further, diatoms have a silicate shell, which makes silicate another limiting nutrient.” 

Line 183, eq. 7: According to this equation, the fraction of active biomass decreases with r and INCREASES with delta, which is the opposite of what one would expect. Is this not the INACTIVE fraction of cell volume? This also appears much simpler than one would expect based on the equations (in Appendix B) for the membrane volume. If this equation is correct, please provide an explanation in the text, or otherwise correct it. 

The equation is incorrect. Yes, the active fraction is indeed 1-3*delta/r. This is fixed now. 

Line 202, eq. 10: Jnet is not defined here, nor in Table B1. Therefore, it was not clear, until looking carefully at Table B1, whether growth rate was limited by nutrient uptake. This should be made more clear. 

We have now defined the equations B1.14 and B1.15 and made it explicit that division rate is indeed dependent on light and resources. 

Lines 310-317: It is worth mentioning here that in subtropical sites, DOC (adsporbed to the filters) can contribute a variable but substantial portion (up to about half) of PP measured by typical incubation experiments (Karl et al. 1998), and that studies focusing on subtropical sites may benefit from accounting for the contribution of DOC to PP. 

That’s interesting, we have not been aware of that. We have added “The DOC part of primary production can contribute a substantial portion of measured primary production \citep{Karl1998role}.” … “However, we refrain from including DOC NPP, which means that we may underestimate measured NPP at oligotrophic sites.” 

Karl, D. M., D. V. Hebel, K. Bjorkman and R. M. Letelier (1998): The role of dissolved organic matter release in the productivity of the oligotrophic North Pacific Ocean. Limnol. Oceanogr., 43(6), 1270–1286. 

Figure 9: A legend defining the lines and symbols would help to make the figure easier to understand. 

We have added descriptions on panels b-c: “spring/summer”, “AMT fall”, “Spring” on panels b-d 

Line 391: “an” should be “and”. 

Fixed 

Figure 14: This rather complicated figure was difficult to interpret at first. Labels indicating which panels / lines are from the chemostat (0D), water column (1D) and global (3D) models would help to make this figure much easier to understand. 

We have included such labels on the right side of the figure now. 

Line 441: “underway” should be either “under development” or if already developed and published, “available” or similar. 

Fixed. 

Line 465: It is not clear what “Including this effect” means. Re-word to specify the meaning: 

For example, “Imposing a realistic lower size limit for diatoms would…” 

In the updated simulations we have imposed a realistic lower limit. This sentence has therefore been removed.