This is an interesting study in which the authors use a suite of ESMs to evaluate future changes in NPP and their associated drivers in the EBUSs. The study highlights that the traditionally assumed key mechanism, upwelling-favorable winds, can only explain about 25% of NPP changes in the EBUS, while the remaining three-quarters are influenced by diverse drivers.

However, I do not recommend this manuscript for publication at this stage, as major revision is required. In my view, the current analysis of the relationship between upwelling-favorable winds and future NPP changes is somewhat simplistic and would benefit from a deeper investigation of model biases, inter-model differences, and regional variability. In addition, the discussion of the “diverse drivers” is uneven: analyses of physical processes such as geostrophic transport and wind-stress curl are extremely limited, whereas stratification and biogeochemical control of subsurface nutrient reservoirs are discussed in much greater detail. I therefore suggest that the authors thoroughly reconsider and improve the structure and balance of the manuscript.

Please find my detailed comments in the attached file.

General comments:

The paper by Cioffi et al. uses model projections from the CMIP6 ensemble to investigate future changes in net primary production (NPP) in eastern boundary upwelling systems (EBUS), and their driving mechanisms. This is a very important topic that has been mostly analyzed within individual EBUS until now, with studies across EBUS focusing on upwelling itself rather than primary production. The paper was interesting and well written and nicely demonstrates that the future of NPP within EBUS is highly uncertain.

However, I believe that the topic of future NPP changes deserves very careful evaluation due to the large uncertainties in biogeochemical (BGC) modeling (as pointed out by the authors). This is particularly problematic in climate models in EBUS that are known to be biased, notably due to a lack of model resolution (e.g., Varela et al., 2022, https://doi.org/10.3390/jmse10121970 or papers cited in the conclusion). As such, several aspects need to be improved before publishing the paper.

First, the “model evaluation” section (2.4) is wholly insufficient. The only evaluation presented compares mean NPP from all CMIP6 models to mean NPP from five data products. This is problematic on several levels: (1) each CMIP6 model should be evaluated separately instead of evaluating their mean (this will be important for the comment in the next paragraph); (2) the data NPP products should be used separately and not averaged; (3) other metrics beyond mean bias should be considered. Regarding (2), NPP satellite algorithms are known to display contrasting patterns and trends, and averaging products of likely varying quality is difficult to justify. I suggest the authors read the comment by Toby Westberry on the paper by Ryan-Keogh et al. (2023, https://doi.org/10.5194/essd-15-4829-2023) at https://essd.copernicus.org/articles/15/4829/2023/essd-15-4829-2023-discussion.html. Also note that the authors of that paper themselves analyzed the NPP products separately in their subsequent paper (Ryan-Keogh et al. 2025, https://doi.org/10.1038/s43247-025-02051-4 - I have no affiliation with either paper) with a discussion of this in the section “Assessing the merits of the different remote sensing algorithms”. Regarding (3), I would suggest looking at other metrics beyond the mean, in particular the 1998-2014 trend – before assessing the model trends in the future, it would be useful to evaluate their trend in the reference period. Other metrics could be used (eg standard deviation) if the authors think they are pertinent to evaluating how likely models are to be “right” in the future.

Second, I think there is a missed opportunity to narrow the projection uncertainty by considering the accuracy of each CMIP6 model in representing historical NPP (related to point (1) in the previous paragraph). I realize that most CMIP analyses only consider model ensembles, and this approach is valid for physical parameters that are well constrained by first principles equations. However, BGC model formulation can be highly variable and, even if BGC models similarly perform globally, their performance can be variable in specific regions (here EBUS) and for specific variables (here NPP). Moreover, the different models have an horizontal resolution varying from 0.25° to 1.5° in the 13-model subset, and resolution has a strong impact on upwelling representation and model accuracy in EBUS (see above). It would be very useful to define a “score” for each CMIP model that defines how well they perform during the historical period, and to see if the standard deviation of future projection decreases by only keeping the top X% performing models (half?). This is conceptually similar to the method used by Ryan-Keogh et al. (2025) although a simpler ranking scheme could be used based on model evaluation.

Third, some of the hypotheses in the discussion could be easily tested or at least explored further. For instance, if geostrophic transport and wind-stress curl are important contributor to NPP changes, there should still be agreement between w60, surfNO3 and NPP in regions where there is no agreement between tau, surfNO3 and NPP (section 4.2). Similarly, agreement between subsurface nitrate, surfNO3 and NPP could be tested (section 4.3).

Minor comments

Abstract “This uncertainty partially results from compensating positive and negative NPP anomalies within individual systems, with consistent multi-model responses only emerging at subsystem scales.” and similar statements elsewhere (l. 159-172): I don’t really understand what you mean. Are you referring to anomalies compensating between models (some models are positive, other negative) or regions (some subregions are positive, other negative)? In either case, that suggests higher variability in subregions than when considering entire EBUS, so I am unclear why this explains part of the uncertainty in projections. Also consistent multi-model responses are not very common according to Fig. 3 (1/3 of the coastal area perhaps? Should be quantified) so “consistent patterns projected across models” (caption) is contradicted by the figure, from what I see.

Statements l. 36 “Most past studies have focused on (…) wind-induced upwelling” and l. 75 “we reevaluate the standard perspective that perturbations to upwelling favorable winds (…)” are somewhat contradicted by l. 45-66 that highlight a rich literature supporting other drivers (stratification, geostrophic transport, subsurface nitrate etc). These statements should be rephrased.

l. 125-127: did you use the products by Ryan-Keogh et al. (2023), or calculated your own? In the first case this paper needs to be cited. In the second case, more details are needed since the authors only cite OC-CCI which, as far as I know (and found on the website) is only chlorophyll.

l. 280: do all CMIP6 models include iron limitation? If not it could be useful to see if model performance in regions where it may matter most (CanCS) depends on iron inclusion in BGC models.

l. 120: the 30m approximation has only been shown to be valid in the CalCS system; is it valid for other systems too? Also please note that this is defined from the mixed layer depth (as described in Jacox et al., 2018) rather than Ekman layer depth.

Technical comments

l. 105 what do you mean by “equatorward and downward component of surface wind stress”? Did you just mean equatorward?

l. 110, 207, 209 and elsewhere “nitrate” is a more common spelling than nitrates.

l. 227-229 it would be useful to refer to Fig. 5 right column (that supports this statement)

l. 354 cite the Biogeoscience paper rather than Biogeoscience Discussions