Review comments for the manuscript “Simulating vertical phytoplankton dynamics in a stratified ocean using a two-layered ecosystem model” by Zheng et al.

Reviewer: Camila Serra Pompei (DTU).

In this study, the authors developed a 2-layered NPZ model to better understand the mechanisms driving surface and sub-surface plankton dynamics in a stratified system. The BATS Ocean station is used as a case study, and the model is used to better understand the diverging phytoplankton trends observed in the surface vs subsurface layers in this location. The model provides an interesting avenue to better represent stratified systems with simple NPZ models, where the good model performance compared to field data highlights the strength of the approach. The article is relevant as it shows that different processes drive the dynamics of these two plankton communities, and that changes seen at the surface do not necessarily reflect changes observed in the sub-surface layer.

The manuscript is well written, clear, and model assumptions are well stated. My comments mostly center on the choice of euphotic zone depth and how this affects the representation of the sub-surface community and model comparison with data. There are also some aspects of the sensitivity analysis that need to be clarified.

Comments regarding the assumption of the euphotic zone depth:

The assumption of 0.001% might need to be discussed in the discussion section. Although phytoplankton might have been seen growing at these low light levels, is the choice of such a light level representative of where most phytoplankton in the subsurface layer reside? The choice of this low-light level results in the sub-surface layer always having a zeu deeper than 300m, and a depth range (zeu-zm) varying between 150 and 300m. On the other hand, when I look at BATS Chl data, the bulk Chl does not seem to go deeper than 180m, and anything deeper than that has an extremely low Chl concentration. So, it seems to me that the depth range of the sub-surface layer encompasses two layers in the “real world”: the deep Chl maximum and whatever is below it. So, in short, could the authors discuss how this assumption affects the model compared to the real world and how it affects its comparison to field data for the sub-surface layer?

The above comment specially applies in regard to the method chosen to compare the subsurface layer with field data: first field Chl is averaged between zm and zk and then multiplied by the thickness of the sub-surface layer in the model (zeu-zm). Can the authors mention if the sub-surface layer is representative of the zk-zm depth range? The value of zk is missing in the text and should be provided.

If possible, it could have been nice to have some additional figures (e.g. in an appendix) showing the vertical distribution of the field data used. This will help the reader have a better notion of the effect of the data transformation for model comparison.

Does the zeu assumption affect the light experienced in the sub-surface layer? How would Chlc,d and Chld change when choosing a different zeu value? Given the exponential decay of PAR and its averaging to use as input in the model, I expect Chlc,d not to change much, but Chld will certainly change (simply due to its multiplication with zeu). How would this affect the data comparison? In other words, it would be good to do a sensitivity analysis on the 0.001% value assumed to define zeu.

To sum up: I suggest that (i) an additional paragraph in the discussion section is added addressing the zeu choice and how does this compare to the observed depth range of the subsurface plankton community, (ii) a sensitivity analysis is performed on the 0.001% value, and (iii) some figures showing the vertical distribution of the Chl data are provided in an appendix.

Comments regarding the sensitivity analysis:

I am not sure I understand how z_e0 has been varied in the sensitivity analysis. Z_e0 is the initial value of the euphotic zone, but zeu is not a fixed parameter, it varies with surface light Kd, Chl etc. So, is this value only different in the first time-step and afterwards just follows the light forcing and kd and amount of Chl?

The same goes for z_m0, since zm is externally forced, does this mean that only the value in the first time step is changed and afterwards all is the same as the MLD forcing? Or does this initial condition somehow affect zm over the entire time-series? L188, why is the initial value of zm set as the mean of the ML? shouldn’t it be set as the initial value of the MLD provided by the data at the starting date of the model? Please provide an explanation on the choice of initial conditions for these forcing variables, why are they not the value provided by the forcing input variables on the corresponding day of the time series?

Other comments:

L 225 Just to clarify, is the “methodology outlined by Viljoen et al. (2024)” the one explained in the rest of this paragraph? Or is it something different?

L118  “but also involves the phytoplankton and zooplankton excursion from the surface layer (Eq. (14))” this makes it sound as if there was a vertical migration process. Perhaps rephrase this sentence as some kind of remineralization of dead plankton and sloppy feeding.

L124 mention in this part of the text that table 1 shows a description of the parameters.

L147, please indicate in this part of the text the units of N, P, Z and Chl in the model.

Equation 11, remove the multiplication signs (*) for consistency with the rest of equations. The (C:N) ratio could be given a parameter name (e.g. Q_{C:N}, or whatever you want to call it), since phytoplankton C:N ratio is fixed in the model. I also suggest changing “MWc” to contain only one capital letter and have the rest of letters as subscripts, as not to be mistaken as two different parameters.

Mention earlier in the text that the model has not been spinned-up but that it converges quickly. This is mentioned in the results section, but I was wondering about it while reading the methods.

L269, I would not use the word “histogram”, as the figure is not a histogram. Perhaps “bar-plot” or simply refer to the figure.

L283 typo in “notebaly”

L334 “modelling” should probably be “modelled”

L335 “is very similar to in observations” should probably be rephrased.

The authors developed a novel two-layered NPZ model to simulate the dynamics of surface and subsurface phytoplankton. Using BATS as the research object and adopting reasonable two-layered different parameter settings, they found opposite trends in surface and subsurface phytoplankton after 2011, and analyzed the reasons for this difference. The characteristics of phytoplankton are very consistent with observations in both seasonal changes and long-term trends, further enhancing the credibility of the research results. By establishing and validating this model, the authors provides support for studying the changing characteristics of ecosystems under current conditions. The design and results of this model may also have great applicability and can be further used to predict ecosystem changes under continued warming conditions in the future.

The manuscript is clearly structured and well written. The authors present a compelling study and vividly illustrate their hypotheses and conclusions using well-designed figures and clear explanations. I have some comments and suggestions that I hope will be helpful.

(Text directly from the article is enclosed in “quotation marks”.)

1. Questions

(1) The two layers of ecosystems are interconnected by dynamic transport of nutrients and light attenuation. The authors used a fixed mixing coefficient (μm) in this study. In equations (5, 17), there is a linear relationship between the nutrients interacting between the two layers and the depth of the mixed layer. But in the real ocean, this relationship is influenced by wind-driven, temperature gradients, and ocean circulation, and is not necessarily linear. I suggest that the authors further explain the effects of the dynamic transport of nutrients in the model on the change of phytoplankton in the two layers in the discussion section. This not only helps to understand the applicability of the two-layered model under real-world conditions, but also provides a reference for future improvements in the model to better simulate actual ocean processes.

(2) The time series of surface Ns observed in Figure 4(c) has many abrupt changes, while the model seems to follow some kind of seasonal cycle. Besides, in Figure 4(d), the standard deviation of the observed results appears to be much larger than the model, although the time-averaged results (Figure 4h) don’t show this (You also mentioned this in the article). Why the observation and models are so different on these? What are the potential consequences of this. (Is it related to the daily time resolution?)

(3) For the model results, the seasonal cycles in Figure 4(e) and Figure 4(g) looks highly correlated. In contrast, the seasonal cycles in Figure 4(f) and Figure 4(h) looks anti-correlated. While, for the subsurface, there seem to be other processes affecting their relationship in winter and spring (Jan-May). Can you explain these further?

(4) Line 397 “This correlation becomes strong over 2011–2022, reflected by an increase in Corr. to 0.75 (Table 4)”: The cause of the improved correlation isn't explored. Why does the correlation improve after 2011?

2. Suggestions for improving the content description

Abstract

(1) “reproducing trends post 2011 caused by ocean warming”: I suggest going straight to the specific main trends observed at the surface and in the subsurface.

(2) “simulating the ecosystem in the subsurface layer was more challenging than the ecosystem in the surface mixed-layer”: I would suggest briefly mentioning why simulating the subsurface is more challenging, or specifically stating which aspects (e.g. nutrient dynamics) are not accurate.

Introduction

In the fourth paragraph of the introduction section “A wide range of ecosystem models are available…”: This part seems a bit information-dense, which may make it difficult for readers to grasp the main points. I propose to categorize the models listed and describe the ecological problems that each type can solve.

Results

Page 18: For this paragraph “To understand the drivers of the decreasing trend in Chla_s over 2011–2022, we first show…”, the authors explain that the vertical-average result is also decreasing, indicating that the decreasing trend is not purely caused by the decrease in surface layer water volume. This explanation is very good. However, at the beginning of this paragraph, authors mentioned “we first show the interannual variability of observational mixed layer depth…”, readers would therefore expect this paragraph to describe the relationship between mixed layer depth and chlorophyll, and/or whether the trend of mixed layer depth affects the trend of chlorophyll. It is recommended to adjust/add appropriate content to make the content complete.

Discussion 

Line 515-517: It was mentioned that “The processes designed in this model do not incorporate all the key biogeochemical processes in stratified systems, such as nitrogen fixation and iron limitation”. Could you further explain how the lack of this part will affect this study?

3. Minor error

Informal expressions like “our ocean” and “less is known about it” should be replaced.

Figure 3 caption: “The integration of chlorophyll and nutrients from model and observation, and their relationship" can be revised to: The integration of chlorophyll and nutrients from model and observation, and their relationships.

Figure 4 caption: “(d) Daily nitrogen stocks at the surface layer”: This should be from subsurface.

Line 499: “Integration between 2011 to 2022” can be revised to: Integration between 2011 and 2022