Reviewer’s comments on the manuscript “Anthropogenic climate change drives non-stationary phytoplankton variance” by Elsworth et al. (2022)

General comments:

The authors investigated future changes in interannual variability of phytoplankton carbon biomass by using the CESM1 large ensemble simulation. Their results seem to indicate highly spatially heterogenous response of interannual variability in the biomass to the global warming by the end of 21^st century and relatively important contribution from changes in “top-down” control of the phytoplankton growth.

I totally agreed with the authors’ initial point that, especially in the context of ocean biogeochemistry, future changes in variabilities have not been paid much attentions compared to those in the climate mean states, although these are critically important on decisions of mitigation and adaptation policy. I don’t think that this study has no potential for being a step to help our understanding about the ocean ecosystem (including from lower to higher trophic levels) response to climate changes. However, I can not recommend the editor to publish the current manuscript from the BioGeosciences, because of the following two concerns: (1) Model validity and (2) authors’ conceptual misunderstanding about MLR analysis.

(1) Model validity: The model ability to represent observed variability is critical on judging if projected future changes are valid. The author must show 1) “additional model-observation comparisons” with choosing the variables which are relevant to this study’s focus (i.e., phytoplankton biomass) and 2) “evidence” on which results projected from the model with biases can be considered conclusive.

The authors showed the model-observation comparison of variability of annual mean phytoplankton carbon biomass (main target of this study) in Figure 1 and mentioned “Similar spatial patterns (to observation) are apparent (in the model)” in L139. But, for me, obviously, the model special pattern has different spatial characteristics from the observations. In the high latitudes, the observation shows the maximum variance in the most pole-side latitudes in the both hemispheres, while the model shows the variance maximum in somehow equator latitude around 50-60N and 50-60S. In the equator, although there is a strong latitudal maximum along the equator in the model, no such structure can not be seen in the observation, rather higher variabilities are observed in the off-equatorial regions. Moreover, model overestimations of the observed variability can exceed 200% in the equator and the subpolar North Atlantic.

The author also showed the model validity by comparing global internal variance in chlorophyll between the model ensemble and observational ensemble (Figure S2). However, this study’s focus is the phytoplankton carbon biomass, not chlorophyll, and these two can have very different spectrums. I think that model-observation comparisons in the biomass are more suitable for the purpose and the author should assess the model in the regional scale (not global), given the spatial heterogenous response of the biomass.

(2) MLR analysis: The methodology is unclear and it seems wrong.

The author tried to reconstruct the contribution of each driver variable to phytoplankton biomass using the MLR coefficients (equation 3 and 4). However, it obviously failed. As shown in Figure S6, the reconstructed “Total Carbon” is not equal to the sum of other terms (i.e., equation 3 and 4 are not correct), maybe because of inaccurate MLR coefficients, neglecting offset term or strong multicollinearity between variables (e.g., MLD and SST, SST and Solar).

Linear decompositions should be applied for “change/anomaly”, not for “climatology (10-year mean)”.

Given a function F(X,Y,Z), in general, the first order Taylor expansion is robust only for a small change in the F (ΔF),

ΔF = (∂F/∂X)ΔX+(∂F/∂Y)ΔY+(∂F/∂Z)ΔZ + (Residuals from high-order and cross terms).

The author should apply such analysis for “change” (not “climatology”) with considering residual terms. As the authors also mentioned, the partial differential coefficients are time-varying. The authors should be able to calculate the coefficients analytically using the model equations of phytoplankton carbon biomass.

Specific comments:

L29–: Any reference? And, does this mean that the CESM1 shows the opposite response of the high-latitude biomass to the global warming? (Figure 3a shows increase in biomass only in the sea-ice biome)

L49–: Please elaborate “Clarifying how variance in phytoplankton biomass may be changing over long time scales with climate change is important for fisheries management, especially at regional scales.”. What kind of impacts on fisheries by changing in variance in Phytoplankton biomass can one expect?

L82–85: I could not understand clearly. Please clarify with showing equations.

L94–97: The author’s description of the experimental setting of CESM1 large ensemble is inaccurate. Please describe it correctly.

L99–101: Show figure as an example.

L118–120: Please provide the map of the aggregated biological provinces used in this study as supplementally figure or superpose the biome boundary on the main figures (e.g., Figure 3).

Figure 1: Please use same colormap and same value range for fair comparison. And, it is better to show the ensemble mean of the σtemporal with a rank analysis (to show whether the observational σ is inside of the ensemble spread at grid by grid).

L179: Figure 2d?

L213–216: Which regions did the author chose? Please show these on map.

Technical corrections:

I don’t list any small technical/editorial corrections at this time. Above-mentioned conceptual/major comments should be addressed or fixed by the authors before going into the detail.

he manuscript, Anthropogenic climate change drives non-stationary phytoplankton variance, summarizes projected changes in global and regional phytoplankton variability using the NCAR CESM1 Large Ensemble under a high emissions scenario. The authors explore the key drivers of declining phytoplankton variability, highlighting the importance of top-down, zooplankton grazing in potentially driving future phytoplankton response. 

Generally, the article concisely represents its findings but there are several points of clarification I would recommend. In particular, the use of specific statistical terminology could be more accurate. Multiple times throughout the text, the term “variance” is used when, I think, “variability” is intended. In many cases this “variability” is being assessed via the standard deviation of the large ensemble members which is similar to the variance but not the same. Additionally, I am not proficient in MLR, but the comments made in the prior review are troubling especially considering the results are key to the paper's conclusions regarding top-down controls but these results seem underrepresented in the primary manuscript text. I’ve included several additional minor comments and suggestions below pertaining to clarity and organization. 

Specific Comments and Suggestions:

Lines 49-52: Clarifying how variance in phytoplankton biomass may be changing over long time scales with climate change is important for fisheries management, especially at regional scales. Near- term predictions of phytoplankton biomass may also benefit from knowledge of the projected magnitude of internal variability, as the chaotic nature of internal variability hampers near-term predictions (Meehl et al., 2009, 2014).

I think it’s worth noting that the internal variability quantified using a large ensemble is Internal variability specific to the model and indicative of our uncertainty that results from its simplified representations of the real world processes and numerics. It doesn’t necessarily have any bearing on real world manifestations of variability. Its primary utility to management and fisheries is in guiding our level of confidence in disentangling model signals from the noise. 

Lines 103-104: Six CESM1-LE members had corrupted ocean biogeochemistry

I’m curious, what does “corrupted ocean biochemistry” mean? it might help to explain what makes an ensemble member usable versus not. 

Figure 1. Add units: standard deviation should have the same units as the variable being assessed (i.e., phytoplankton carbon) but none appear in figure 1. 

Lines 121-122: Internal variability at each location (x,y) is approximated as the standard deviation across ensemble members (EMs) at a given time (t)

The method described here indicates that the standard deviation is being used to quantify variability. However, throughout the paper, the authors reference the “variance” when I think they mean “variability”. This is problematic because “variance” and “standard deviation”, while related, are two different values and the way they are interchanged throughout the text is confusing. Please check all instances of “variance” in the paper for intended meaning and replace with “variability” where appropriate. I suggest including a description of the “coefficient of variance” method here, too.

Lines 142-143: However, while the model ensemble captures regional patterns of observed variability, the CESM1-LE overestimates the magnitude of observed interannual variability.

I may be mistaken but it seems this was determined using only a single ensemble member - is it appropriate for conclusions to be drawn for the full ensemble when only considering one ensemble member?

Lines 147: A synthetic ensemble is a novel technique

I don’t think this technique can be called “novel” if it appears in two prior references

Lines 149-151: Compared to the internal variability over the observational period (2002 to 2020) (purple circle, (Figure S2), the model ensemble slightly overestimates the magnitude of internal variability in chlorophyll observed in the real world.

This seems like a result/ should appear in the result section. Also, it makes an assessment of the ensemble as a whole, but isn’t it still based on the results from the single ensemble member? If not, this was a point of confusion for me, and I suggest clarifying.

Lines 153-154: Annually averaged, global mean, upper-ocean (top 150m) integrated phytoplankton biomass across the model ensemble decreases from 76.1 mmol C m-2 to 66.2 mmol C m-2.

It’s not clear what timeframes these values represent. Is it 2006 vs. 2100? If so, it seems that such a narrow, 1-year window would risk aliasing higher frequency variability and potentially under- or overestimate the change in mean state. This is somewhat compensated for by the size of the ensemble but differs from the 10-year averaging described later in Line 223 Also, I suggest reporting the standard deviations for these numbers. 

Lines 177-178: we calculated the coefficient of variance as the standard deviation across the ensemble members for a given year (ensemble spread) divided by the ensemble mean.

I suggest including this description in the methods section rather than the results. 

Lines 178-180: Figure 2b illustrates the change in the coefficient of variance from the historical period through the RCP8.5 forcing scenario (1920 to 2100).

The results seem to jump from Figure 2a, to Figure 3, then back to 2b which is a bit confusing.

Line 180: The coefficient of variance is relatively constant across the historical period (1920 to 2005), and then significantly declines by ~20% from 2006-2100.

I’m not sure I agree with the assessment that the coefficient of variance is relatively constant across the historical period. 1920-1980 appears to have a positive trend with a range of about 6.1 to 7.3, which appears similar to the range of the time period covered by the dashed line in Figure 2b. I suggest testing the significance of the 1920-1980 trend. Also, could the drop in coefficient of variance instead be explained by temporal distance from the perturbation that differentiates the ensemble members? If the 34 ensemble members differ in initial air temperature conditions, would the spread perhaps be expected to decrease as the simulation integrates further away from that initial discrepancy (i.e., solutions start to converge)?

Lines 190-193: From 2006 to 2100, the coefficient of variance decreases by 3.3 x 10-5 yr-1 in the CESM1-LE, 2.0x10-4 yr in the MPI-ESM-LR1, 5.2x10-5 yr-1 in the CanESM2, and 3.9 x10-4 yr-1 in the GFDL-ESM2M. These declines are statistically significant in all model ensembles with the exception of the MPI-ESM-LR1 (Figure S2).

It’s not clear how these values across models are calculated, whether the end points of the time series or a range of years - the latter would be more appropriate (as done in Line 223) to avoid higher frequency variability and thus under- or overestimating the nature of the change. I also suggest reporting the specific statistical testing methods in the text if stating that the changes are significant. 

Line 201: We observe the largest magnitude decline in total phytoplankton carbon variance

The table is reporting change in standard deviation, not variance. Standard deviation is expressed in the same units as the analyzed variable while variance is reported in the square of those units.

Figure 4: It’s not clear what this figure adds to the discussion - it seems to be redundant with information in Figure 5. Perhaps if the outlines of the ecological regions were included? 

Lines 219-221: We quantified the relationship between phytoplankton carbon and the variables which contribute to changing phytoplankton biomass and its internal variability by performing a multiple linear regression (MLR) analysis. The MLR analysis was per- formed on linearly detrended annual anomalies using the ordinary least squares function of the Python package statsmodels.api

This and the associated equations seem to belong in the methods section. 

Line 274: …and important global biogeochemical regions…

What is considered an important biogeochemical region? This seems somewhat vague - I suggest elaborating to be a bit more specific. 

Lines 278-280: As such, the magnitude of changes in phytoplankton internal variance derived from the model ensemble should be interpreted as an overestimate when considering changes in phytoplankton internal variance driven by anthropogenic warming.

Again, my impression was that this conclusion was derived from analyzing a single ensemble member which seems insufficient for assessing the entire ensemble.

In this manuscript, the authors use the Community Earth System Model I Large Ensemble to evaluate the impacts of anthropogenic climate change on long-term variability in phytoplankton distributions within the global ocean. The authors additionally use a multiple linear regression to evaluate the ecological drivers of this change, reporting zooplankton grazing as being a major factor in reducing variability in phytoplankton biomass.

The analysis of earth systems models is well outside my area of expertise. So while the authors' main finding that variance in phytoplankton biomass is anticipated to decrease in the future ocean seems informative from my perspective, I defer to the first reviewer's comments regarding best practices in model interpretation. I was interested to see the multiple linear regression results, which seem to highlight a particularly strong coupling between phytoplankton biomass and grazing in model results. However, by the authors' admission on L265, it does not seem possible to establish cause and effect regarding the nature of this interaction. With this, it seems like an overstatement to suggest (as in the abstract and elsewhere) that these results provide evidence for grazing-driven declines in phytoplankton biomass.

More importantly, insufficient documentation is provided for the reader to interpret the MLR results. Critically, it is not immediately clear from the text how contributions to phytoplankton/diatom variance were calculated. Equations should be provided, and associated details on the MLR analysis should be moved to the methods section to make this information easier to locate in the manuscript. Moreover, the MLR results themselves seem insufficiently documented. No details are provided on the overall model fit nor on uncertainties associated with the MLR coefficients. The relationship between the parameters in equations 3 and 4 and the larger set of parameters included in figure 5 is unclear as well.

The discussion should also be expanded to provide more context on the authors' interpretation of these results. Altogether, even after reading the manuscript several times, I'm not sure why the results shouldn't be interpreted as a weakening of top-down control in the future ocean (with the decrease in contributions to phytoplankton biomass variance due to grazing in Figure 5 reflecting a reduced coupling of phytoplankton biomass and grazing and, by extension, a strengthening of bottom-up controls). If this interpretation is beyond what can be determined based on the analysis (for instance because of large uncertainties in coefficient errors), this is not evident from the information provided.

Without this information on the MLR results, it is impossible to critically evaluate some of the the manuscript's main conclusions. With this, and in light of the comments made by the first reviewer regarding issues with the authors' analysis of the CESM results, I cannot recommend this manuscript for publication without major revisions. A few specific comments are provided below.

Specific comments

L114 – 115 - A quick review of the method used in Tagliabue et al. would be useful here. What were the multivariate statistical methods used? How were they applied? A map of the biomes would be informative as well.

L159 – 161 - This text feels more appropriate for conclusion/discussion.

L215 - FAO citation and the associated reference seem to be improperly formatted

L220 — 221 - More information on how and why this transformation was performed would be useful.

L219 – 234 - This text feels more appropriate for the methods section

L289 – 291 - Is this conclusion inconsistent with the disclaimer provided at L264 – 266?

Equations 3 & 4 — Why are the terms in the equations (e.g., Solar, SST, Nutrient, etc.), different from those included in figure 5? Were the equations in the text just providing a summary of the actual equations used? If so, this should be made explicitly clear, with some description of all the variables included.

Figure 4 — Minor tick marks not necessary on color scale; difficult to see regions dominated by diazotrophs. Maybe use color palette with more contrast?

Figure 5 —Note inconsistent capitalization of biomes in subplots; Are units correctly labeled? Are the units for "contribution to phytoplankton/diatom variance" really mmol C m^-2? On a related note, where did the values on the Y axis come from? Based on the axis label they don't correspond to the MLR coefficients, but I didn't see any details in the text.