Articles | Volume 20, issue 21
Research article
10 Nov 2023
Research article |  | 10 Nov 2023

Anthropogenic climate change drives non-stationary phytoplankton internal variability

Geneviève W. Elsworth, Nicole S. Lovenduski, Kristen M. Krumhardt, Thomas M. Marchitto, and Sarah Schlunegger

Earth system models suggest that anthropogenic climate change will influence marine phytoplankton over the coming century with light-limited regions becoming more productive and nutrient-limited regions less productive. Anthropogenic climate change can influence not only the mean state but also the internal variability around the mean state, yet little is known about how internal variability in marine phytoplankton will change with time. Here, we quantify the influence of anthropogenic climate change on internal variability in marine phytoplankton biomass from 1920 to 2100 using the Community Earth System Model 1 Large Ensemble (CESM1-LE). We find a significant decrease in the internal variability of global phytoplankton carbon biomass under a high emission (RCP8.5) scenario and heterogeneous regional trends. Decreasing internal variability in biomass is most apparent in the subpolar North Atlantic and North Pacific. In these high-latitude regions, bottom-up controls (e.g., nutrient supply, temperature) influence changes in biomass internal variability. In the biogeochemically critical regions of the Southern Ocean and the equatorial Pacific, bottom-up controls (e.g., light, nutrients) and top-down controls (e.g., grazer biomass) affect changes in phytoplankton carbon internal variability, respectively. Our results suggest that climate mitigation and adaptation efforts that account for marine phytoplankton changes (e.g., fisheries, marine carbon cycling) should also consider changes in phytoplankton internal variability driven by anthropogenic warming, particularly on regional scales.

1 Introduction

Anthropogenic climate change significantly impacts marine ecosystems from phytoplankton (Bopp et al.2001, 2013; Laufkötter et al.2015; Kwiatkowski et al.2020) to fish (Perry et al.2005; Cheung et al.2009, 2010; Mills et al.2013; Wernberg et al.2016; Flanagan et al.2018; Staudinger et al.2019). As the base of the marine food web, phytoplankton support diverse marine ecosystems by providing food for higher trophic levels (Falkowski2012). Constraining future changes in phytoplankton with anthropogenic warming is important at regional scales for fishery adaptation (Pauly and Christensen1995; Chassot et al.2010; Link and Marshak2019; Marshak and Link2021), particularly as phytoplankton biomass is incorporated into offline fishery models to predict changing catch potential (Christensen and Walters2004; Travers-Trolet et al.2009; Lehodey et al.2010; Maury2010; Blanchard et al.2012; Christensen et al.2015; Jennings and Collingridge2015; Tittensor et al.2018; Petrik et al.2019; Heneghan et al.2021). In this context, understanding changes in both phytoplankton biomass and its internal variability is essential in reducing uncertainty in marine ecosystem projections.

The abundance and distribution of phytoplankton, the base of the marine food web and an important component of the marine carbon cycle, will likely change with anthropogenic warming. Future projections of climate change impacts reveal a global loss of marine net primary production (NPP) and phytoplankton biomass, particularly at middle and low latitudes (Steinacher et al.2010; Bopp et al.2013; Lotze et al.2019; Tittensor et al.2021). A majority of Earth system models (ESMs) project an increase in phytoplankton abundance in the high-latitude ocean as light limitation is alleviated by stratification, increasing temperature stimulates photosynthesis, and sea ice cover declines (Steinacher et al.2010; Bopp et al.2013). In contrast, a decrease in phytoplankton abundance in the low-latitude oceans is projected as nutrient limitation from thermal stratification is enhanced (Steinacher et al.2010; Kwiatkowski et al.2020). While bottom-up controls (e.g., nutrient flux, light availability) have been shown to affect phytoplankton growth in a changing climate, top-down controls (i.e., zooplankton grazing) also play a role. For example, analysis across a suite of models forced under climate change scenarios reveals grazing pressure as a driver of biomass decline in low- to intermediate-latitude regions (Laufkötter et al.2015). Additionally, top-down controls have been shown to affect regional changes in NPP and export production (Bopp et al.2001), as well as the timing of phytoplankton bloom onset (Yamaguchi et al.2022). Regional redistributions of phytoplankton biomass have consequences for fishery management and conservation (Blanchard et al.2017; Stock et al.2017) and may have implications for economics and policy making decisions (Moore et al.2021).

While climate change is known to impact the mean state of phytoplankton biomass or NPP (Bopp et al.2013; Kwiatkowski et al.2020), less is known about how climate change will affect internal variability in these quantities. One recent modeling study found that climate change alters the timing of seasonal blooms in many regions of the global ocean, an effect that could be realized by the end of the century (Yamaguchi et al., 2022). Several other recent studies have demonstrated how other aspects of the coupled atmosphere–ocean climate system are projected to experience changes in internal variability in a changing climate (Resplandy et al.2015; Landschützer et al.2018; Kwiatkowski and Orr2018; Rodgers et al.2021). For example, Resplandy et al. (2015) examined the contribution of internal variability to air–sea CO2 and O2 fluxes with climate change using a suite of six ESMs. Their analyses revealed distinct regional differences in internal variability of air–sea gas fluxes, with the Southern Ocean and the tropical Pacific playing a significant role. Other studies have revealed increases in the frequency of modes of internal variability such as El Niño and La Niña events in response to greenhouse warming (Timmermann et al.1999; Cai et al.2014, 2015, 2022). Clarifying how internal variability in phytoplankton biomass may be changing over long timescales with climate change is important for fishery management, especially at regional scales, as it affects our ability to make accurate near-term predictions of fishery production. 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).

Here, we quantify changes in the internal variability (ensemble spread) of phytoplankton biomass over the next century by using a large ensemble of an ESM in which each ensemble member experiences a different phasing of internal climate variability but is forced with a common emissions scenario. We illustrate the drivers of these changes in internal variability via statistical analysis of physical and biogeochemical model output and demonstrate their relative importance in key fishery regions.

2 Methods

2.1 Community Earth System Model 1 Large Ensemble

2.1.1 Model description

We evaluate changes in phytoplankton biomass internal variability using output from the Community Earth System Model 1 Large Ensemble (CESM1-LE) (Kay et al.2015). CESM1 is a fully coupled climate model that simulates Earth's climate under historical and Representative Concentration Pathway (RCP) 8.5 external forcing by simulating the evolution of coupled atmosphere, ocean, land, and sea ice component models (Hurrell et al.2013). The ocean physical model is the ocean component of the Community Climate System Model version 4 (Danabasoglu et al.2012) and has a nominal 1 resolution and 60 vertical levels. The Parallel Ocean Program version 2 (POP2) ocean model consists of an upper-ocean ecological module that incorporates multi-nutrient co-limitation of nitrate, ammonium, phosphate, dissolved iron, and silicate into phytoplankton growth and dynamic iron cycling (Moore et al.2004; Doney et al.2006; Moore and Braucher2008). The Biogeochemical Elemental Cycle (BEC) model simulates three phytoplankton functional types (PFTs): diatoms, diazotrophs, and small phytoplankton (i.e., cyanobacteria, nanophytoplankton, picoeukaryotes). Each PFT plays a unique role in the marine ecosystem and occupies a distinct ecological niche. For example, diatoms grow faster in cool, high-nutrient environments, while small phytoplankton thrive in warmer, low-nutrient environments. In contrast, diazotrophs are not limited by nitrogen availability due to their ability to biologically fix nitrogen from the atmosphere. Each PFT has a maximum growth rate, which is dictated by temperature (scaled by a temperature function with a Q10 of 2.0) and limited by nutrient and light availability (Moore et al.2004, 2013). Anthropogenic warming can alter these environmental variables and, in turn, affect phytoplankton abundance and productivity. Photoadaptation (variable chlorophyll to carbon ratios) occurs in response to variations in irradiance and nutrient availability (Geider et al.1998; Moore et al.2004). In addition to these bottom-up controls, top-down controls, such as zooplankton grazing, can also affect phytoplankton biomass. The ecosystem model simulates a single generic zooplankton functional type (ZFT) with different grazing rates and half-saturation constants prescribed for different PFTs (e.g., slower zooplankton grazing rates for larger phytoplankton, i.e., diatoms). Grazing rate is computed using a Holling Type III (sigmoidal) relationship and is a function of both prey density and temperature (Fig. S1, Eq. 5). Zooplankton loss is a function of a linear mortality term, which represents natural mortality, and a non-linear predation term, which represents losses from predation. Both of these loss terms scale with temperature. While zooplankton growth and loss terms both scale with temperature, a non-linear parameterization of the loss term results in a relatively larger increase in loss than increase in production with warming.

Large ensembles of ESMs are a recently developed research tool that allows us to disentangle fluctuations due to internal climate variability from those imposed by externally forced anthropogenic trends. Internal variability refers to variability in the climate system that occurs in the absence of external forcing and includes processes related to the coupled ocean–atmosphere system (e.g., Pacific Decadal Oscillation, El Niño–Southern Oscillation – ENSO) (Santer et al.2011; Deser et al.2010; Meehl et al.2013). In contrast, external forcing refers to the signal imposed by processes external to the climate system, such as solar variability, volcanic eruptions, and rising greenhouse gases from fossil fuel combustion (Deser et al.2012, 2010; Schneider and Deser2018). The CESM1-LE simulates the evolution of the climate system with multiple ensemble members, each initiated with slightly different atmospheric temperature fields and branched from a multi-century 1850 control simulation with constant pre-industrial forcing (Lamarque et al.2010; Danabasoglu et al.2012). The CESM1-LE simulates the evolution of the climate system from 1920 to 2100 with multiple ensemble members, each expressing a different phasing of internal climate variability while responding to a shared external forcing prescription (Kay et al.2015). Variable phasing of internal climate variability (e.g., ENSO) across ensemble members can influence phytoplankton biomass variability through the propagation of physical climate variability to biologically relevant environmental variables. RCP8.5 forcing was applied from 2006 to 2100 (Meinshausen et al.2011) with well-mixed greenhouse gases and short-lived aerosols projected by four different integrated assessment models (Lamarque et al.2010). A total of 40 ensemble members were generated for the CESM1-LE experiment. Of the 40, 6 CESM1-LE members had corrupted ocean biogeochemistry; therefore, we used the 34 CESM1-LE members with valid ocean biogeochemistry.

2.1.2 Statistical analysis of model output

Analyses were conducted using annual mean output at 1 resolution from 1920 to 2100. Changes in CESM1 phytoplankton internal variability can be assessed via statistical analysis of chlorophyll concentration, net primary productivity (NPP), or phytoplankton carbon concentration (an indicator of total biomass). In our analysis we focus on biomass (phytoplankton carbon concentration) because it is an important predictor variable in offline fishery models (Christensen and Walters2004; Travers-Trolet et al.2009; Lehodey et al.2010; Maury2010; Blanchard et al.2012; Christensen et al.2015; Jennings and Collingridge2015; Tittensor et al.2018). Additionally, under climate change scenarios, phytoplankton biomass may be a more reliable indicator of climate change impacts than NPP (Bopp et al.2022). Vertical integrals (top 150 m) of biomass carbon concentration from each PFT were calculated and then summed to create maps of total phytoplankton biomass.

Figure 1Comparison between observed and modeled phytoplankton biomass interannual variability. (a) Temporal standard deviation in annual mean phytoplankton carbon concentration (mg m−3) reconstructed from remotely sensed chlorophyll concentrations, backscattering coefficients, and phytoplankton absorption (1998 to 2019) (Bellacicco et al.2020). (b) Temporal standard deviation in annual mean phytoplankton carbon concentration (mg m−3) simulated by ensemble member 1 of the CESM1-LE over the same observational period (1998 to 2019). Note the different magnitudes on the color bars.

We classified the marine environment into 11 ecologically cohesive biomes as in Tagliabue et al. (2021) and Vichi et al. (2011) (Fig. S2), which are a consolidation of the 38 ecological regions defined in Longhurst (2007). The provinces were aggregated using multivariate statistical analysis of physical (i.e., salinity, temperature, mixed layer depth) and biological (i.e., chlorophyll concentration) ocean parameters to group ocean regions with similar physical and environmental conditions (Vichi et al.2011). The ocean provinces were defined by randomly selecting from a combination of model and observational datasets and testing for statistical significance using analysis of similarities (ANOSIM) (Vichi et al.2011). Although we consider all 11 biomes in our analysis, we analyze drivers in 4 biomes that are particularly relevant for fishery production and/or of high biogeochemical interest: the subpolar Atlantic (ASP), the subarctic Pacific (SAP), the equatorial Pacific (EQP), and the Southern Ocean (SOC) (Fig. S2). ASP is a consolidation of aggregated biogeochemical provinces 4, 11, and 15; SAP a consolidation of 50 and 51; EQP a consolidation of 61, 62, and 63; and SOC a consolidation of 21, 81, 82, and 83 (Longhurst2007; Vichi et al.2011) (Fig. S2). Important biogeochemical regions are those characterized by coherent physical and environmental conditions, which support unique marine ecosystems and play an outsized role in global ocean biogeochemistry.

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

(1) σ ( x , y , t ) = σ ( EM ( x , y , t ) ) .

The coefficient of variation (CoV) is calculated as the standard deviation across the ensemble members divided by the ensemble mean,

(2) CoV ( x , y , t ) = σ ( EM ( x , y , t ) ) LE EM .

The forced response of the large ensemble is calculated as the mean of ensemble members at a given location and time,

(3) LE ( x , y , t ) = 1 n EM ( x , y , t ) n ,

where n is the number of ensemble members.

We quantified the drivers of phytoplankton carbon biomass CoV in key ocean regions by generating an ensemble of boosted regression trees. Unlike linear models, boosted trees are able to capture non-linear interaction between the predictors and the response and have been used in a number of ecological applications (Elith et al.2008; Roberts et al.2016; Lamb et al.2021; Dannouf et al.2022; Denvil-Sommer et al.2023). A regression tree ensemble is a predictive model composed of a weighted combination of multiple regression trees. At every step, the ensemble fits a new learner to the difference between the observed response and the aggregated prediction of all learners grown previously, aiming to minimize mean-squared error. We generate an ensemble of boosted regression trees (maximum tree depth of 10) using the MATLAB function fitrensemble. Our predictor variables are the regional mean, ensemble mean temperature, mixed layer depth, incoming shortwave radiation, physically mediated iron, physically mediated phosphate, zooplankton carbon, and zooplankton grazing (diatom, small phytoplankton, or their sum) annually resolved from 2006 to 2100, while our response variable is CoV of phytoplankton carbon (diatom, small phytoplankton, or their sum) annually resolved from 2006 to 2100. We use the MATLAB function predictorImportance to estimate the importance of the predictors for each tree learner in the ensemble; it computes the importance of the predictors in a tree by summing changes due to splits on every predictor and dividing the sum by the total number of branches. The machine learning model was tuned to a learning rate of 1 and a tree depth of 10, generating 100 trees. We tuned several hyperparameters to generate the highest-quality predictive results with the least computational expense. While learning rate can affect the quality of the solution, we experimented with a range of learning rates (0.1–1) with no change in the predictive results. Similarly, we tuned the tree depth using a range of 1 to 10 splits, and tree depths less than 10 produced a higher RMSE of the testing dataset.

2.2 Model evaluation

We used remotely sensed estimates of phytoplankton carbon to evaluate the representation of phytoplankton interannual variability in the CESM1-LE. In other words, we evaluate the temporal variability in modeled phytoplankton biomass from year to year. We note that this interannual variability is different than the internal variability (ensemble spread) that we discuss at length in this study, but is nevertheless a target for model validation. Although phytoplankton carbon concentrations cannot be measured directly by satellites, they can be reconstructed using algorithms that incorporate remotely sensed chlorophyll concentrations, detrital backscattering coefficients, and phytoplankton absorption (Kostadinov et al.2016; Martinez-Vicente et al.2017; Roy et al.2017; Sathyendranath et al.2020; Brewin et al.2021). We use the observational phytoplankton carbon dataset of Bellacicco et al. (2020), annually averaged and interpolated onto a 1 grid, to evaluate interannual variability in phytoplankton biomass in a single model ensemble member. Figure 1a shows satellite-derived estimates of interannual variability in phytoplankton carbon with regions of relatively low phytoplankton variability shown in yellow and regions of relatively high variability in purple. Remotely sensed observations capture areas of high interannual variability in the subpolar North Atlantic, North Pacific, and Southern Ocean and areas of low interannual variability in the subtropical gyre regions. Similar spatial patterns are apparent when compared to the range of phytoplankton interannual variability in ensemble member 1 of the CESM1-LE over the observational period (1998 to 2019) (Fig. 1b). However, while the model ensemble captures regional patterns of observed variability, the CESM1-LE overestimates the magnitude of observed interannual variability. Some regions of the global ocean display a substantial mismatch in interannual variability between the model and that estimated from observations (Fig. 1, Table S1). While the differences can be quite large in some regions, we note that this is an evaluation of interannual variability (rather than internal variability, the focus of this study) and that estimates from the satellite product are derived from a collection of data products that may also display biases (Table S1).

As an evaluation of the model's ability to represent internal variability (ensemble spread), we compare the internal variability in chlorophyll simulated in the CESM1-LE to a synthetic ensemble generated from observed surface chlorophyll concentrations over the MODIS remote sensing record (Elsworth et al.2020, 2021) (Fig. S3; chlorophyll was readily available in the CESM1-LE and can be directly compared with our synthetic ensemble of observed surface chlorophyll). A synthetic ensemble is a technique that allows the observational record to be statistically emulated to create multiple possible evolutions of the observed record, each with a unique sampling of internal climate variability (McKinnon et al.2017; McKinnon and Deser2018). Compared to the internal variability over the observational period (2002 to 2020) (purple circle, Fig. S3), the model ensemble slightly overestimates the magnitude of internal variability in chlorophyll observed in the real world.

Taken together, our model validation exercises demonstrate that the model tends to overestimate both the temporal (interannual) variability and the internal variability in phytoplankton, as compared to satellite observations on both global and regional scales. Thus, we must interpret our findings with this caveat in mind.

Figure 2(a) Global change in annual mean total phytoplankton carbon concentration simulated by the CESM1-LE (mmol C m−2) from the historical period through the RCP8.5 forcing scenario (1920 to 2100). The ensemble mean is shown in the black curve and the 34 individual ensemble members are shown in the gray curve. (b) Global change in the coefficient of variation in annual mean total phytoplankton carbon concentration over the same period, smoothed using a 5-year window. The trend in the coefficient of variation over the RCP8.5 forcing scenario is shown with the black dashed line.


Figure 3(a) Percentage change in annual total phytoplankton carbon concentration over the RCP8.5 forcing scenario (2006 to 2100) simulated by the CESM1-LE. (b) Percentage change in annual total phytoplankton internal variability over the same period. The change in the mean and the variability is calculated using averages across the first (2006 to 2016) and last (2090 to 2100) decades of the RCP8.5 forcing scenario. Hatched areas indicate regions of trend insignificance determined by a t test with a p value greater than 0.05. Summary statistics for the t test are available in the Supplement (Table S2).

3 Results

We evaluate the change in mean phytoplankton biomass and its internal variability across the CESM1-LE globally and regionally. Annually averaged, global mean, upper-ocean (top 150 m) integrated phytoplankton biomass across the model ensemble decreases from 76.1 mmol C m−2 to 66.2 mmol C m−2 from the historical period through the RCP8.5 forcing scenario (1920 to 2100), a decline of 13 % (black curve; Fig. 2a). The change in the mean is calculated as the difference between the first (1920 to 1930) and last (2090 to 2100) decades across the historical and RCP8.5 forcing scenario. Phytoplankton biomass declines globally, except in polar regions (Fig. 3a). Regional changes in mean phytoplankton biomass across the RCP8.5 forcing scenario (2006 to 2100) display increasing biomass in portions of the Arctic and the Southern Ocean that gradually become ice free over the century (on the order of 20 %–40 % of the mean biomass across the century) and decreasing biomass across the subtropical gyres (on the order of 15 %–30 % of the mean biomass across the century; Figs. 3a, S4a). In the North Atlantic subpolar gyre, the phytoplankton biomass declines by 40 %–50 % of its mean (Figs. 3a, S4a). This result is consistent with previous modeling studies that identified a 50 % reduction in North Atlantic primary production associated with AMOC (Atlantic Meridional Overturning Circulation) weakening during the last glacial period (Schmittner2005). A weakening of the AMOC is also projected with anthropogenic warming (Manabe and Ronald1993; Stocker and Schmittner1997).

Regional changes in phytoplankton biomass are dominated by changes in diatom and small phytoplankton (Table 1). We aggregate biomass across 11 ecological provinces (Vichi et al.2011; Tagliabue et al.2021) and present changes in total and PFT biomass over the RCP8.5 scenario in Table 1. The CESM1-LE simulates the largest decline in total phytoplankton carbon concentration in the Atlantic subpolar (ASP) region, where diatom biomass declines by  80 mmol C m−2 and small phytoplankton biomass increases slightly ( 8 mmol C m−2). We observe moderate decreases in the subpolar Pacific (SAP) region that are again driven by declines in diatom carbon concentration, with minor contributions from changes in small phytoplankton carbon concentration (Table 1). The CESM1-LE simulates a smaller decline in total carbon concentration in the EQP region, where diatom biomass declines by  7 mmol C m−2 and small phytoplankton biomass declines by  5 mmol C m−2. The smallest decline in total carbon concentration occurs in the South Pacific subtropical gyre (SPS) region, where diatom biomass declines by  4.3 mmol C m−2 and small phytoplankton biomass declines by  4.6 mmol C m−2.

Internal variability in global phytoplankton biomass, which is indicated by the spread across the individual ensemble members (gray lines; Fig. 2a), declines over the RCP8.5 forcing period from 2006 to 2100. To quantify how the range of internal variability in phytoplankton biomass is changing with anthropogenic warming, we calculated the coefficient of variation as the standard deviation across the ensemble members for a given year (ensemble spread) divided by the ensemble mean. Figure 2b illustrates the change in the coefficient of variation from the historical period through the RCP8.5 forcing scenario (1920 to 2100). The coefficient of variation is relatively constant across the historical period (1920 to 2005) and then significantly declines by  20 % from 2006 to 2100.

A decrease in global phytoplankton internal variability with anthropogenic warming is not unique to the CESM1-LE. We illustrate this by analyzing surface phytoplankton chlorophyll (rather than biomass as surface chlorophyll was readily available in the CMIP5 archive) from three other CMIP5 ESM large ensembles that include representation of ocean biogeochemistry: the GFDL-ESM2M from the Geophysical Fluid Dynamics Laboratory (GFDL; Dunne et al.2012, 2013), the CanESM2 from the Canadian Centre for Climate Modelling and Analysis (Christian et al.2010; Arora et al.2011), and the MPI-ESM-LR from the Max Planck Institute for Meteorology (MPI; Giorgetta et al.2013; Ilyina et al.2013). These ensembles consist of 30, 50, and 100 ensemble members, respectively. Similarly to the CESM1-LE, historical forcing was applied through 2005, followed by RCP8.5 forcing through 2100. While there is substantial spread in the mean coefficient of variation across the four models, a similar decline in the coefficient of variation can be observed across each of the four ESM ensembles (Fig. S3). From 2006 to 2100, the coefficient of variation decreases by 3.3 × 10−5 yr−1 in the CESM1-LE, 2.0 × 10−4 yr−1 in the MPI-ESM-LR1, 5.2 × 10−5 yr−1 in the CanESM2, and 3.9 × 10−4 yr−1 in the GFDL-ESM2M. The change in the coefficient of variation is calculated using averages across the first (2006 to 2016) and last (2090 to 2100) decades of the RCP8.5 forcing scenario. These declines are statistically significant in all model ensembles with the exception of the MPI-ESM-LR1 (Fig. S3).

In comparison to the mean change in phytoplankton biomass, changes in phytoplankton internal variability with time are spatially more heterogeneous across the global ocean (Fig. 3b). The largest decreases in internal variability are apparent in the North Atlantic and North Pacific subpolar regions (on the order of 50 %–70 % of the mean biomass internal variability), with smaller declines in the equatorial Pacific and Southern oceans (on the order of 30 %–50 % of the mean biomass internal variability) (Figs. 3b, S4b). Changes in internal variability in the subtropical regions are characterized by mixed trends with areas of both increasing and decreasing internal variability across the RCP8.5 forcing scenario (Figs. 3b, S4b).

Global changes in total phytoplankton biomass standard deviation are a manifestation of changes in diatom and small phytoplankton variability (Table 1). We observe the largest-magnitude decline in total phytoplankton carbon standard deviation in the subpolar Atlantic (ASP) region, where diatom standard deviation declines by  10 mmol C m−2 and small phytoplankton standard deviation declines by  2 mmol C m−2 (Table 1). The CESM1-LE simulates a moderate-magnitude decline in total phytoplankton standard deviation in the subarctic Pacific (SAP) region, driven by a decrease in small phytoplankton standard deviation ( 2 mmol C m−2) with minor contributions from declines in diatom standard deviation ( 1 mmol C m−2) (Table 1). Moderate declines in standard deviation are also simulated in the Arctic (ARC), North Atlantic subtropical gyre (NAS), SOC, and EQP regions, driven by declines in diatom carbon standard deviation in the SOC region and declines in small phytoplankton internal variability in the EQP region (Table 1).

Table 1Changes in phytoplankton biomass and its internal variability in the CESM1-LE from 2006 to 2100 for the 11 ecological provinces defined in Vichi et al. (2011) and Tagliabue et al. (2021) (mmol C m−2). The change in the mean and standard deviation is calculated using averages across the first (2006 to 2016) and last (2090 to 2100) decades of the RCP8.5 forcing scenario.

Download Print Version | Download XLSX

To guide our analysis of changing phytoplankton biomass internal variability, we considered the dominant ecological assemblage across different regions of the global ocean. The CESM1-LE simulates three phytoplankton functional types, each of which thrive in distinct regions of the global ocean. Diatoms dominate in the subpolar Atlantic and Pacific, the eastern equatorial upwelling zone, and portions of the Southern Ocean, while small phytoplankton dominate across the subtropical gyres and portions of the Southern Ocean (Fig. 4). In contrast, diazotrophs, a minor contributor to total carbon biomass, are present at such low concentrations that they do not dominate anywhere in the global ocean (Fig. 4). Using the ecologically cohesive regions defined by Tagliabue et al. (2021) and Vichi et al. (2011), we selected areas that align with the most productive fishery regions by catch in the Atlantic and Pacific basins (FAO2020), as well as regions of global biogeochemical importance for further analysis. In each ecological region we identified the dominant phytoplankton functional type to include in our analysis. In regions where multiple phytoplankton functional types dominated, we used total carbon concentrations to reflect the mixed ecological assemblage.

Figure 4Distribution of the dominant phytoplankton functional type in biomass carbon averaged across the RCP8.5 forcing scenario (2006 to 2100). The CESM1-LE simulates three phytoplankton functional types: diatoms, diazotrophs, and small phytoplankton. Regions where diatoms dominate are shown in yellow, and regions where small phytoplankton dominate are shown in purple. Diazotrophs do not dominate in any region of the global ocean. The four ecological provinces shown are the subpolar Pacific (SAP), the subpolar Atlantic (ASP), the equatorial Pacific (EQP), and the Southern Ocean (SOC).

We identify the importance of different predictors to changing phytoplankton biomass CoV in four distinct ecological regions using a machine learning (boosted regression tree) approach. In the subpolar Atlantic (ASP) and subpolar Pacific (SAP) ecological provinces (Fig. 4), diatom biomass CoV declines between the beginning and end of the century (Table 1). In the Atlantic subpolar region, the most important predictor of diatom biomass CoV is phosphate advection, with smaller contributions from zooplankton carbon (Fig. 5a). In the subarctic Pacific region, sea surface temperature is the most important predictor of diatom biomass CoV, with phosphate advection playing a secondary role (Fig. 5b).

As the SOC and EQP ecological provinces are characterized by mixed phytoplankton assemblages where both diatoms and small phytoplankton dominate (Fig. 4), we identify the predictors of total phytoplankton CoV here. In contrast to the subpolar Atlantic and subpolar Pacific provinces, we observe a relatively smaller decline in phytoplankton CoV between the beginning and end of the century in the Southern Ocean (Table 1). The most important predictors of phytoplankton CoV in the SOC region derive from solar flux, with more minor contributions from iron and phosphate advection (Fig. 5c). In the equatorial Pacific region, zooplankton carbon is the most important predictor of total phytoplankton CoV, while iron and phosphate advection play less of a predictive role (Fig. 5d).

In all four ecological provinces, a combination of bottom-up controls (e.g., nutrient supply, light availability) and top-down controls (e.g., grazer biomass) predicts the decline in phytoplankton biomass CoV with anthropogenic warming. Our statistical analysis reveals that phosphate advection is an important predictor in the high-latitude regions of both the subpolar Atlantic and Pacific, with sea surface temperature playing an important role in the subpolar Pacific. However, in the Southern Ocean and the equatorial Pacific, solar flux and grazer biomass dominate the predictive skill in phytoplankton biomass CoV.

Figure 5Relative importance of predictor variables for phytoplankton biomass coefficient of variation across the RCP8.5 forcing scenario (2006 to 2100). Marine ecological regions are defined in Tagliabue et al. (2021). Regions were selected that aligned with the highest fishery catch in the (a) Atlantic and (b) Pacific basins and the biogeochemically important (c) Southern Ocean and (d) equatorial Pacific regions. The dominant phytoplankton functional type is considered in each region. In regions with a mixed ecological assemblage, total phytoplankton carbon is considered. The RMSE (mmol C m−2) for the testing dataset of each machine learning analysis is included in the upper right corner of each panel.


4 Conclusions and discussion

We quantify both global and regional changes in phytoplankton internal variability across the RCP8.5, a business-as-usual forcing scenario, in the CESM1-LE. We observe a global decline in phytoplankton internal variability in the model ensemble, which is reflected in similar declines in phytoplankton internal variability across a suite of CMIP5 models (Fig. S3). Regional changes in phytoplankton variability with anthropogenic climate change in the model ensemble are spatially heterogeneous, with highly productive fishery regions and important global biogeochemical regions experiencing large changes in internal variability. Using a machine learning approach, we identify the importance of different predictors to changing phytoplankton biomass internal variability. In all four ecological provinces, a combination of bottom-up controls (e.g., nutrient supply, light availability) and top-down controls (e.g., grazer biomass) predicts the decline in phytoplankton biomass CoV with anthropogenic warming.

While the CESM1-LE represents the overall spatial pattern of observed interannual variability in phytoplankton carbon, the model overestimates the magnitude of observed interannual and internal variability in phytoplankton on regional scales. This caveat is particularly important to consider when interpreting projections from offline fishery models in the context of fishery adaptation and planning in a warming climate.

Our statistical analysis approach has inherent limitations, especially in the context of a attributing changes in an inherently coupled system (i.e., one in which predictor variables co-vary). In a coupled system such as this, it is difficult to definitively identify cause and effect. In this context, the statistical method can be used as an effective tool to provide a first-order approximation of contributions to changing phytoplankton CoV.

While many studies attribute bottom-up controls to changing phytoplankton with anthropogenic warming (Steinacher et al.2010; Bopp et al.2013; Lotze et al.2019; Tittensor et al.2021), top-down controls may also play an important role, particularly in our understanding of changing phytoplankton biomass and its internal variability. Our study demonstrates a connection between phytoplankton internal variability and zooplankton carbon in the subpolar North Pacific and equatorial Pacific. Previous studies of phytoplankton change with climatic warming have demonstrated that grazing pressure, the fraction of phytoplankton biomass grazed, is a contributor to biomass decline in low- to intermediate-latitude regions across a suite of model simulations with different marine ecosystem models (Laufkötter et al.2015) and that top-down controls can affect regional changes in NPP and export production (Bopp et al.2001) and are a contributor to future shifts in bloom timing (Yamaguchi et al.2022). While grazing pressure has been shown to increase in response to climate change, several ecosystem models have also identified zooplankton grazing as a dominant contributor to phytoplankton assemblage succession during blooms (Hashioka et al.2013; Prowe et al.2012a). Additionally, top-down controls have also been observed to affect the onset of the spring bloom (Behrenfeld2010; Behrenfeld et al.2013) and to influence primary production in a trait-based ecosystem model (Prowe et al.2012b).

The relative simplicity of the ocean biogeochemical ecosystem model in CESM1 (e.g., representation of a single zooplankton functional type with multiple grazing rates) may limit a more detailed evaluation of changing grazing pressures with climate change. While the recent parameterization of the biogeochemical ecosystem model in CESM2 (MARBL – Marine Biogeochemistry Laboratory) includes similar representation of three PFTs and a single adaptive ZFT (Long et al.2021), more complex configurations of MARBL include explicit representation of additional PFTs such as coccolithophores (Krumhardt et al.2019) and ZFTs. Additional insights into contributions to internal variability may be gained using more complex models. Additionally, the use of an ecosystem model of higher complexity may provide more realistic projections of the marine ecosystem with climate change considering change in phytoplankton and zooplankton species diversity with anthropogenic warming (Benedetti et al.2021).

The magnitude and direction of regional changes in phytoplankton internal variability are an essential constraint for near-term (subseasonal to decadal) predictions of the local marine ecosystem, particularly in important fishery regions such as the subpolar Atlantic (ASP) and the subpolar Pacific (SAP) ecological provinces (FAO2020). Accurate near-term predictions require foreknowledge of both internal climate variability and external climate change signals. On subseasonal to decadal timescales, the magnitude of internal climate variability is often stronger than forced climate change signals (Meehl et al.2009, 2014). In this context, a decline in phytoplankton internal variability with anthropogenic climate change may improve the accuracy of near-term predictions of phytoplankton biomass, producing more reliable forecasts of fishery productivity and marine carbon cycling. Future work can utilize these constraints on phytoplankton internal variability, particularly on regional scales, to inform climate mitigation and adaptation efforts.

Data availability

CESM1-LE output is available from the Earth System Grid at (last access: 6 November 2023,, Kay and Deser2016).


The supplement related to this article is available online at:

Author contributions

GWE and NSL conceptualized the study. GWE analyzed simulation results and prepared the manuscript. NSL assisted in preparing the manuscript. KMK assisted in analysis of simulation results. TMM, KMK, and SS assisted in reviewing the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


Computational resources were provided by the Computational and Information Systems Laboratory (CISL) at the National Center for Atmospheric Research (NCAR), through a resource allocation to Geneviève W. Elsworth and Nicole S. Lovenduski.

Financial support

This research has been supported by the National Science Foundation (grant nos. OCE-1558225 and OCE-1752724).

Review statement

This paper was edited by Ciavatta Stefano and reviewed by Nicholas Bock and five anonymous referees.


Arora, V., Scinocca, J., Boer, G., Christian, J., Denman, K., Flato, G., Kharin, V., Lee, W., and Merryfield, W.: Carbon emission limits required to satisfy future representative concentration pathways of greenhouse gases, Geophys. Res. Lett., 38, L05805,, 2011. a

Behrenfeld, M.: Abandoning Sverdrup's Critical Depth Hypothesis on phytoplankton blooms, Ecology, 91, 977–89,, 2010. a

Behrenfeld, M., Doney, S., Lima, I., Boss, E., and Siegel, D.: Annual cycles of ecological disturbance and recovery underlying the subarctic Atlantic spring plankton bloom: PHYTOPLANKTON BLOOMS, Global Biogeochem. Cy., 27, 526–540,, 2013. a

Bellacicco, M., Pitarch, J., Organelli, E., Martinez-Vicente, V., Volpe, G., and Marullo, S.: Improving the Retrieval of Carbon-Based Phytoplankton Biomass from Satellite Ocean Colour Observations, Remote Sens., 12, 3640,, 2020. a, b

Benedetti, F., Vogt, M., Hofmann Elizondo, U., Righetti, D., Zimmermann, N., and Gruber, N.: Major restructuring of marine plankton assemblages under global warming, Nat. Commun., 12, 5226,, 2021. a

Blanchard, J., Jennings, S., Holmes, R., Harle, J., Merino, G., Allen, I., Holt, J., Dulvy, N., and Barange, M.: Potential consequences of climate change on primary production and fish production in large marine ecosystems, Philos. T. R. Soc. B, 367, 2979–2989,, 2012. a, b

Blanchard, J., Watson, R., Fulton, E., Cottrell, R., Nash, K., Bryndum-Buchholz, A., Büchner, M., Carozza, D., Cheung, W., Elliott, J., Davidson, L., Dulvy, N., Dunne, J., Eddy, T., Galbraith, E., Lotze, H., Maury, O., Müller, C., Tittensor, D., and Jennings, S.: Linked sustainability challenges and trade-offs among fisheries, aquaculture and agriculture, Nature, 1, 1240–1249,, 2017. a

Bopp, L., Monfray, P., Aumont, O., Dufresne, J.-L., Treut, H., Madec, G., Terray, L., and Orr, J.: Potential impact of climate change on marine export production, Global Biogeochem. Cy., 15, 81–100,, 2001. a, b, c

Bopp, L., Resplandy, L., Orr, J., Doney, S., Dunne, J., Gehlen, M., Halloran, P., Heinze, C., Ilyina, T., Séférian, R., Tjiputra, J., and Vichi, M.: Multiple stressors of ocean ecosystems in the 21st century: Projections with CMIP5 models, Biogeosciences, 10, 6225–6245,, 2013. a, b, c, d, e

Bopp, L., Aumont, O., Kwiatkowski, L., Clerc, C., Dupont, L., Ethé, C., Gorgues, T., Séférian, R., and Tagliabue, A.: Diazotrophy as a key driver of the response of marine net primary productivity to climate change, Biogeosciences, 19, 4267–4285,, 2022. a

Brewin, B., Sathyendranath, S., Platt, T., Bouman, H., Ciavatta, S., Dall'Olmo, G., Dingle, J., Groom, S., Jönsson, B., Kostadinov, T., Kulk, G., Laine, M., Martinez-Vicente, V., Psarra, S., Raitsos, D., Richardson, K., Rio, M.-H., Rousseaux, C., Salisbury, J., and Walker, P.: Sensing the ocean biological carbon pump from space: A review of capabilities, concepts, research gaps and future developments, Earth-Sci. Rev., 217, 103604,, 2021. a

Cai, W., Borlace, S., Lengaigne, M., Rensch, P., Collins, M., Vecchi, G., Timmermann, A., Santoso, A., McPhaden, M., Wu, L., England, M., Wang, G., Guilyardi, E., and Jin, F.-F.: Increasing frequency of extreme El Niño Events due to greenhouse warming, Nat. Clim. Change, 4, 111–116,, 2014. a

Cai, W., Wang, G., Santoso, A., McPhaden, M., Wu, L., Jin, F.-F., Timmermann, A., Collins, M., Vecchi, G., Lengaigne, M., England, M., Dommenget, D., Takahashi, K., and Guilyardi, E.: Increased frequency of extreme La Niña Events under greenhouse warming, Nat. Clim. Change, 5, 132–137,, 2015. a

Cai, W., Ng, B., Wang, G., Santoso, A., Wu, L., and Yang, K.: Increased ENSO sea surface temperature variability under four IPCC emission scenarios, Nat. Clim. Change, 12, 228–231,, 2022. a

Chassot, E., Bonhommeau, S., Dulvy, N., Mélin, F., Watson, R., and le Pape, O.: Global marine primary production constrains fisheries catches, Ecol. Lett., 13, 495–505,, 2010. a

Cheung, W., Lam, V., Sarmiento, J., Kearney, K., Watson, R., and Pauly, D.: Projecting Global Marine Biodiversity Impacts under Climate Change Scenarios, Fish Fish., 10, 235–251,, 2009. a

Cheung, W., Lam, V., Sarmiento, J., Kearney, K., Watson, R., Zeller, D., and Pauly, D.: Large-scale Redistribution of Maximum Fisheries Catch Potential in the Global Ocean under Climate Change, Glob. Change Biol., 16, 24–35,, 2010. a

Christensen, V. and Walters, C.: Ecopath With Ecosim: Methods, Capabilities and Limitations, Ecol. Model., 172, 109–139,, 2004. a, b

Christensen, V., Coll, M., Buszowski, J., Cheung, W., Frölicher, T., Steenbeek, J., Stock, C., Watson, R., and Walters, C.: The global ocean is an ecosystem: Simulating marine life and fisheries, Glob. Ecol. Biogeogr., 24, 507–517,, 2015. a, b

Christian, J., Arora, V., Boer, G., Curry, C., Zahariev, K., Denman, K., Flato, G., Lee, W., Merryfield, W., Roulet, N., and Scinocca, J.: The global carbon cycle in the Canadian Earth system model (CanESM1): Preindustrial control simulation, J. Geophys. Res.-Biogeo., 115, G03014,, 2010. a

Danabasoglu, G., Bates, S. C., Briegleb, B. P., Jayne, S. R., Jochum, M., Large, W. G., Peacock, S., and Yeager, S. G.: The CCSM4 ocean component, J. Clim., 25, 1361–1389,, 2012. a, b

Dannouf, R., Yong, B., Ndehedehe, C., Correa, F., and Ferreira, V.: Boosted Regression Tree Algorithm for the Reconstruction of GRACE-Based Terrestrial Water Storage Anomalies in the Yangtze River Basin, Front. Environ. Sci., 10, 917545,, 2022. a

Denvil-Sommer, A., Buitenhuis, E., Kiko, R., Fabien, L., Guidi, L., and Le Quéré, C.: Testing the reconstruction of modelled particulate organic carbon from surface ecosystem components using PlankTOM12 and machine learning, Geosci. Model Dev., 16, 2995–3012,, 2023. a

Deser, C., Phillips, A., Bourdette, V., and Teng, H.: Uncertainty in climate change projections: The role of internal variability, Clim. Dynam., 38, 527–546,, 2010. a, b

Deser, C., Knutti, R., Solomon, S., and Phillips, A.: Communication of the role of natural variability in future North American climate, Nat. Clim. Change, 2, 775–779,, 2012. a

Doney, S., Lindsay, K., Fung, I., and John, J.: Natural variability in a stable, 1000-Yr global coupled climate–carbon cycle simulation, J. Clim., 19, 3033–3054,, 2006. a

Dunne, J., John, J., Adcroft, A., Griffies, S., Hallberg, R., Shevliakova, E., Ronald, S., Cooke, W., Dunne, K., Harrison, M., Krasting, J., Malyshev, S., Milly, P., Phillips, P., Sentman, L., Samuels, B., Spelman, M., Winton, M., Wittenberg, A., and Zadeh, N.: GFDL's ESM2 Global Coupled Climate–Carbon Earth System Models, Part I: Physical Formulation and Baseline Simulation Characteristics, J. Clim., 25, 6646–6665,, 2012. a

Dunne, J., John, J., Shevliakova, E., Ronald, S., Krasting, J., Malyshev, S., Milly, P., Sentman, L., Adcroft, A., Cooke, W., Dunne, K., Griffies, S., Hallberg, R., Harrison, M., Levy, H., Wittenberg, A., Phillips, P., and Zadeh, N.: GFDL's ESM2 Global Coupled Climate–Carbon Earth System Models, Part II: Carbon System Formulation and Baseline Simulation Characteristics, J. Clim., 26, 2247–2267,, 2013. a

Elith, J., Leathwick, J., and Hastie, T.: A Working Guide to Boosted Regression Trees, J. Anim. Ecol., 77, 802–13,, 2008. a

Elsworth, G., Lovenduski, N., McKinnon, K., Krumhardt, K., and Brady, R.: Finding the Fingerprint of Anthropogenic Climate Change in Marine Phytoplankton Abundance, Curr. Clim. Change Rep., 6 37–46,, 2020. a

Elsworth, G., Lovenduski, N., and McKinnon, K.: Alternate History: A Synthetic Ensemble of Ocean Chlorophyll Concentrations, Global Biogeochem. Cy., 35, e2020GB006924,, 2021. a

Falkowski, P.: Ocean Science: The power of plankton, Nature, 483, S17–S20,, 2012. a

FAO: The State of World Fisheries and Aquaculture 2020, Sustainability in action, Rome, The United Nations,, 2020. a, b

Flanagan, P., Jensen, O., Morley, J., and Pinsky, M.: Response of marine communities to local temperature changes, Ecography, 42, 214–224,, 2018. a

Geider, R., Macintyre, H., and Kana, T.: A dynamic regulatory model of phytoplanktonic acclimation to light, nutrients, and temperature, Limnol. Oceanogr., 43, 679–694,, 1998. a

Giorgetta, M., Jungclaus, J., Reick, C., Legutke, S., Bader, J., Böttinger, M., Brovkin, V., Crueger, T., Esch, M., Fieg, K., Gorges, K., Gayler, V., Haak, H., Hollweg, H.-D., Ilyina, T., Kinne, S., Kornblueh, L., Matei, D., Mauritsen, T., and Stevens, B.: Climate and carbon cycle changes from 1850 to 2100 in MPI-ESM simulations for the Coupled Model Intercomparison Project Phase 5, J. Adv. Model. Earth Syst., 5, 572–597,, 2013. a

Hashioka, T., Vogt, M., Yamanaka, Y., Le Quéré, C., Buitenhuis, E. T., Aita, M. N., Alvain, S., Bopp, L., Hirata, T., Lima, I., Sailley, S., and Doney, S. C.: Phytoplankton competition during the spring bloom in four plankton functional type models, Biogeosciences, 10, 6833–6850,, 2013. a

Heneghan, R., Galbraith, E., Blanchard, J., Harrison, C., Barrier, N., Bulman, C., Cheung, W., Coll, M., Eddy, T., Erauskin-Extramiana, M., Everett, J., Fernandes, J., Guiet, J., Maury, O., Palacios Abrantes, J., Petrik, C., Du Pontavice, H., Richardson, A., and Tittensor, D.: Disentangling diverse responses to climate change among global marine ecosystem models, Prog. Oceanogr., 198, 102659,, 2021. a

Hurrell, J. W., Holland, M. M., Gent, P. R., Ghan, S., Kay, J. E., Kushner, P. J., Lamarque, J.-F., Large, W. G., Lawrence, D., Lindsay, K., Lipscomb, W. H., Long, M. C., Mahowald, N., Marsh, D. R., Neale, R. B., Rasch, P., Vavrus, S., Vertenstein, M., Bader, D., Collins, W. D., Hack, J. J., Kiehl, J., and Marshall, S.: The Community Earth System Model: A framework for collaborative research, Bull. Am. Meteorol. Soc., 94, 1339–1360,, 2013. a

Ilyina, T., Six, K., Segschneider, J., Maier-Reimer, E., Li, H., and Núñez-Riboni, I.: Global ocean biogeochemistry model HAMOCC: Model architecture and performance as component of the MPI-Earth System Model in different CMIP5 experimental realizations, J. Adv. Model. Earth Syst., 5, 287–315,, 2013. a

Jennings, S. and Collingridge, K.: Predicting Consumer Biomass, Size-Structure, Production, Catch Potential, Responses to Fishing and Associated Uncertainties in the World's Marine Ecosystems, PLoS ONE, 10, e0133794,, 2015. a, b

Kay, J. E., Deser, C., Phillips, A., Mai, A., Hannay, C., Strand, G., Arblaster, J. M., Bates, S. C., Danabasoglu, G., Edwards, J., Holland, M., Kushner, P., Lamarque, J.-F., Lawrence, D., Lindsay, K., Middleton, A., Munoz, E., Neale, R., Oleson, K., Polvani, L., and Vertenstein, M.: The Community Earth System Model (CESM) Large Ensemble project: A community resource for studying climate change in the presence of internal climate variability, Bull. Am. Meteorol. Soc., 96, 1333–1349,, 2015. a, b

Kay, J. and Deser, C.: The Community Earth System Model (CESM) Large Ensemble Project, UCAR/NCAR Climate Data Gateway [data set], 2016. a

Kostadinov, T., Milutinovic, S., Marinov, I., and Cabre, A.: Carbon-based phytoplankton size classes retrieved via ocean color estimates of the particle size distribution, Ocean Sci., 12, 561–575,, 2016. a

Krumhardt, K., Lovenduski, N., Long, M., Levy, M., Lindsay, K., Moore, J., and Nissen, C.: Coccolithophore Growth and Calcification in an Acidified Ocean: Insights From Community Earth System Model Simulations, J. Adv. Model. Earth Syst., 11, 1418–1437,, 2019. a

Kwiatkowski, L. and Orr, J.: Diverging seasonal extremes for ocean acidification during the twenty-first century, Nat. Clim. Change, 8, 141–146,, 2018. a

Kwiatkowski, L., Torres, O., Bopp, L., Aumont, O., Chamberlain, M., Christian, J., Dunne, J., Gehlen, M., Ilyina, T., John, J., Lenton, A., Li, H., Lovenduski, N., Orr, J., Palmiéri, J., Santana-Falcón, Y., Schwinger, J., Séférian, R., Stock, C., and Ziehn, T.: Twenty-first century ocean warming, acidification, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model projections, Biogeosciences, 17, 3439–3470,, 2020. a, b, c

Lamarque, J.-F., Bond, T. C., Eyring, V., Granier, C., Heil, A., Klimont, Z., Lee, D., Liousse, C., Mieville, A., Owen, B., Schultz, M. G., Shindell, D., Smith, S. J., Stehfest, E., Van Aardenne, J., Cooper, O. R., Kainuma, M., Mahowald, N., McConnell, J. R., Naik, V., Riahi, K., and van Vuuren, D. P.: Historical (1850–2000) gridded anthropogenic and biomass burning emissions of reactive gases and aerosols: methodology and application, Atmos. Chem. Phys., 10, 7017–7039,, 2010. a, b

Lamb, S., Haacker, E., and Smidt, S.: Influence of Irrigation Drivers Using Boosted Regression Trees: Kansas High Plains, Water Resour. Res., 57, e2020WR028867,, 2021. a

Landschützer, P., Gruber, N., Bakker, D., Stemmler, I., and Six, K.: Strengthening seasonal marine CO2 variations due to increasing atmospheric CO2, Nat. Clim. Change, 8, 146–150,, 2018. a

Laufkötter, C., Vogt, M., Gruber, N., Aita-Noguchi, M., Aumont, O., Bopp, L., Buitenhuis, E., Doney, S. C., Dunne, J., Hashioka, T., Hauck, J., Hirata, T., John, J., Le Quéré, C., Lima, I. D., Nakano, H., Seferian, R., Totterdell, I., Vichi, M., and Völker, C.: Drivers and uncertainties of future global marine primary production in marine ecosystem models, Biogeosciences, 12, 6955–6984,, 2015. a, b, c

Lehodey, P., Murtugudde, R., and Senina, I.: Bridging the gap from ocean models to population dynamics of large marine predators: A model of mid-trophic functional groups, Prog. Oceanogr., 84, 69–84,, 2010. a, b

Link, J. and Marshak, A.: Characterizing and comparing marine fisheries ecosystems in the United States: determinants of success in moving toward ecosystem-based fisheries management, Rev. Fish Biol. Fish., 29, 23–70,, 2019. a

Long, M., Moore, J., Lindsay, K., Levy, M., Doney, S., Luo, J., Krumhardt, K., Letscher, R., Grover, M., and Sylvester, Z.: Simulations with the Marine Biogeochemistry Library (MARBL), J. Adv. Model. Earth Syst., 13, e2021MS002647,, 2021. a

Longhurst, A.: Ecological Geography of the Sea, Academic Press,, 2007. a, b

Lotze, H., Tittensor, D., Bryndum-Buchholz, A., Eddy, T., Cheung, W., Galbraith, E., Barange, M., Barrier, N., Bianchi, D., Blanchard, J., Bopp, L., Büchner, M., Bulman, C., Carozza, D., Christensen, V., Coll, M., Dunne, J., Fulton, E., Jennings, S., and Worm, B.: Global ensemble projections reveal trophic amplification of ocean biomass declines with climate change, P. Natl. Acad. Sci. USA, 10, 12907–12912,, 2019. a, b

Manabe, S. and Ronald, S.: Century-Scale Effects of Increased Atmospheric CO2 on the Ocean-Atmosphere System, Nature, 364, 215–218,, 1993. a

Marshak, A. and Link, J.: Primary production ultimately limits fisheries economic performance, Sci. Rep., 11, 12154,, 2021. a

Martinez-Vicente, V., Evers-King, H., Roy, S., Kostadinov, T., Tarran, G., Graff, J., Brewin, B., Dall'Olmo, G., Jackson, T., Hickman, A., Röttgers, R., Krasemann, H., Maranon, E., Platt, T., and Sathyendranath, S.: Intercomparison of Ocean Color Algorithms for Picophytoplankton Carbon in the Ocean, Front. Mar. Sci., 4, 378,, 2017. a

Maury, O.: An overview of APECOSM, a spatialized mass balanced “Apex Predators ECOSystem Model” to study physiologically structured tuna population dynamics in their ecosystem, Prog. Oceanogr., 84, 113–117,, 2010. a, b

McKinnon, K. and Deser, C.: Internal variability and regional climate trends in an observational large ensemble, J. Clim., 31, 6783–6802,, 2018. a

McKinnon, K., Poppick, A., Dunn-Sigouin, E., and Deser, C.: An “observational large ensemble” to compare observed and modeled temperature trend uncertainty due to internal variability, J. Clim., 30, 7585–7598,, 2017. a

Meehl, G., Goddard, L., Murphy, J., Ronald, S., Boer, G., Danabasoglu, G., Dixon, K., Giorgetta, M., Greene, A., Hawkins, E., Hegerl, G., Karoly, D., Keenlyside, N., Kimoto, M., Kirtman, B., Navarra, A., Pulwarty, R., Smith, D., Stammer, D., and Stockdale, T.: Decadal Prediction. Can It Be Skillful?, Bull. Am. Meteorol. Soc., 90, 1467–1485,, 2009. a, b

Meehl, G., Hu, A., Arblaster, J., Fasullo, J., and Trenberth, K.: Externally Forced and Internally Generated Decadal Climate Variability Associated with the Interdecadal Pacific Oscillation, J. Clim., 26, 7298–7310,, 2013. a

Meehl, G., Goddard, L., Boer, G., Burgman, R., Branstator, G., Cassou, C., Corti, S., Danabasoglu, G., Doblas-Reyes, F., Hawkins, E., Karspeck, A., Kimoto, M., Kumar, A., Matei, D., Mignot, J., Msadek, R., Pohlmann, H., Rienecker, M., Rosati, T., and Yeager, S.: Decadal Climate Prediction: An Update from the Trenches, Bull. Am. Meteorol. Soc., 95, 243–267,, 2014. a, b

Meinshausen, M., Smith, S., Calvin, K., Daniel, J., Kainuma, M., Lamarque, J.-F., Matsumoto, K., Montzka, S., Raper, S., Riahi, K., Thomson, A., Velders, G. J. M., and Vuuren, D.: The RCP greenhouse gas concentrations and their extensions from 1765 to 2300, Climatic Change, 109, 213–241,, 2011. a

Mills, K., Pershing, A., Brown, C., Chen, Y., Chiang, F.-S., Holland, D., Lehuta, S., Nye, J., Sun, J., Thomas, A., and Wahle, R.: Fisheries Management in a Changing Climate: Lessons From the 2012 Ocean Heat Wave in the Northwest Atlantic, Oceanography, 26, 191–195,, 2013. a

Moore, C., Morley, J., Morrison, B., Kolian, M., Horsch, E., Frolicher, T., Pinsky, M., and Griffis, R.: Estimating the Economic Impacts of Climate Change on 16 Major U.S. Fisheries, Clim. Change Econ., 12, 2150002,, 2021. a

Moore, J., Lindsay, K., Doney, S., Long, M., and Misumi, K.: Marine Ecosystem Dynamics and Biogeochemical Cycling in the Community Earth System Model [CESM1(BGC)]: Comparison of the 1990s with the 2090s under the RCP4.5 and RCP8.5 Scenarios, J. Clim., 26, 9291–9312,, 2013. a

Moore, J. K. and Braucher, O.: Sedimentary and mineral dust sources of dissolved iron to the world ocean, Biogeosciences, 5, 631–656,, 2008. a

Moore, K., Doney, S., and Lindsay, K.: Upper ocean ecosystem dynamics and iron cycling in a global three-dimensional model, Global Biogeochem. Cy., 18, GB4028,, 2004. a, b, c

Pauly, D. and Christensen, V.: Primary production required to sustain global fisheries, Nature, 374, 255–257,, 1995. a

Perry, A., Low, P., Ellis, J., and Reynolds, J.: Climate Change and Distribution Shifts in Marine Fishes, Science, 308, 1912–1915,, 2005. a

Petrik, C., Stock, C., Andersen, K., van Denderen, D., and Watson, J.: Bottom-up drivers of global patterns of demersal, forage, and pelagic fishes, Prog. Oceanogr., 176, 102124,, 2019. a

Prowe, A. F., Pahlow, M., Dutkiewicz, S., Follows, M., and Oschlies, A.: Top-down control of marine phytoplankton diversity in a global ecosystem model, Prog. Oceanogr., 101, 1–13,, 2012a. a

Prowe, A. F., Pahlow, M., and Oschlies, A.: Controls on the diversity–productivity relationship in a marine ecosystem model, Ecol. Model., 225, 167–176,, 2012b. a

Resplandy, L., Séférian, R., and Bopp, L.: Natural variability of CO2 and O2 fluxes: What can we learn from centuries-long climate models simulations?, J. Geophys. Res.-Ocean., 120, 384–404,, 2015. a, b

Roberts, D., Bahn, V., Ciuti, S., Boyce, M., Elith, J., Guillera-Arroita, G., Hauenstein, S., Lahoz-Monfort, J., Schröder, B., Thuiller, W., Warton, D., Wintle, B., Hartig, F., and Dormann, C.: Cross-validation strategies for data with temporal, spatial, hierarchical, or phylogenetic structure, Ecography, 40, 913–929,, 2016. a

Rodgers, K. B., Lee, S.-S., Rosenbloom, N., Timmermann, A., Danabasoglu, G., Deser, C., Edwards, J., Kim, J.-E., Simpson, I. R., Stein, K., Stuecker, M. F., Yamaguchi, R., Bódai, T., Chung, E.-S., Huang, L., Kim, W. M., Lamarque, J.-F., Lombardozzi, D. L., Wieder, W. R., and Yeager, S. G.: Ubiquity of human-induced changes in climate variability, Earth Syst. Dynam., 12, 1393–1411,, 2021. a

Roy, S., Sathyendranath, S., and Platt, T.: Size-partitioned phytoplankton carbon and carbon-to-chlorophyll ratio from ocean colour by an absorption-based bio-optical algorithm, Remote Sens. Environ., 194, 177–189,, 2017. a

Santer, B., Mears, C., Doutriaux, C., Caldwell, P., Gleckler, P., Wigley, T., Solomon, S., Gillett, N., Ivanova, D., Karl, T., Lanzante, J., Meehl, G., Stott, P., Taylor, K., Thorne, P., Wehner, M., and Wentz, F.: Separating signal and noise in atmospheric temperature changes: The importance of timescale, J. Geophys. Res.-Atmos., 116, D22105,, 2011. a

Sathyendranath, S., Platt, T., Kovac, Z., Dingle, J., Jackson, T., Brewin, B., Franks, P., Maranon, E., Kulk, G., and Bouman, H.: Reconciling models of primary production and photoacclimation, Appl. Optics, 59, C100–C113,, 2020. a

Schmittner, A.: Decline of the marine ecosystem caused by a reduction in the Atlantic overturning circulation, Nature, 434, 628–33,, 2005. a

Schneider, D. and Deser, C.: Tropically driven and externally forced patterns of Antarctic sea ice change: reconciling observed and modeled trends, Clim. Dynam., 50, 4599–4618,, 2018. a

Staudinger, M., Mills, K., Stamieszkin, K., Record, N., Hudak, C., Allyn, A., Diamond, T., Friedland, K., Golet, W., Henderson, M., Hernandez, C., Huntington, T., Ji, R., Johnson, C., Johnson, D., Jordaan, A., Kocik, J., Li, Y., Liebman, M., and Yakola, K.: It's about time: A synthesis of changing phenology in the Gulf of Maine ecosystem, Fish. Oceanogr., 28, 532–566,, 2019. a

Steinacher, M., Joos, F., Frölicher, T., Bopp, L., Cadule, P., Cocco, V., Doney, S., Gehlen, M., Lindsay, K., Moore, J., Schneider, B., and Segschneider, J.: Projected 21st century decrease in marine productivity: A multi-model analysis, Biogeosciences, 7, 979–1005,, 2010. a, b, c, d

Stock, C., John, J., Rykaczewski, R., Asch, R., Cheung, W., Dunne, J., Friedland, K., Lam, V., Sarmiento, J., and Watson, R.: Reconciling fisheries catch and ocean productivity, P. Natl. Acad. Sci. USA, 114, E1441–E1449,, 2017. a

Stocker, T. and Schmittner, A.: Influence of CO2 emission rates on the stability of the thermohaline circulation, Nature, 388, 862–865,, 1997. a

Tagliabue, A., Kwiatkowski, L., Bopp, L., Butenschön, M., Cheung, W., Lengaigne, M., and Vialard, J.: Persistent Uncertainties in Ocean Net Primary Production Climate Change Projections at Regional Scales Raise Challenges for Assessing Impacts on Ecosystem Services, Front. Clim., 3, 738224,, 2021. a, b, c, d, e

Timmermann, A., Oberhuber, J., Bacher, A., Esch, M., Latif, M., and Roeckner, E.: Increased El Niño frequency in a climate model forced by future greenhouse warming, Nature, 398, 694–697,, 1999. a

Tittensor, D., Eddy, T., Lotze, H., Galbraith, E., Cheung, W., Barange, M., Blanchard, J., Bopp, L., Bryndum-Buchholz, A., Büchner, M., Bulman, C., Carozza, D., Christensen, V., Coll, M., Dunne, J., Fernandes, J., Fulton, E., Hobday, A., Huber, V., and Walker, N.: A protocol for the intercomparison of marine fishery and ecosystem models: Fish-MIP v1.0, Geosci. Model Dev., 11, 1421–1442,, 2018. a, b

Tittensor, D., Blanchard, J., Fulton, E., Cheung, W., Novaglio, C., Harrison, C., Heneghan, R., Barrier, N., Bianchi, D., Bopp, L., Bryndum-Buchholz, A., Britten, G., Büchner, M., Christensen, V., Coll, M., Dunne, J., Eddy, T., Everett, J., Fernandes, J., and Stock, C.: Next-generation ensemble projections reveal higher climate risks for marine ecosystems, Nat. Clim. Change, 11, 973–981,, 2021. a, b

Travers-Trolet, M., Shin, Y.-J., Jennings, S., Machu, E., Huggett, J., Field, J., and Cury, P.: Two-way coupling versus one-way forcing of plankton and fish models to predict ecosystem changes in the Benguela, Ecol. Model., 220, 3089–3099,, 2009.  a, b

Vichi, M., Allen, I., Masina, S., and Hardman-Mountford, N.: The emergence of ocean biogeochemical provinces: A quantitative assessment and a diagnostic for model evaluation, Global Biogeochem. Cy., 25, GB2005,, 2011. a, b, c, d, e, f, g

Wernberg, T., Bennett, S., Babcock, R., de Bettignies, T., Cure, K., Depczynski, M., Dufois, F., Fromont, J., Fulton, C., Hovey, R., Harvey, E., Holmes, T., Kendrick, G., Radford, B., Santana-Garcon, J., Saunders, B., Smale, D., Thomsen, M., Tuckett, C., and Wilson, S.: Climate-driven regime shift of a temperate marine ecosystem, Science, 353, 169–172,, 2016. a

Yamaguchi, R., Rodgers, K., Timmermann, A., Stein, K., Schlunegger, S., Bianchi, D., Dunne, J., and Slater, R.: Trophic level decoupling drives future changes in phytoplankton bloom phenology, Nat. Clim. Change, 12, 1–8,, 2022. a, b

Short summary
Anthropogenic climate change will influence marine phytoplankton over the coming century. Here, we quantify the influence of anthropogenic climate change on marine phytoplankton internal variability using an Earth system model ensemble and identify a decline in global phytoplankton biomass variance with warming. Our results suggest that climate mitigation efforts that account for marine phytoplankton changes should also consider changes in phytoplankton variance driven by anthropogenic warming.
Final-revised paper