General comments:

The addressed question of this paper is very interesting and the applied statistical package to detect shifts in phenology seems to be adequate. However, I am not convinced that the strategy of using an area-based average gives meaningful results.  As mentioned by the authors, the North Atlantic shows a very heterogenous regional pattern between 30°-60°N.  Longhurst classified at least 4-5 distinct oceanic biogeographical provinces.  One might expect each biogeographic province to exhibit different temporal behaviour; e.g. Figure 5 clearly shows this for NorESM (a dipole structure for 2070-1850 and only an 1 day change of the mean peak NPP day over the entire domain, as stated in line 219). From my perspective, the analysis of peak NPP day averaged over the entire domain is therefore of little informational value.   

I recommend to repeat the analysis with the kernel based model for smaller domains (Longhurst provinces or regions aligned with common characteristics of the ESMs?). Then it might also be possible to identify local drivers of the change of peak NPP day. I assume that the intention of the investigation on the "Day of MLD <= 40m" was to identify one of these potential drivers. This MLD analysis came as quite a surprise and its findings should be mentioned in the abstract.  Again, I’m not convinced that the approach of using an average value of MLD day for the entire domain is reasonable. The relatively high cross correlation values with negative lag (in Fig. 9) might be related to this averaging.    

However, I appreciate this data analysis and therefore recommend publication of this manuscript with major revisions. Please see also my specific comments.

Specific comments:

L75: Please motivate the analysis of the change in MLD day already in the introduction.

L96: There is a difference between “primary production” and “net primary production”. The latter is the daily growth of phytoplankton minus respiratory demand. Please use net primary production throughout the manuscript.    

L101: Skyalls et al (2019) analysis did only cover a north-south transect east of 20°W from two cruises. Thus, the good agreement with observations is only shown for a very small part of the domain. Please add this information here.   

L119: Please add that phytoplankton growth is also a function of light and temperature in iHAMOCC.

L 124: Please replace the subtitle “Observations”. CAFE is a model based on satellite data, which strictly speaking are also only results of an algorithm (i.e. a model) and not observations.  Same for the subtitle 3.1

L163 : typo, capitalize Gulf

L171: typo, delete “the” in “For the the latter…”

L195:  Statements about the multidecadal variability of an 18-year time series should be avoided. They are meaningless.    

L214ff: It is difficult to reconcile the results of Fig. 4, which shows a shift towards an earlier peak NPP day in NorESM at the end of the simulation and Fig. 5, which presents an extended area with a later peak NPP day in the Gulf stream region.  In L219, you give a shift of only 1 day for NorESM between 2070 and 1850. From my point of view, the weakness of area-based averaging is clearly evident here. Could you also please provide the significance of the results of Fig. 5?  Is ±10 days significant, especially for EC-Earth?  It might be useful to also show the peak NPP day for the CAFE data set as a longitude-latitude plot. (see also comment L259)

L252: Please find different symbols for “l1” and “l2” for a better readability (e.g. capitalize L1 ?)

L255: Could you please give an explanation why L1 and L2 do not identify the 2070 change point in NorESM?  Would the results be more consistent if you decrease the penalty for L1 and L2 ?  

L259:  The finding of the kernel based model for EC-Earth is consistent with the findings of Fig. 5 (12 days change in peak NPP day compared to 11 ). However, the result for NorESM is quite different (10 days instead of 1). In addition, Fig.6 seems to give a much larger change for the end of the simulation than 1 or 10 days. Could you please elaborate this issue for NorESM?  Furthermore, why didn’t you already show the period 2010-2029 in Fig. 5?  It seems to contain a strong change pattern.  Why do you use 30-year averages for Fig. 5 and here you use a 50-year average?  

L273:  Please give this motivation to look at the MLD already in the introduction

 I have to admit that I’m not familiar with the used kernel based model. However, results shown for day of peak NPP in Fig.6 seem plausible: the decreasing penalty increases the numbers of change points that are found, but the change point with highest penalty is identical with one of the two changes points (there is an overlap of the pink line with one red lines in both models).  I became confused when I looked at the results of “Day of MLD <= 40 m”. There is no alignment between pink/red lines.  Why does the kernel based model miss the largest change point with a decreasing penalty?  

L292: Could you explain why both models have relatively high cross correlations for negative lag ( -3 to -1) years?

L313:  Here just to summarize that the change in peak NPP day is 1 day for NorESM is not convincing. For a different time period (1960-2010) versus (2010-2060) it was 10. It would be helpful to provide a little more discussion.  

L316 Please share your view on the results of the cross correlation analysis. Does it support the theory of Behrenfeld and Boss?  

Figure captions:

Fig. 1 Give period for the data averaging; is it 2003-2021 as in Fig 2 or 2002-2021 as stated in L135?  Did you also masked the model data as in Fig 2? Please indicate in the text if you always use masked ESM data for comparison with CAFE. How large is the difference between masked and unmasked ESM data ?

Fig. 3 The overlapping colour coding is difficult. Please find a better solution (e.g. instead of full time series show only mean/min/max for each ESM model.)  

Fig. 5 Please add significance level

Fig. 6 I assume that the centre of the symbols (circle and triangle) corresponds to the year of the change point. Please add this information. Please use the same names for I1 and l2 through the manuscript and in the figures. E.g. L1 and L2, not “model l1” or “model=l1”

Fig. 7  I assume you used the highest penalty level (corresponding to the pink line). Please note.  Please also change order of figures – usually EC-Earth is the top figure.

Fig. 8 The figure caption is incomplete ; add …in the median (model l1). Or just write : same as Figure 6, but for first day of ….

See attached review.