Cascading effects augment the direct impact of CO₂ on phytoplankton growth in a biogeochemical model

Marine primary production, mostly performed by unicellular algae (Middelburg, 2019), contributes to the global carbon cycle by the fixation of CO₂ in organic biomass via photosynthesis in the surface ocean and subsequent sinking and remineralization of biomass at depth (organic carbon pump). By the process of calcification, on the other hand, CO₂ is released into the surrounding water (carbonate pump). About half of the pelagic calcification is performed by the phytoplankton group of coccolithophores (Brownlee and Taylor, 2004; Beardall and Raven, 2013), which therefore contributes to both the organic carbon and the carbonate pump. Although their contribution to global primary production is relatively small (<10%; Poulton et al., 2007), some coccolithophore species can form large blooms that are observable from satellites (Hopkins and Balch, 2018) and autonomous ocean profilers (Terrats et al., 2020). Locally, calcification can thereby considerably dampen the net oceanic uptake of atmospheric CO₂ (Shutler et al., 2013).

Due to increasing anthropogenic CO₂ emissions, environmental conditions in the oceans are currently changing at an unprecedented rate, resulting in oceans that are warmer, more acidic, less oxygenated, and poorer in surface-ocean nutrients (Gruber, 2011; Kwiatkowski et al., 2020). In response to these changes, a decrease in global marine primary production (Moore et al., 2018) and shifts in species composition (e.g., Barton et al., 2016; Poloczanska et al., 2016) are expected in the future. Effects on the organic carbon and carbonate pump can already be observed locally (e.g., Bindoff et al., 2019; Pinkerton et al., 2021). For instance, observations in the North Atlantic suggest up to a 10-fold increase in coccolithophore abundance over the last decades (1965-2010, Rivero-Calle et al., 2015; 1990-2012, Krumhardt et al., 2016). Yet, the underlying reasons for this increase are not understood conclusively, with different studies suggesting either warming (Beaugrand et al., 2013) or a positive effect of dissolved CO₂ (CO_2(aq); Rivero-Calle et al., 2015; Krumhardt et al., 2016) as the primary cause for higher coccolithophore abundances.

Temperature, nutrient concentrations, and light availability have long been considered as the most important environmental drivers that control the magnitude of phytoplankton growth and coccolithophore calcification. However, even the carbonate system itself influences growth and calcification and becomes increasingly relevant given the current pace of ocean acidification. The impact of the carbonate system can be described with an optimum function: CO_2(aq) and bicarbonate (⁠$HCO3−$⁠) are the primary substrates of growth and calcification (e.g., Kottmeier et al., 2014) and rising substrate levels (carbonation) can enhance growth and calcification. Rising CO_2(aq) concentrations result, however, in a concomitant increase in proton concentrations and thus a lower pH (acidification), which in turn eventually dampens growth and calcification (e.g., Rost and Riebesell, 2004; Gao et al., 2019). Phytoplankton species and groups differ considerably in their ability to make use of carbonation and in their sensitivity to acidification (e.g., Rost et al., 2008; Dutkiewicz et al., 2015). For instance, diatoms are affected only moderately by changes in CO_2(aq) and $HCO3−$ concentrations due to their exceptionally effective CO₂ concentrating mechanisms (CCMs; Reinfelder, 2011; Pierella Karlusich et al., 2021). By diversely changing the competitive fitness of species and groups, CO₂ has also an indirect effect on phytoplankton community structure (e.g., Doney et al., 2009; Paul and Bach, 2020). Calcification, however, is particularly sensitive to ocean acidification, and laboratory experiments show that coccolithophore growth is also impacted disproportionately by increasing CO_2(aq) levels (Hoppe et al., 2011; Meyer and Riebesell, 2015; Zhang et al., 2019; Seifert et al., 2020). The opposing effects of carbonation and acidification lead to a transition of beneficial to detrimental impacts with increasing partial pressure of CO₂ (pCO₂), which has been described mechanistically by several authors (Bach et al., 2011; Bach et al., 2015; Gafar et al., 2018; Paul and Bach, 2020).

Marine biogeochemical models are important tools with predictive capabilities to understand connections and variations of biological, chemical, and biogeochemical components in space and time (e.g., Holt et al., 2014). In these models, temperature, light, and nutrients are traditionally considered as the environmental drivers to control phytoplankton growth, and very few models to date account for the effects of CO₂ on growth and calcification. Dutkiewicz et al. (2015) assessed the pCO₂ impact on all phytoplankton functional groups in the ecosystem model DARWIN from preindustrial to 2100 levels under the high emission scenario (RCP8.5) using different conceptual approaches based on linear, Michaelis-Menten type, and Hill type functions. They found that the picocyanobacteria Synechococcus and diazotrophs benefit most from increasing pCO₂, while the biomass of other groups, including coccolithophores and diatoms, tend to decrease. In Krumhardt et al. (2019), only coccolithophore growth was considered to be limited by low CO_2(aq) concentrations using a Monod function, similar to other nutrient limitations. The ratio of particulate inorganic (PIC) to organic carbon (POC) production by coccolithophores (PIC:POC)_cocco and, thus, the calcite production for a given coccolithophore biomass, was defined to decrease linearly with increasing CO_2(aq). Globally, model runs yielded a decrease in coccolithophore calcification by 17%, but an increase in coccolithophore growth at 900 μatm atmospheric CO₂ relative to present-day CO₂ levels. Gangstø et al. (2011) did not account for CO₂ effects on growth, but described a Michaelis-Menten-like dependence of the PIC:POC ratio on the calcite saturation state. In their model, implicit total calcification decreases by 20–60% in a high-emission scenario by the year 2100. Thus, while the CO₂ impacts on coccolithophores are not consistent across previous studies, models generally agree that increasing CO₂ can considerably alter the community structure and lead to a decrease in calcification.

None of these models, however, explicitly considers the impact of carbon speciation and the resulting two-sided carbonation-acidification effect on both growth and calcification of phytoplankton. Further, CO₂ effects on phytoplankton growth and calcification are currently largely omitted in the models of the Coupled Model Intercomparison Project (CMIP) Phase 6 (Eyring et al., 2016), which are used for the Sixth Assessment Report (AR6) of the IPCC. Besides, due to their considerable global contribution to the carbonate pump and the particular sensitivity of calcification to ocean acidification, a thorough assessment of CO₂ effects requires the implementation of coccolithophores as a separate phytoplankton group. So far, they have only been included in a few global and regional models, mainly representing the bloom-forming species Emiliania huxleyi (Le Quéré et al., 2005; Follows et al., 2007; Gregg and Casey, 2007; Dutkiewicz et al., 2015; Kvale et al., 2015; Nissen et al., 2018; Krumhardt et al., 2019), whereas most biogeochemical models do not represent coccolithophores or calcification explicitly.

In the present study, we aim to quantify the responses of marine phytoplankton biomass and coccolithophore calcification to different atmospheric CO₂ levels and distinguish between the direct and the indirect effects of CO₂ using the physical-biogeochemical model FESOM-REcoM. First, we describe the implementation of coccolithophores and their calcification into the model. Based on a literature review of laboratory experiments, we then developed optimum functions describing the sensitivity to carbonation and acidification of growth of all modelled phytoplankton groups (coccolithophores, diatoms, small phytoplankton) and for calcification. These functions were subsequently applied in idealized simulations at preindustrial, present-day, and future atmospheric CO₂ levels. Finally, motivated by the observational studies of Beaugrand et al. (2013), Rivero-Calle et al. (2015), and Krumhardt et al. (2016), we have attempted to disentangle the effects of warming and ocean acidification on changes in global and North Atlantic coccolithophore biomass during the last decades.

We used the global Regulated Ecosystem Model version 2 (REcoM-2-M) coupled to the Finite Element Sea Ice-Ocean Model (FESOM 1.4, Schourup-Kristensen et al., 2014; Wang et al., 2014). The ocean model FESOM applies a finite element method that solves primitive hydrostatic equations on an unstructured mesh, allowing for flexible multi-resolution meshes (Sidorenko et al., 2011; Wang et al., 2014). REcoM-2-M describes the biogeochemical cycling of carbon, nitrogen, silicon, iron, and oxygen (Hauck et al., 2013; Karakuş et al., 2021). The CO₂ flux from the atmosphere to the ocean and the 3D carbonate system are computed by the mocsy 2.0 routine (Orr and Epitalon, 2015). In the previous version of the model (Karakuş et al., 2021), the ecosystem is described by two phytoplankton groups (diatoms and small-sized phytoplankton), two zooplankton groups (small, fast-growing zooplankton and slow-growing polar macrozooplankton), and two detritus groups (small, slow-sinking and large, fast-sinking particles). While the classification of diatoms is taxonomic, the small phytoplankton comprises a wide range of taxa, such as non-silicifying and non-calcifying haptophytes and green algae. In this study, we have extended the phytoplankton groups by coccolithophores, which we consider as a taxonomically separate group, similar to diatoms.

Changes in biomass result from growth and loss terms. We describe environmental impacts on growth of all phytoplankton groups, with a special focus on the parameter choice for coccolithophores. Parameters were chosen to lie within ranges based on qualitative literature reviews, and further tuned to reproduce observed coccolithophore biomass distributions (MAREDAT dataset; O’Brien et al., 2013). While the bloom-forming coccolithophore species E. huxleyi is numerically the most important coccolithophore species (Paasche, 2001) and most frequently used in laboratory studies, it is not a typical representative of coccolithophores (Rost and Riebesell, 2004). Therefore, we aimed to represent the average characteristics of the diverse group of coccolithophores in our model in the parameter tuning. Subsequently, we describe the explicit implementation of coccolithophore calcification and CO₂-dependent calcite dissolution.

The growth rate GR of phytoplankton group p (coccolithophores, diatoms, or small phytoplankton) depends on the group-specific constant maximum growth rate $μmax$ (Table 1), nutrient and light limitation terms (⁠$fN$ and $fL$⁠), and a temperature function (⁠$fT$⁠):

According to Margalef’s mandala (Margalef, 1978; Wyatt, 2014), diatoms usually follow the productive r-strategy with higher growth rates, and coccolithophores the efficient K-strategy with lower growth rates. While our group of small phytoplankton includes various species that presumably cover a wide range of maximum growth rates, we assume here that the majority of coccolithophore species, other than the bloom-forming E. huxleyi, are growing rather slowly and are specialized to low nutrient regions and seasons. Thus, we chose a smaller $μmax$ for coccolithophores than for diatoms and small phytoplankton in the model (Table 1).

Nutrient limitation f^N in Equation 1 is determined by the most limiting nutrient (dissolved inorganic nitrogen, DIN; dissolved iron, DFe; and, for diatoms, also dissolved silicic acid, DSi), and is defined as 0 ≤ $fN$ ≤ 1 with 0 denoting strongest limitation and 1 denoting no limitation. Limitations by DIN and DSi (⁠$lDIN$ and $lDSi$⁠) depend on the intracellular nitrogen or silicate to carbon quotas, and the group-specific half-saturation constants $kDIN$ and $kDSi$ (Table 1) for DIN and DSi, respectively (Geider et al., 1998; Hauck et al., 2013). Limitation by DFe (⁠$lDFe$⁠) is described by a Michaelis-Menten function, which is dependent on the group-specific half-saturation constant $kDFe$ (Table 1) and the concentration of DFe in the watercolumn (i.e., $[DFe]/([DFe]+kDFe)$⁠). For coccolithophores and small phytoplankton, nutrient limitation is then defined as:

Diatoms are additionally affected by DSi limitation (⁠$lDSi$⁠). They have a high uptake rate at high nutrient concentrations (e.g., luxury consumption) and are less efficient under low nutrient concentrations (Sarthou et al., 2005; Litchman et al., 2007; Marañón et al., 2013). Small phytoplankton, by contrast, have the greatest ability to thrive under low nutrient conditions, e.g. due to their small cell size (Irwin et al., 2006; Marañón et al., 2013; Irion et al., 2021). Half-saturation constants $kDIN$ and $kDFe$ of coccolithophores were chosen to be smaller than those of diatoms and higher than those of small phytoplankton because their luxury consumption for nitrate is low compared to diatoms, and they are less competitive under low nutrient concentrations than other phytoplankton groups (Riegman et al., 2000; Rost and Riebesell, 2004). Literature values for maximum and minimum nitrate-to-carbon (N:C) ratios (⁠$qmax$ and $qmin$⁠) do not show a robust difference between phytoplankton groups (Finkel et al., 2009). For diatoms and small phytoplankton, $qmax$ and $qmin$ (Table 1) translate into C:N ratios between 5 and 20, which are well in line with values reported from laboratory studies of diatoms (C:N range: 2.7–29.7; Sarthou et al., 2005). For coccolithophores, lower values result from model tuning, and $qmax$ and $qmin$ translate into C:N ratios between 6.6 and 25. For the temperature sensitivity $fT$ in Equation 1, we used a power function proposed by Fielding (2013) for coccolithophores (⁠$fcoccoT$⁠), which is based on the experimental growth rate-temperature relations of E. huxleyi (Figure S1):

with $Tdeg$ being the temperature ≥0°C. At temperatures <0°C, the function was set to a small value (2.23 × 10⁻¹⁶) rather than zero for numerical reasons. For diatoms and small phytoplankton, $fT$ is defined by an Arrhenius function (Figure S1):

with T_K being the temperature in K, and $TK,ref$ being the reference temperature of 288.15 K (15°C). Light limitation $fL$ in Equation 1 is calculated as a function of the group-specific maximum light harvesting efficiency α (Table 1), which represents the initial slope of the photosynthesis-irradiance curve, $μmax$ (Table 1), photosynthetically active radiation (PAR), the variable chlorophyll-to-carbon ratio $qChl:C$⁠, as well as nutrient limitation $fN$ (Equation 2) and the temperature function $fT$ (Equations 3 and 4; Geider et al., 1998):

Consequently, changes in either the chlorophyll synthesis or the biomass production and the resulting alteration of $qChl:C$ affect the light limitation term. While coccolithophores are tolerant of high light levels, their ability to grow under low light conditions is lower compared to diatoms (Nielsen, 1997; Paasche, 2001; Le Quéré et al., 2005; Zondervan, 2007). In our model, coccolithophores therefore have a small maximum light harvesting efficiency α and a lower $qmaxChl:N$ than diatoms (Table 1).

Loss terms of phytoplankton biomass are composed of grazing, aggregation, respiration, and excretion of dissolved organic matter. Grazing is parametrized with the Fasham formulation (Fasham et al., 1990; Vallina et al., 2014) and a variable grazing preference depending on a group-specific initial preference term and the relative contribution of each phytoplankton group to total phytoplankton biomass (Karakuş et al., 2021). We assumed that the calcite shell of coccolithophores does not serve as grazing protection (Mayers et al., 2020) and therefore chose a similar effective zooplankton grazing on coccolithophores as on small phytoplankton. We acknowledge that a protective function of the calcite shell is still under debate (e.g., Monteiro et al., 2016) and that a coccolithophore diet might reduce zooplankton growth rates by reduced digestion rates or increased swimming efforts due to the mass of the indigestible coccoliths (Haunost et al., 2021). Because the effective grazing in the model is dependent on the biomass of the respective phytoplankton group and because coccolithophore biomass is an order of magnitude smaller than small-phytoplankton biomass (see Section Representation of coccolithophores and biogeochemical fluxes in FESOM-REcoM), we used a higher grazing preference of the small zooplankton group on coccolithophores than on small phytoplankton (Table 1) to obtain a similar grazing loss of both groups. We chose the same grazing preference of the polar macrozooplankton group on coccolithophores as on small phytoplankton as its biomass is generally highest in areas with low coccolithophore biomass.

Deviating from the model version described in Hauck et al. (2013) in which calcification is performed by a constant fraction of the small phytoplankton with a fixed $(PIC:POC)sphy$ ratio of 1, our model calculates calcification (⁠$Calccocco$⁠) depending on the specific growth rate of coccolithophores (⁠$GRcocco$⁠; Equation 1), coccolithophore biomass (C_cocco), and a reference PIC:POC ratio, $(PIC:POC)coccoref$⁠, that is modified by temperature (⁠$fCaCO3T$⁠) and DIN limitation (⁠$fCaCO3N$⁠):

Here, we use a $(PIC:POC)coccoref$ ratio of 1, which is motivated by a compilation of PIC:POC ratios of different coccolithophore species showing that PIC:POC ratios >1 can occur where highly calcified species like Calcidiscus leptoporus dominate, and PIC:POC ratios <1 can occur in the Southern Ocean if the lightly calcified E. huxleyi morphotype B/C dominates (Krumhardt et al., 2017). Temperature dependence of calcification $fCaCO3T$ with decreasing (PIC:POC)_cocco ratio at temperatures <10.6°C follows Krumhardt et al. (2017):

with $Tdeg$ being the temperature in °C. Krumhardt et al. (2017) describe the modification of the (PIC:POC)_cocco ratio by nutrient availability using an equation of the following form:

In their study, [N] is the PO₄ concentration and $kcoccoN$ is the coccolithophore half-saturation constant for phosphate uptake. In their literature review, Krumhardt et al. (2017) found a 37% increase in (PIC:POC)_cocco from phosphate-repleted to phosphate-limited condition, from which they derived x = –0.48 and y = 1.48 (both unitless). Because our model does not include phosphate, we adapted Equation 8 to (PIC:POC)_cocco modification under DIN limitation, with [N] and $kcoccoN$ being the DIN concentration and coccolithophore half-saturation constant for nitrate uptake (Table 1), respectively. Krumhardt et al. (2017) also assessed the effect of nitrate limitation on (PIC:POC)_cocco, with a 25% higher ratio under nitrate-limited compared to nitrate-repleted conditions from which we derived x = –0.31 and y = 1.31.

Instead of an exponential increase with depth with a vertical length scale of 3500 m following Yamanaka and Tajika (1996) as in the previous model version (Hauck et al., 2013), calcite dissolution λ now depends on the carbonate ion concentration following Aumont et al. (2015):

with Ca²⁺ being the calcium ion concentration, kspc the stoichiometric solubility product (Mucci, 1983), $CO32−$ the carbonate ion concentration, and $CO3 sat2−$ the saturated carbonate ion concentration. The carbonate system was computed once per week where PAR>1% of surface PAR, and once per month where PAR≤1% of surface PAR to increase computational efficiency. To account for calcium carbonate dissolution above the carbonate saturation horizon, we introduced an additional dissolution in zooplankton guts. The acidic environment in guts of starving copepods can dissolve up to 38% of the calcite taken up by grazing (Pond et al., 1995; Jansen and Wolf-Gladrow, 2001; White et al., 2018). In addition, aragonite and high-magnesium forms of calcium carbonate have a shallower saturation horizon than calcite and contribute to upper-ocean calcium carbonate dissolution (Sabine et al., 2002; Feely et al., 2004; Barrett et al., 2014; Battaglia et al., 2016; Sulpis et al., 2021). Besides, high-CO₂ microenvironments in aggregates caused by metabolic activity of microbes can lead to calcium carbonate dissolution even above the saturation horizon (Milliman et al., 1999; Sulpis et al., 2021). Hence, we have deliberately overestimated gut dissolution with 50% of ingested calcite to account for aragonite and high-magnesium calcite dissolution that is not covered by the calcite-based dissolution formulation of Aumont et al. (2015).

For the implementation of CO₂ dependencies, we added limitation functions $fCO2$ similar to the dependencies of other environmental drivers to Equations 1 and 6:

Multiple concepts describing growth and calcification dependencies on the carbonate system exist in the literature (Bach et al., 2011; Bach, 2015; Bach et al., 2015; Gafar et al., 2018; Paul and Bach, 2020). Here we used the concept of Bach et al. (2015), which has been developed for coccolithophore calcification, as it can be well accommodated to our model structure. Furthermore, in contrast to the concepts of Paul and Bach (2020, modified Monod equation based on a generic proton to CO_2(aq) relationship) and Bach et al. (2011, modified Michaelis-Menten function depending only on fCO_2(aq)), it allows non-linear changes between carbon species concentrations. Because they are less well suited for our model structure, we refrained from using the concepts of Gafar et al. (2018, CO₂ term including temperature and light effects and interactions therein) and Bach (2015, bicarbonate to proton ratio). The concept of Bach et al. (2015) accounts for changes in the bicarbonate (⁠$HCO3−$⁠), CO_2(aq), and proton concentrations in a modified Michaelis-Menten function:

with p indicating the CO₂ effect on phytoplankton growth (diatoms, small phytoplankton, or coccolithophores), and CaCO₃ indicating the CO₂ effect on coccolithophore calcification. Although the effect of CO₂ on phytoplankton growth and coccolithophore calcification can be either limiting (on the carbonation side) or inhibiting (on the acidification side) depending on the state of the carbonate system, we use the term “CO₂ limitation” for both in the remaining text equivalent to nutrient and light limitation.

To derive parameters in Equation 13 we used published laboratory experiments describing phytoplankton growth rates and PIC:POC ratios of coccolithophores at different pCO₂ levels. The collection builds on the literature review by Paul and Bach (2020), complemented by experimental data of other publications that fit our selection criteria (Table 2; for more details see Tables S1–S4). In particular, the data had to i) belong to one of our modelled phytoplankton groups, ii) report the response variable (growth rate and/or (PIC:POC)_cocco ratio) for at least three pCO₂ treatment levels, and iii) indicate the concentrations of all relevant carbonate system parameters (CO_2(aq), $HCO3−$⁠, and proton concentrations or pH). By including as many species, ecotypes, and strains as possible, we acknowledge intraspecific plasticity, which can alleviate the sensitivities towards changing environmental conditions (e.g., Wolf et al., 2018) and, thus, lead to more conservative estimates of CO₂-driven changes in phytoplankton growth and coccolithophore calcification.

We used the R function nls to fit Equation 13 with a non-linear least square estimate to the laboratory data of group-specific phytoplankton growth rates and (PIC:POC)_cocco ratios. We forced the functions through the origin to constrain the lower end of the curve to reasonable values (i.e., preventing the fit to yield positive growth or calcification when concentrations of inorganic carbon are zero). Prior to that, experimental growth rates and (PIC:POC)_cocco ratios were normalized, with the highest value within each dataset being scaled to 1. Because these values are reached over a wide range of carbonate species concentrations, fitted functions yielding a maximum value of 1 is not possible mathematically. To obtain functions that are non-limiting (i.e., reach a value of 1) at ideal conditions, we forced the functions through the median that was obtained from the inorganic carbon species concentrations at the maximum growth rate or (PIC:POC)_cocco ratio of each dataset. To test whether the fitted functions would be biased by extreme pCO₂ levels tested in some experiments (up to 5500 μatm), we also fitted the function only through data points below 2000 and 1000 μatm, respectively (Figure S2). In both cases, growth and (PIC:POC)_cocco were more sensitive to increasing pCO₂, showing that including extreme pCO₂ levels does not exacerbate but softens the impact of CO₂ on growth and calcification. Hence, we selected the functions that were fitted over the entire pCO₂ range. Parameter values resulting from the function fits are listed in Table 2. Because the fits were deducted from experimental carbonate system environments, extreme carbonate system environments in the model, e.g. in the Arctic Ocean with freshwater influence, resulted in carbonate system sensitivities >1. For the implementation in REcoM-2-M, we limited the resulting term to be ≤1 to prevent the carbonate system from enhancing growth and calcification to more than the pre-set maximum value.

Figure 1 displays the function fits for a sample matrix of bicarbonate, CO_2(aq), and proton concentrations. Calcification and diatom growth benefit strongest from increasing pCO₂ at very low levels due to the carbonation effect (Figure 1). The steep increase in diatom growth rate already at low pCO₂ is well in line with the presence of effective and highly regulated CCMs in the group of diatoms (Reinfelder, 2011; Pierella Karlusich et al., 2021). The relatively stable diatom growth rate over a wide pCO₂ range after the steep initial increase displayed in our function can mirror the reallocation of energy, which is saved due to downregulation of CCMs under higher CO_2(aq) and bicarbonate availability (Shi et al., 2019). The increment from preindustrial (280 μatm) to present-day (420 μatm) pCO₂ levels increases the growth rate of small phytoplankton, whereas the growth rate of the other phytoplankton groups and the calcification rate decrease in this pCO₂ range. Besides, growth rates of small phytoplankton depict a relatively broad plateau close to maximum growth rates between present-day and future (750 μatm) pCO₂ levels (Figure 1). Studies show that small phytoplankton species such as Micromonas spp. and single trichomes of Trichodesmium spp. benefit from pCO₂ levels higher than preindustrial (Boatman et al., 2018), whereas pCO₂ levels higher than present-day must not necessarily trigger further stimulation in growth (Maat et al., 2014; Boatman et al., 2018; White et al., 2020). Calcification decreases most steeply after reaching an optimum at relatively low pCO₂ in our function (Figure 1), agreeing with experimental studies that reported coccolithophore calcification unambiguously to be the process most sensitive to any variability in the carbonate system, with an overall strong decrease in (PIC:POC)_cocco with increasing pCO₂ (e.g., Findlay et al., 2011; Meyer and Riebesell, 2015). In natural communities, this effect is caused by not only a decrease in calcification per cell, but also a shift from highly to less calcified coccolithophore species and morphotypes (Beaufort et al., 2011). Coccolithophore growth rates in our function initially decrease more steeply with increasing pCO₂ than small-phytoplankton growth rates, and are only outpaced by small-phytoplankton growth rates at high pCO₂ levels (Figure 1). Our fits are in line with experimental studies which found coccolithophore growth to decrease most strongly with increasing pCO₂ (Riebesell et al., 2017; Seifert et al., 2020), with the exception of some coccolithophore strains (Meyer and Riebesell, 2015). In summary, we are confident that our fitted functions based on the concept of Bach et al. (2015) and constrained with laboratory data simulate the general response patterns of phytoplankton groups to different pCO₂ levels reasonably well.

We performed two sets of model simulations on the CORE-II mesh with higher resolved dynamic areas (up to 20 km, e.g. coastlines and equatorial belt) and coarser resolution of less dynamic areas (up to 150 km, open ocean; Sidorenko et al., 2011; Wang et al., 2014). We used the JRA-55-do reanalysis data set as atmospheric forcing in all simulations (version 1.3.1; Tsujino et al., 2018). Temperature and salinity were initialized with hydrographic data from the Polar Science Center Hydrographic Climatology (PHC 3.0; updated from Steele et al., 2001), DIN and DSi from the World Ocean Atlas (Garcia et al., 2013), and DFe from PISCES output (Aumont et al., 2003) which was corrected using observed profiles (de Baar et al., 1999; Boye et al., 2001). DIC and alkalinity were initialized from GLODAPv2 (Lauvset et al., 2016), and biomass from small values (2.23 × 10⁻⁵ for biomass in nitrogen units and 14.77 × 10⁻⁵ for biomass in carbon units, i.e. converted with the Redfield ratio).

In the first set of simulations, we aimed to investigate the effect of different atmospheric CO₂ concentrations on phytoplankton growth and biomass. We performed in total four model simulations without and with CO₂ dependence at present-day atmospheric CO₂ levels (PRESENT and PRESENT_CO2, respectively, 420 μatm; Table 3) and with CO₂ dependence under preindustrial and future CO₂ levels (PREIND_CO2, 280 μatm, and FUTURE_CO2, 750 μatm; Table 3). The simulations were forced with repeated year forcing, looping over the year 1961. We initialized DIC with preindustrial DIC for the PREIND_CO2 simulation, and with present-day DIC in all others (Lauvset et al., 2016). Each simulation was run for 32 years, and a mean over the last five years was analyzed for this study. Two additional simulations were performed to test the sensitivity of phytoplankton biomass towards the parameter choice of the CO₂ dependence by increasing and decreasing each parameter by 10% (Figure S3).

In a second set of four simulations, we aimed to disentangle the effects of warming and CO₂ on phytoplankton biomass over the period 1958 to 2019. In one simulation, both climatological and atmospheric CO₂ forcing varied inter-annually (VARCLI_VARCO2; Table 3) which accounts for both the radiative and the geochemical climate change effect. Further, a second simulation was forced by repeating the year 1961 for all atmospheric variables except CO₂ (CONSTCLI_VARCO2; Table 3), while in a third, atmospheric CO₂ was held constant at 280 μatm and all other variables varied inter-annually (VARCLI_CONSTCO2; Table 3). With these two simulations, we can separate the effects of the atmospheric CO₂ increase and other changing atmospheric forcing fields (air temperature, radiation, humidity, freshwater fluxes, wind). In a control simulation (CTRL; Table 3), both repeated year forcing of the year 1961 and constant atmospheric CO₂ were used. All simulations of the second set were started from previously spun-up model states, where the spin-ups covered the period 1850 to 1957 using JRA repeated year forcing of the year 1961 (Tsujino et al., 2018). The spin-up simulation used to initialize the VARCO2 simulations was forced with increasing atmospheric CO₂ concentration as in Friedlingstein et al. (2020), and the spin-up simulation used to initialize the CONSTCO2 simulations was forced with constant atmospheric CO₂ of 278 ppm. As these spin-ups were conducted with the prior model version without coccolithophores and CO₂ dependence, we reset all biological tracers (phytoplankton and zooplankton biomass) to small initialization values, while taking advantage of the spun-up state for physical fields, carbonate chemistry and nutrients. This approach was chosen to ensure that the simulated carbonate chemistry represents the historical evolution of atmospheric CO₂. We started our analysis in year 4 of the simulations when the upper ocean ecosystem reached a quasi-equilibrium.

Simulated coccolithophore biomass concentrations of the PRESENT and the PRESENT_CO2 simulations were evaluated with a compilation of observations of global coccolithophore biomass concentrations (MARine Ecosystem DATa dataset, MAREDAT; O’Brien et al., 2013) and literature estimates of global phytoplankton groups biomass, net primary production, and biogeochemical fluxes. Additionally, we evaluated the diatom representation in our model with the corresponding MAREDAT observational dataset (Leblanc et al., 2012). For a quantitative comparison of the model with the MAREDAT datasets, the modelled coccolithophore and diatom biomass was subsampled only for grid points with a corresponding observational grid point and, thus, does not represent the entire modelled biomass of each group. Subsequently, biomass estimates of the model and the MAREDAT datasets were separated into oceanic biomes (Fay and McKinley, 2014) to assess the regional representation of simulated coccolithophores and diatoms compared to observations. An additional analysis of ecological niche occupations in the model compared to a statistical model of observations is described in Text S1.

To evaluate the phytoplankton representation in our model, we focus on the PRESENT simulation, which does not include the CO₂ dependence. To assess the impact of the CO₂ dependence at present-day atmospheric CO₂ levels, we then compare the PRESENT with the PRESENT_CO2 simulation. Simulated annual mean coccolithophore biomass is highest in the North Atlantic, the equatorial regions, upwelling regions, and the northern boundary of the Southern Ocean (Figure 2). It thereby resembles the distribution of diatoms except in the polar regions, where diatom biomass is highest of all three phytoplankton groups, whereas coccolithophore biomass is very small in high latitudes. Small phytoplankton biomass is more homogeneously distributed, with slightly lower biomass in the subtropical gyres and the polar regions compared to the remaining ocean. Global coccolithophore biomass is about one order of magnitude smaller than that of diatom or small phytoplankton (Figure 2).

The comparison to the MAREDAT dataset (O’Brien et al., 2013) reveals that our model overestimates coccolithophore biomass in the northern hemisphere and the equatorial region, especially in the subtropical seasonally stratified (ST-SS) north region that covers major parts of our coccolithophore patch in the North Atlantic (0.3 × 10⁻³ Tg C in the observations versus 5 × 10⁻³ Tg C in the model; Figure 3A). Modelled coccolithophore biomass in the ST-SS (south) region, covering the coccolithophore patch at the northern border of the Southern Ocean, is fairly close to observational data with about 2 × 10⁻³ Tg C. Because the model does not depict the coccolithophore biomass of 0.1 × 10⁻³ Tg C from observations further south in the subpolar seasonally stratified (SP-SS) region (Figure 3A), high coccolithophore biomass usually occurring between 40 and 60°S in the Southern Ocean (Balch et al., 2011) is shifted northwards in our model. The coccolithophore niche analysis (see Text S1 and Figure S4) reveals a good model representation of E. huxleyi and coccolithophore species in warm regions, while coccolithophore species in cold regions are not represented. Diatom biomass is in the range of observational data for most regions, except for an overestimation of modelled biomass in the ST-SS (north) region and a strong underestimation in the equatorial region (Figure 3B; Leblanc et al., 2012). We attribute this underestimation to a patch with very high observed diatom biomass in the equatorial Atlantic which is not captured by the model, as well as a mismatch between high diatom biomass on the model grid and the MAREDAT grid points for which the model was sampled. For both coccolithophores and diatoms, trends in inter-biome biomass differences are similar in the grid subsampled for MAREDAT datapoints (Figure 3) and the entire model grid (Figure 2; data not shown as barplots).

Globally integrated annual mean coccolithophore biomass amounts to 0.03 Pg C, corresponding to the upper end of literature estimates (0.001–0.032 Pg C), and to 2.3% of the total simulated phytoplankton biomass (Table 4). The contribution of coccolithophores to total NPP is slightly higher in our model (3.3% of total NPP) than in literature estimates (2% of total NPP; Table 4). Diatom biomass (0.49 Pg C), silicate export (84.3 Tmol yr⁻¹), and the relative share of diatoms to total NPP (29.9%) fit well into the fairly wide range of literature estimates (0.01–0.94 Pg C, 69–185 Tmol yr⁻¹, and 15–35%, respectively; Table 4). With approximately 30 Pg C yr⁻¹, total NPP in our model is at the lower end of literature estimates (Table 4), which can be attributed at least partly to an underestimation of small-phytoplankton productivity, given that both coccolithophore and diatom NPP are likely at the upper end, and to the absence of nitrogen fixers in our model that usually increase primary productivity in the low latitudes of the real ocean. The (PIC:POC)_cocco of 1.2 and the calcite export of 1.2 Pg C yr⁻¹ are well within the range of literature values (0.7–1.4 and 0.1–4.7 Pg C yr⁻¹, respectively; Table 4), indicating a good representation of coccolithophore calcification and calcite dissolution in the model. With 0.3, the rain ratio (PIC:POC export at 100 m) from the model is well within the very broad range of the literature estimates (0.06–1.0; Table 4).

In the PRESENT_CO2 simulation, coccolithophore biomass is lower in all regions compared to the PRESENT simulation, especially in the ST-SS (north), the subtropical permanently stratified (north), and the equatorial region (Figure 3A). This lower amount of biomass results in a smaller over-estimation of coccolithophore biomass compared to the MAREDAT dataset and, thus, in a better agreement with observations. Differences in diatom biomass between the PRESENT and the PRESENT_CO2 simulation are very small (Figure 3B). Global estimates of biomass, production, and export fluxes decrease slightly compared to the PRESENT simulation, but are still in the range of literature estimates (Table 4). Varying the parameters of the CO₂ dependencies (Table 2) by ±10% can change phytoplankton biomass locally by up to 10% (Figure S3). However, our analysis shows that the parameter set used in the PRESENT_CO2 simulation (Table 2), which is based on function fits to laboratory data, modifies phytoplankton biomass only little (less than 10%) compared to the PRESENT simulation. Hence, we are confident that our approach of the parameter estimation results in a good representation of present-day phytoplankton. Overall, our model represents phytoplankton biomass, distribution, and export fluxes reasonably well in both the PRESENT and the PRESENT_CO2 simulation.

In this section, we first describe spatial patterns of changes in carbonate chemistry (Figure 4) and phytoplankton biomass (Figure 5). Direct and indirect CO₂ effects are then analyzed in four selected regions (Figure 6), namely the three coccolithophore key regions SP-SS + ST-SS (north), equatorial, and ST-SS (south) because of the high sensitivity of coccolithophores to CO₂ changes (Figure 1), and the Southern Ocean region, SP-SS + Ice (south) where the spatial comparison between the simulations reveals a considerable decrease of small phytoplankton and a concomitant increase of diatoms (Figure 5B and C). Temperature effects remain the same in all simulations and, thus, are not discussed. Global estimates of biomass and NPP are integrated over the upper 150 m of the water column. We refer to surface maps and, for the separate regions, to absolute and relative changes at the surface, noting that relative changes are mostly in the same order of magnitude for depth-integrated properties (see Figure S5 for depth distribution of phytoplankton biomass). If the relative changes of surface and depth-integrated biomass differ significantly, we mention them separately.

Phytoplankton biomass and NPP at preindustrial pCO₂ (280 μatm) are compared to values at present-day pCO₂ (420 μatm) and related to changes in carbonate chemistry. On the global scale, largest changes are observed for small phytoplankton with 2.1% lower total biomass and NPP (Table 4), driven primarily by changes in the equatorial and temperate zones (Figure 5C). This coincides with the regions of strongest reduction in surface $HCO3−$ concentrations of up to 8% (about 160 mmol m⁻³; Figure 4B). An exception is the South Pacific subtropical gyre, where considerably lower surface $HCO3−$ concentrations barely affect small-phytoplankton biomass. Total biomass and NPP of diatoms are 2.0% and 2.1% higher, respectively, with lower atmospheric CO₂ concentrations (Table 4), resulting mainly from changes in the same regions as for small phytoplankton (Figure 5B). Globally, biomass and NPP of coccolithophores differ only slightly between preindustrial and present-day pCO₂ (Table 4; Figure 5A). Surface CO_2(aq) is by about 5–10% lower in the preindustrial compared to the present-day simulation, mostly in the Arctic (up to 10 mmol m⁻³) and least in the equatorial regions, and pH is up to 0.2 units higher in the polar regions (Figure 4A and C). Alterations in total phytoplankton biomass largely follow the slightly different pattern of change in surface $HCO3−$ as described above, while CO_2(aq) and pH seemingly are not the main factors for differences in phytoplankton biomass in the two simulations. The (PIC:POC)_cocco ratio is about the same in both simulations (Table 4) because spatial changes in calcite concentration follow changes in coccolithophore biomass (Figure 5D). POC and calcite export are slightly lower by 0.1 Pg C yr⁻¹, resulting in an unaltered rain ratio, and silicate export is 0.3 Tmol yr⁻¹ higher (Table 4).

In the four selected regions, direct growth responses to lower atmospheric CO₂ concentrations are not sufficient to explain the rearrangement of phytoplankton groups. Growth rates of small phytoplankton are as much as 4% lower in the SP-SS + ST-SS (north), the equatorial, and the ST-SS (south) regions at preindustrial atmospheric CO₂ (Figure 6A, C, and E). This reduction of the growth rates is caused by a stronger growth limitation by CO₂ and light, which is partly counteracted by weakened nutrient limitation (Figure 4D and E). Hence, small-phytoplankton biomass in these regions is 3–4% lower (0.1–0.3 Tg C). In the same regions, diatom biomass is 2–5% (0.05–0.1 Tg C) higher. Changes in their growth rates, however, are not caused by lower atmospheric CO₂ directly (Figure 6A, C, and E) but by the modulation of other factors. Coccolithophores are affected most in the SP-SS + ST-SS (north) region, where biomass is reduced by 14% (5% integrated over the entire watercolumn), which does not correspond to the comparatively small change in their growth rate (Figure 6A) and must have been caused by other factors. In comparison to the SP-SS + ST-SS (north) region, coccolithophore growth rates and biomass vary only little in the equatorial and the ST-SS (south) region (Figure 6C and E). As calcification is less inhibited by protons at preindustrial atmospheric CO₂, it can partly counteract the reduction in calcite concentrations caused by lower coccolithophore biomass in the SP-SS + ST-SS (north) region (Figure 6A). In the Southern Ocean, this disproportional effect of higher pH on calcification and on growth also explains why the difference in calcite concentration between present and preindustrial atmospheric CO₂ is stronger than in coccolithophore biomass (Figure 6G). Absolute differences in coccolithophore biomass and calcite, however, are negligible in this region due to the low initial concentrations, and changes in diatom and small-phytoplankton biomass are small (Figure 6G). Generally, differences in NPP follow differences in biomass, with slight discrepancies of up to 5% caused by changes in loss rates (Figure 6A and E). Altogether, small phytoplankton is most negatively and diatoms are most positively affected by preindustrial compared to present-day atmospheric CO₂ concentrations, while the global effects on coccolithophores and calcite are smaller and changes are confined to the regional scale.

Between future atmospheric CO₂ levels (750 μatm) and present-day CO₂ levels (420 μatm), globally averaged phytoplankton biomass and NPP change only marginally (less than 0.4%; Table 4). Except for slightly higher silicate exports at future atmospheric CO₂ concentrations (from 84 to 85.7 Tmol yr⁻¹; Table 4), these marginal changes imply only a small CO₂ effect on globally averaged quantities.

In contrast to the global averages, CO₂ has pronounced regional and group-specific effects (Figures 5 and 6B, D, F, and H). On a regional scale, future coccolithophore growth rates are lower compared to present-day rates because of altered sensitivities to CO₂ and light, with up to 5% stronger limitations in both cases. Weakened nutrient limitation by up to 9%, primarily caused by DIN (Figure 4D) but also by DFe (Figure 4E), as well as the modulation of other factors (see Section Unraveling the cascading effects of changing atmospheric CO₂ levels) regionally compensates the growth-decreasing CO₂ and light effects (Figure 6B, D, and F). As a result, coccolithophore biomass is 12% and 14% higher (2% and 5% in the entire water column), respectively, in the equatorial and the SP-SS + ST-SS (north) region, and remains rather unchanged in the ST-SS (south) region. In all regions, the negative impact of CO₂ on calcification increases by 8–10% (Figure 6B, D, and F). The biomass-related higher calcification and calcite concentrations are thereby dampened in the SP-SS + ST-SS (north) and equatorial regions, and both calcification and calcite concentrations are 8% lower in the ST-SS (south) region where coccolithophore biomass barely changes (Figure 6B, D, and F). Diatoms and small-phytoplankton biomass are barely affected by higher atmospheric CO₂ concentrations in the SP-SS + ST-SS (north), equatorial, and ST-SS (south) regions (Figure 6B, D, and F).

In the Southern Ocean (SP-SS + Ice (south) region), changes in CO_2(aq) and pH are strongest with 10–20% higher surface CO_2(aq) (up to 15–20 mmol m⁻³) and lower pH (up to 0.3 units; Figure 4A and C). Here, more growth-decreasing effects of CO₂ and light, partly counteracted by weakened nutrient limitation (Figure 4D and E), lower small-phytoplankton growth rates by about 3% and biomass by about 5% (Figure 6H). CO₂, light, and nutrient effects on diatoms and the resulting growth rate changes are less pronounced than for small phytoplankton, but have the same sign. Despite that similarity, diatom biomass increases by 8%, implying modulations of other factors (see Section Unraveling the cascading effects of changing atmospheric CO₂ levels). Relative changes in coccolithophore biomass and calcite in the Southern Ocean are negligible due to the generally low biomass and calcite concentrations.

As in the PREIND_CO2 simulations, changes in NPP between the PRESENT_CO2 and the FUTURE_CO2 simulation follow changes in biomass, but vary slightly in magnitude with deviations up to 5% compared to biomass differences (Figure 6B). In summary, the most pronounced changes compared to present-day CO₂ can be seen in the regionally higher coccolithophore biomass and its impact on calcite, as well as a shift to fewer diatoms and more small phytoplankton in the Southern Ocean at future atmospheric CO₂.

We identify three cascading effects triggered by changes in atmospheric CO₂, which can modulate or even counteract the direct CO₂ effect on phytoplankton growth and biomass. CO₂ affects the nutrient and light limitation terms which feed back on the growth rate. The seasonal variation of the nutrient and light feedbacks determine the translation of direct CO₂ effects into biomass. Finally, we found an additional role of cascading effects via top-down factors. In the following, we explain each of the cascading effects for selected phytoplankton groups and regions (Figure 7A).

Coccolithophore growth rates in the FUTURE_CO2 simulation are strongly modified by alleviated nutrient limitation, stronger light limitation (Figure 6B, D, and F), and lower chlorophyll-to-carbon ratios (Figure S6) compared to the PRESENT_CO2 simulation. Nutrient feedbacks are caused by a shift in the community structure and total biomass, leading to changes in available nutrients (Figure 7B). Light feedbacks are more complex (Figure 7C). Firstly, the nutrient limitation term $fN$ modifies the light limitation term $fL$ directly, with a stronger light limitation under more replete nutrient conditions (Equation 5). Secondly, $fN$ modifies the chlorophyll-to-carbon ratio $qChl:C$ in $fL$ (Figure S6), because chlorophyll synthesis is dependent on biomass production of a phytoplankton group, but also varies with nitrogen assimilation and PAR (Hauck et al., 2013) which causes non-proportional changes in cellular chlorophyll and carbon concentrations. Phytoplankton can make less use of a prevailing light condition with a lower $qChl:C$⁠. Finally, PAR has a direct impact on $fL$⁠, and is itself modified by the prevailing chlorophyll concentration in the watercolumn. Thus, lower nutrient limitation comes at the expense of stronger light limitation. Hence, the negative direct impact of high CO₂ levels on coccolithophores is modulated by these two feedbacks.

The nutrient and light feedbacks imposed by CO₂ dependencies vary seasonally, and manifest themselves most strongly in changes in biomass and NPP during the growing season (Figure 7B and C). For instance, while the negative effect of CO₂ and light on coccolithophore growth in the FUTURE_CO2 simulation can be counteracted partly by decreasing nutrient limitation in the annual mean, this effect cannot fully explain the considerably higher coccolithophore biomass in the SP-SS + ST-SS (north) region (Figure 6B). Monthly differences between the FUTURE_CO2 and the PRESENT_CO2 simulations show that coccolithophore growth rates are up to 5% higher between May and September at high atmospheric CO₂ (Figure 8A), caused by lower light and nutrient limitation. As this period roughly coincides with the growing season with high initial biomass, it causes an annual mean increase in biomass and NPP. Lower annual mean growth rates at future compared to present-day CO₂ levels are caused by lower growth rates in winter (10–12% between October and April) due to a stronger light limitation and a relatively smaller alleviation in nutrient limitation compared to the remaining year, a time where initial biomass levels are low and changes in biomass and NPP are negligible. Annual means of growth rates are therefore often insufficient to explain changes in NPP and biomass, and the seasonality of cascading effects needs to be considered.

Finally, we identified modifications in the grazing pressure as another indirect CO₂ effect on phytoplankton biomass. Due to the variable grazing preference in our model (Fasham et al., 1990; Vallina et al., 2014), grazing fluxes change non-linearly with variations in the share of each phytoplankton group to total biomass as well as with changes in the total biomass of zooplankton (Figure 7D). Despite the higher diatom biomass in the FUTURE_CO2 compared to the PRESENT_CO2 simulation (about 7%) due to increasing growth rates during summer in the SP-SS + Ice (south) region and the corresponding higher prey availability, grazing on diatoms decreases by about 6% in the respective months (Figure 8B). This decrease is caused by a lower zooplankton biomass, related to the decreasing small-phytoplankton biomass (Figure 6H), a food source that cannot be replaced entirely by diatoms (note grazing preferences in Section Implementing coccolithophores into REcoM). Lower grazing on diatoms in austral winter (June–August; Figure 8B) is a result of less diatom biomass. Grazing feedbacks also explain higher diatom biomass concentrations in the PREIND_CO2 simulation in all regions except from the Southern Ocean (Figure 6A, C, and E). Another example from the preindustrial SP-SS + ST-SS (north) region is the increasing grazing pressure on coccolithophores and their consequently lower biomass, which results from a combination of relatively high coccolithophore biomass in that region and biomass-decreasing effects of lower CO₂ on small phytoplankton which are then less available as prey. Other loss terms of phytoplankton biomass (respiration, excretion, aggregation) depend directly on the total phytoplankton and detritus biomass and/or biomass of the respective phytoplankton group and are not or negligibly affected by shifts in the share between phytoplankton groups. Hence, grazing on a phytoplankton group can be alleviated if the biomass of a zooplankton group decreases because its preferred prey is not present in sufficient abundance. In summary, cascading effects of CO₂ on light and nutrient limitations, on their seasonality, and on the grazing pressure are pivotal to the biomass and NPP responses of the phytoplankton groups.

In our second set of simulations (Table 3) we aimed to disentangle the effects of warming and CO₂ on historical changes in phytoplankton biomass. Globally, depth-integrated coccolithophore biomass levels in the VARCLI_VARCO2 simulation are highly variable in the simulated period (between 24 and 27 Tg C) with lower biomass in the 1990s–2000s compared to the 1960s–1970s (24–25 Tg C and 26–27 Tg C, respectively; Figure 9A). Over this time period, temperatures and CO_2(aq) concentrations increase by about 0.5°C and 0.003 mmol m⁻³, respectively (Figure 9B). Fluctuations in coccolithophore biomass levels are small at constant climatological forcing (CONSTCLI_VARCO2 simulation), implying that interannual variability of coccolithophore biomass is driven mainly by climate variability and, thus, combined alterations in temperature and CO_2(aq) concentrations. The lower CO_2(aq) concentration in the VARCLI_CONSTCO2 simulation compared to the VARCLI_VARCO2 simulation (–0.001 mmol m⁻³ at the beginning and increasing over the course of the time series) increases the coccolithophore biomass by almost 1 Tg C in the entire time period. However, increasing CO_2(aq) concentrations over time at constant temperature have almost no effect on the integrated coccolithophore biomass (CONSTCLI_VARCO2 simulation). Hence, according to our simulations, CO_2(aq) is not the main driver for global biomass changes since the 1960s, and variability in global coccolithophore biomass is driven mainly by climate variability.

In the North Atlantic, mean coccolithophore biomass levels in the 2000s are 4% higher than biomass levels in the 1960s (2.4 and 2.3 Tg C, respectively; Figure 9C). Our model therefore simulates a much smaller increase in coccolithophore biomass in the North Atlantic between 1960 and 2010 than the up to 10-fold increase in coccolithophore occurrence in observations (Beaugrand et al., 2013; Rivero-Calle et al., 2015; Krumhardt et al., 2016). After 2010, biomass levels decrease quickly by 0.2 Tg C, reaching biomass levels of 2.2 Tg C in 2019. Temperatures increase rapidly from about 11.5°C before the 1990s to about 12.25°C by the end of the time series with strong interannual variability, whereas CO_2(aq) concentrations increase steadily from about 12 × 10⁻³ mmol m⁻³ in the 1960s to about 15.5 × 10⁻³ mmol m⁻³ in the 2010s (Figure 9D). Thus, just as on a global scale, interannual biomass variability is caused by climate variability. In contrast to what is simulated on a global scale, however, biomass levels are increasingly higher over the course of the time series at increasing CO_2(aq) concentrations (VARCLI_VARCO2 simulation) than at constant CO_2(aq) concentrations (VARCLI_CONSTCO2 simulation, up to 0.1 Tg C in 2019; Figure 9C). This result is most likely due to a growth enhancing effect of higher CO_2(aq) concentrations, which are about 6 × 10⁻³ mmol m⁻³ higher in the VARCLI_VARCO2 compared to the VARCLI_CONSTCO2 simulation at the beginning of the time series, with increasing difference towards the end of the time series (Figure 9D). However, biomass levels barely change with increasing CO_2(aq) concentrations and constant temperature (CONSTCLI_VARCO2 simulation), pointing towards a small effect of CO_2(aq) on coccolithophore biomass in the period of the time series. We therefore attribute biomass changes predominantly to climate variability.

Our model simulations show that adding a CO₂ dependence to phytoplankton growth leads to a cascade of effects on (i) the growth rate of the individual phytoplankton groups by feedbacks on light and nutrient limitations and the seasonality therein (bottom-up effects), and (ii) the interaction between phytoplankton and zooplankton (top-down effects). In the preindustrial global ocean, small phytoplankton are affected most by direct and indirect CO₂ effects. The lower biomass of small phytoplankton in low to mid-latitudes causes lower zooplankton biomass concentrations and a shift of the grazing pressure to coccolithophores, which results regionally in up to 14% lower coccolithophore biomass. Diatoms, in contrast, benefit from lower zooplankton biomass rather than suffering from higher grazing pressure due to chosen grazing preferences, but effects are rather localized. In the future global ocean, coccolithophores are affected most by direct and indirect CO₂ effects, resulting in higher coccolithophore biomass. Changes in coccolithophore biomass affect the total phytoplankton biomass only marginally and therefore do not feed back on zooplankton biomass and grazing rates. In high latitudes, bottom-up effects cause lower small-phytoplankton biomass with consequently lower zooplankton biomass, which decreases the grazing pressure on diatoms. In low to mid-latitudes, bottom-up effects trigger slightly higher growth rates and biomass concentrations of small phytoplankton.

These results highlight two important aspects that shape top-down effects. Firstly, the standing stock of a group determines whether changes therein matter for the overall zooplankton prey availability (e.g., small phytoplankton versus coccolithophores). Secondly, the grazing preference determines if and how zooplankton can switch to a different prey (e.g., small phytoplankton and coccolithophores versus diatoms). The small zooplankton group in our model grazes preferentially on small phytoplankton and coccolithophores, and the polar macrozooplankton grazes mainly on diatoms. While the spatial distribution of polar macrozooplankton is confined to south of 50°S and north of 50°N, this group out-competes here the small zooplankton group due to its higher grazing efficiency. In the remaining ocean, grazing of the small-phytoplankton group dominates. A variety of grazing formulations are used in biogeochemical models which differ in how they describe prey switching and total grazing pressure, suggesting that using a different grazing parametrization can change the simulated grazing pressure on individual phytoplankton groups (Vallina et al., 2014). Nonetheless, given the good agreement of the simulated biomass fields in our simulations with observational data (Figures 3A and B and S4), we are confident that the qualitative dynamics of the feedbacks triggered by the addition of CO₂ dependencies also hold across grazing parameterizations. The complexity of modelled trophic interactions with partially unexpected feedbacks was also highlighted by Dutkiewicz et al. (2021), who showed that a small reduction in the abundance of one species can increase the total phytoplankton standing stock due to shifts to slower growing grazers and the freeing up of limiting nutrients, which is comparable to the cascading effects in our model. In the following, we discuss regional changes in phytoplankton biomass concentrations.

In the preindustrial global ocean, small-phytoplankton biomass is lower compared to present-day CO₂ levels. Increasing levels of small-phytoplankton biomass over the last decades have also been shown in field and laboratory studies, with the projection of further increasing biomass with ongoing climate change (Daufresne et al., 2009; Peter and Sommer, 2012; Pinkerton et al., 2021). This beneficial effect has been associated with joint warming and increasing CO₂ (Hare et al., 2007; Feng et al., 2009). However, our model results suggest that even the carbonation effect alone, mainly caused by the large growth increase from preindustrial to present-day CO₂ levels, can lead to increasing biomass of small phytoplankton.

In the future Southern Ocean our model reveals a shift to less small phytoplankton and more diatom biomass compared to present-day CO₂ levels (Figures 5B and C and 6H). This advantage of diatoms over small phytoplankton is driven mainly by decreasing grazing pressure on diatoms (Figure 8B). In nature, the advantage of diatoms at future oceanic conditions may be further enhanced by alterations in the diatom species composition. Laboratory and mesocosm studies reveal that under ocean acidification, diatoms likely shift towards larger species (Tortell et al., 2008; Wu et al., 2014; Bach et al., 2019), changing community size structures and composition and having implications on the food web and biogeochemical cycling (Finkel et al., 2009; Alvarez-Fernandez et al., 2018). As currently we do not include cell size structure within a phytoplankton group in our model, accounting for shifts in this important trait (Dutkiewicz et al., 2015) could be a step forward to assess implications on export fluxes, with POC export likely increasing due to the higher biomass and ballasting of large diatom cells. While in parts of the Southern Ocean, a deepening of the mixed layer depth has been observed and modelled for recent years (Hauck et al., 2015; Panassa et al., 2018; Sallée et al., 2021) and future projections (Hauck et al., 2015), other regions may experience a temperature-driven stronger stratification (Constable et al., 2014). Small phytoplankton, because of their lower nutrient requirement, typically follow Southern Ocean diatom blooms that have drawn down nutrient concentrations (Irion et al., 2021). Similarly, a shoaling of the mixed layer would be more advantageous to small phytoplankton than for diatoms (Dutkiewicz et al., 2015; Petrou et al., 2016). Our results indicate, however, that increasing CO₂ levels may dampen the nutrient-driven advantage of small phytoplankton over diatoms in the future Southern Ocean (Figure 6H).

While laboratory CO₂ manipulation experiments tend to show either unaffected or decreasing coccolithophore growth rates between present-day and future CO₂ levels (e.g., Bach et al., 2013; Seifert et al., 2020; but see Meyer and Riebesell, 2015, for the diversity of responses), coccolithophore biomass in our model is increasing due to the unexpected cascading effects that offset the negative effects of CO₂ on growth (Figures 5A and 6B, D, and F). Interactive effects between warming and high CO₂ conditions indicate a shift of the turning point between carbonation and acidification towards higher CO₂ levels at warmer temperatures (e.g., Sett et al., 2014; Gao et al., 2018), which may reduce the negative CO₂ effect on coccolithophores in the future scenario and ultimately lead to higher coccolithophore biomass, just as in our CO₂-only simulation. However, while the CO₂ limitation in our model showed only modest seasonal variations, mesocosm experiments have shown that a more severe CO₂-induced growth reduction at the beginning of the growing season can impede bloom formation and thereby reduce the overall ability to maintain sufficiently large seed populations for the ensuing growth seasons (Riebesell et al., 2017). Besides, a climate change-related increase in upper ocean stratification and the concomitant decrease in nutrient supply could increase the competitive fitness of small phytoplankton, as explained in Section Cascading effects cause fewer small phytoplankton in the preindustrial ocean and more diatoms in the future Southern Ocean, which could outweigh the negative effect of ocean acidification. Less nutrients are then available for diatoms and coccolithophores, which decreases their competitive fitness, an effect that cannot be seen in simulations where only atmospheric CO₂ is changed. Consequently, the negative CO₂ effect on coccolithophores could become more prevalent in a joint warming and high-CO₂ scenario, and decrease coccolithophore biomass in the future.

Our model suggests that the CO₂ impact on calcification, which is positive in the preindustrial and negative in the future ocean compared to present-day conditions, is modified and partly compensated by changes in coccolithophore biomass. Higher calcification rates under preindustrial CO₂ levels compared to the present-day in our study (Figure 6C, E, and G) seem to disagree with findings of Rigual-Hernández et al. (2020) who found that the thickness of the coccoliths barely differs between preindustrial samples, originating from sediment cores, and modern plankton samples. However, coccolith thickness and the sensitivity to CO₂ is highly species-specific (Rigual-Hernández et al., 2020), and shifts in coccolithophore species composition can impact patterns in calcite concentration and export (Beaufort et al., 2011). The number of coccoliths per cell is also variable, at least for E. huxleyi, and can additionally affect calcite concentrations under changing environmental conditions (Paasche, 2001). Whether or not there was a global trend towards lower calcification in all coccolithophore species from preindustrial to present-day remains unresolved.

In our simulations, calcification is the most sensitive process towards future CO₂ levels (Figure 1). Calcite concentrations increase in some regions (SP-SS + ST-SS (north) and equatorial) only because of increasing coccolithophore biomass (Figures 5A and 6B and D). Based on the ballasting hypothesis (e.g., Klaas and Archer, 2002; Cram et al., 2018), this increase would be indicative of an increase of export production, which is currently not represented in our model. Furthermore, while on a global scale, enhanced calcification and the accompanying CO₂ evolution affects the anthropogenic increase in atmospheric CO₂ only little (Gangstø et al., 2011), it can shape the regional air-sea CO₂ flux (Shutler et al., 2013).

Increasing phytoplankton biomass and decreasing calcification at future CO₂ levels in our model compare well with the outcome of other modelling studies that include a CO₂ dependence of growth rates and/or calcification (Gangstø et al., 2011; Dutkiewicz et al., 2015; Krumhardt et al., 2019). Similar to our model, considerable alterations in global phytoplankton communities attributable to different sensitivities towards increasing CO₂ levels were simulated in the ecosystem model DARWIN of Dutkiewicz et al. (2015), whereas warming alone caused a poleward shift in phytoplankton dispersal with no significant impact on the community structure. While we used a mechanistic approach, Dutkiewicz et al. (2015) used a conceptual approach with a linear and positive dependence of growth rates on CO₂. They showed that decreasing nutrient availability due to stronger stratification lowers globally integrated primary production by 5% in 2100, which can be counteracted by CO₂-enhanced growth rates, leading to almost unchanged global primary production levels in 2100 compared to present-day levels. Similarly, Krumhardt et al. (2019) attributed global coccolithophore biomass increase to a direct CO₂ effect in their model. Contrary to these two studies which only account for the carbonation effect, we also modelled the acidification effect on phytoplankton growth whereby the direct CO₂ effect causes a decrease in future growth rates. However, cascading effects through nutrient and light feedbacks (Figure 7) can lead to an overall increase in biomass. While we did not simulate changes in circulation and stratification, we hypothesize that reduced nutrient availability in a more stratified ocean could reinforce rather than counteract the CO₂-driven decrease in growth. Thus, we rather expect a future decrease in phytoplankton biomass, which cannot be balanced by CO₂, different to the idealized simulations of Dutkiewicz et al. (2015).

In our model, increasing future coccolithophore biomass due to cascading effects over-compensates for the pH-driven decrease in calcification (Figure 6B and D), resulting in a net increase of global calcite concentrations that corresponds to the findings of Krumhardt et al. (2019). Both our study and that of Krumhardt et al. (2019) are in line with Gangstø et al. (2011), who projected a negative direct CO₂ effect on calcite production in a high emission scenario of the biogeochemical Bern3D/PISCES model. Going beyond the study of Gangstø et al. (2011) by modelling CO₂ effects on coccolithophore biomass, we have shown that biomass changes can compensate for decreasing calcite production, which is also in line with Krumhardt et al. (2019). In summary, our findings highlight the need to consider both carbonation and acidification effects on phytoplankton growth and coccolithophore calcification.

According to our model, we attribute observed changes in North Atlantic coccolithophore biomass in recent decades primarily to warming and not to CO₂. Thus, our results largely comply with the warming hypothesis of Beaugrand et al. (2013), as changes in the modelled carbonate system during the 60 years of the time series are too small to have a pronounced effect on coccolithophore biomass. We thereby contrast the CO₂ hypotheses of Rivero-Calle et al. (2015) and Krumhardt et al. (2016), who concluded that warming alone did not control changes in coccolithophore biomass. Even with more significant shifts in the carbonate system, the negative impact of CO₂ on coccolithophore growth obtained by our model are over-compensated by more beneficial cascading effects on nutrient limitation (Section Unraveling the cascading effects of changing atmospheric CO₂ levels). Admittedly, none of the temperature functions in our model includes decreasing growth rates caused by temperature stress, which can considerably affect community structures (e.g., D’Amario et al., 2020), but we suppose that thermally stressed species or ecotypes would be replaced by more heat-tolerant species or ecotypes of the same phytoplankton group, resulting in the same changes in coccolithophore biomass as simulated here. We expect thermal stress to become more relevant in polar regions and/or for a more intense temperature increase that outpaces the speed of species migration, but not for the given temperature increase within 60 years in the North Atlantic. We additionally need to bear in mind that our parameterization aims to represent the diverse group of coccolithophores, while the North Atlantic coccolithophore community is dominated by the fast growing, bloom-forming E. huxleyi (Haidar and Thierstein, 2001; Balch et al., 2019), which may respond more strongly to changing CO₂ conditions than coccolithophores in our model. Besides, the CO₂ effect on coccolithophore growth would be different if extreme CO₂ levels of laboratory experiments had been excluded from the parameterization. Based on our model results, however, we conclude that the direct CO₂ impact on coccolithophore biomass in the recent decades, globally and in the North Atlantic, is small compared to the effect of warming, supporting the hypothesis of Beaugrand et al. (2013).

Our description of CO₂ dependencies goes well beyond previous studies by developing mechanistic response functions of phytoplankton functional type growth rates and coccolithophore calcification. We are confident that it captures the first-order effects of CO₂ on phytoplankton biomass and NPP. Nevertheless, we acknowledge that additional processes (i.e., multiple driver interactions, mixotrophy, CO₂ effects on silicification and zooplankton, adaptation and evolution, as well as mesozooplankton grazing) have the potential to further modify the response of the marine ecosystem to high-CO₂ levels in simulations and in the real world.

The response of phytoplankton growth and calcification to changing CO₂ levels can be modified significantly by synergistic and antagonistic interactions with other environmental drivers such as temperature and light (e.g., Harvey et al., 2013; Brandenburg et al., 2019; Seifert et al., 2020), and can result in enhanced or dampened responses than expected from CO₂ alone. However, most biogeochemical models, including ours, do not account for the interaction between environmental drivers. For instance, acidification affects the response to nutrient availability of coccolithophores (Zondervan, 2007; Zhang et al., 2019) and other phytoplankton groups (Li et al., 2012; Spungin et al., 2014; Li et al., 2017). The beneficial effect of nutrients on coccolithophore biomass in the FUTURE_CO2 simulation may thus be weakened by acidification, and therefore dampen the biomass increase compared to present-day conditions. The implementation of multiple driver interactions would therefore be a crucial amendment to phytoplankton modelling.

Mixotrophy, which is the capacity of an organism to live auto- and heterotrophically, can reduce the vulnerability of phytoplankton to environmental changes, for example, under substrate limitation (low carbonation) of photosynthesis. Mixotrophy is very common in species that we summarize as small phytoplankton (Stoecker et al., 2017), but some coccolithophore species can live mixotrophically as well (Godrijan et al., 2020). While our phytoplankton growth rate parameterization is based entirely on photosynthesis, effects of changing CO₂ levels could be diminished if we additionally considered the possibility to switch to mixotrophy when CO₂ or other conditions become unfavorable. Besides, not only calcification, but also silicification by diatoms appears to be CO₂ sensitive (Petrou et al., 2019). Although this sensitivity could alter the response of diatoms towards increasing CO₂ conditions and, for instance, possibly hamper the success of diatoms in the future Southern Ocean, the knowledge about this relationship is still too limited to be incorporated into modelling approaches. Similarly, CO₂ effects on grazers and higher trophic levels could further modify the ecosystem response (Cripps et al., 2014).

Adaptation and evolution can shape substantially the response of phytoplankton towards high or low pCO₂ levels (Brennan et al., 2017), processes that are not usually considered in ocean biogeochemical models. For instance, coccolithophores can adapt within months or years to warming and high CO₂ conditions (Schlüter et al., 2014). However, with the current pace of environmental changes caused by climate change and the limits of adaptation, whether adaptational and evolutionary processes are fast enough to offset negative effects of environmental driver changes is questionable. The assessment of abilities and speed of phytoplankton (and other organism) adaptations is challenging (Kelly and Griffiths, 2021), and the data basis remains highly varied among studies, which currently prevents a robust implementation into biogeochemical models.

While our model represents a polar macrozooplankton as well as a generic small zooplankton group, it does not include mesozooplankton, which is globally distributed (Moriarty and O’Brien, 2013) and an important component of trophic interactions and the biological carbon pump (Buitenhuis et al., 2006). Hence, considering mesozooplankton could alter the grazing feedbacks presented in our study, for instance in the Southern Ocean where the grazing pressure on diatoms is indirectly lowered due to the decrease of small-phytoplankton biomass.

Model projections are the primary tool to understand future shifts in marine productivity and its impact on higher trophic levels and food webs. Even if the resilience of phytoplankton communities to climate change-related impacts might be relatively high because environmental niches can be re-occupied by species that are better adapted to the new environmental conditions (Dutkiewicz et al., 2021), most models project a future decrease in global NPP, with a large spread in the magnitude across different models (Kwiatkowski et al., 2020). Carbonation and acidification may have a relatively smaller effect on phytoplankton growth than other environmental drivers such as warming or changes in nutrient availability, but these processes will become increasingly important with ongoing ocean acidification. Our study highlights that CO₂, especially due to its cascading effects on light and nutrient limitation and grazing, can change regional phytoplankton community compositions and primary production considerably in a simulation with future atmospheric CO₂ concentrations. Modifications in biomass concentrations can locally account for 10% and more. In addition, other studies reveal that CO₂ can tip the scale when multiple driver interactions are involved (e.g., Harvey et al., 2013; Brandenburg et al., 2019; Seifert et al., 2020). Accounting for carbonation and acidification effects on phytoplankton growth and calcification in model projections adds an important physiological constraint on primary production and biogeochemical cycling in the future ocean.

Model data necessary to produce the findings of this study are archived in Zenodo (https://zenodo.org/; Seifert and Hauck, 2022). Full model output can be requested from the corresponding author.

The supplemental files for this article can be found as follows:

Figures S1–S6.pdf

Text S1.pdf

Tables S1–S4.pdf

We thank Meike Vogt from ETH Zurich, Switzerland, and Colleen O’Brien for providing gridded data products of coccolithophore biomass concentrations for our model evaluation. Furthermore, we thank Frédéric Maps and an anonymous reviewer for their constructive comments and suggestions.

This research was supported under the Initiative and Networking Fund of the Helmholtz Association (Helmholtz Young Investigator Group Marine Carbon and Ecosystem Feedbacks in the Earth System [MarESys], grant number VH-NG-1301) as well as the European Union’s Horizon 2020 research and innovation program under grant agreement number 869357 (project OceanNETs: Ocean-based Negative Emission Technologies—analyzing the feasibility, risks, and co-benefits of ocean-based negative emission technologies for stabilizing the climate). MS acknowledges funding from the Helmholtz Graduate School for Polar and Marine Research POLMAR. CN has received funding from the European Union’s Horizon 2020 research and innovation program under grant agreement number 820989 (project COMFORT). The work reflects only the authors’ views; the European Commission and their executive agency are not responsible for any use that may be made of the information the work contains.

The authors have no competing interests, as defined by Elementa, that might be perceived to influence the research presented in this article.

Contributed to conception and design: MS, CN, BR, JH.

Contributed to acquisition of data: MS, CN, BR, JH.

Contributed to analysis and interpretation of data: MS, CN, BR, JH.

Drafted and/or revised the article: MS, CN, BR, JH.

Approved the submitted version for publication: MS, CN, BR, JH.