Complexity in biogeochemical models: consequences for the biological carbon pump

Ocean biogeochemical models underpin projections of future marine ecosystem change, including anticipated shifts in the biological carbon pump (BCP) and broader biogeochemical cycles. However, their outputs remain highly sensitive to model complexity and parameterisation choices. Here, we evaluate five configurations of the Pelagic Interaction Scheme for Carbon and Ecosystem Studies (PISCES) to quantify intramodel variability in net primary productivity (NPP), carbon export (C_exp), and export efficiency (e-ratio) over the 21st century under the high emissions RCP8.5 scenario. The tested PISCES configurations differed from the standard model through distinct modifications to phytoplankton growth processes, but are forced by identical physical variables, representing an ensemble opportunity. All configurations resolve NPP and C_exp within the range of remote-sensing variability. The more complex Quota-based configurations produce 15–21 (10–18) Pg C yr⁻¹ more NPP than the simpler Monod-quota models in the reference (future) period, but this increase, driven by elevated small phytoplankton biomass, does not enhance C_exp, yielding lower e-ratios (0.14–0.17) than in the Monod-quota configurations (∼0.25). The introduction of a picophytoplankton functional type (PFT) emerges as one of the most influential parameterisation choices. It drives opposing future NPP responses between 30–60° N/S, an increase in the Monod-quota configurations versus a decline in the Quota-based ones, as well as contrasting latitudinal trends in C_exp within the same region. Other parameterisations, such as a low-iron scheme, an added diazotroph PFT, and explicit manganese cycling, exert more modest, regionally confined effects under high emissions scenarios, influencing NPP and C_exp primarily at biome scales rather than driving large-scale divergence in model behaviour.

1 Introduction

The biological carbon pump (BCP) plays a critical role in regulating Earth's climate by facilitating the vertical transport of carbon from the surface to the deep ocean, thereby helping to remove atmospheric CO₂ on centennial to millennial timescales (Falkowski et al., 1998, 2000; Buesseler et al., 2007; Visser, 2025). Since pre-industrial times (∼280 ppm), global atmospheric CO₂ concentrations have risen by 51 %, reaching 419 ppm as of 2023 (Friedlingstein et al., 2023). This increase has contributed to a ∼1.5 °C rise in average global surface temperatures, with the previous decade (2014–2024) marking the warmest on record (Ripple et al., 2024; WMO, 2025). The ocean has absorbed 20 %–30 % of anthropogenic CO₂ emissions since the mid-1980s (Sabine et al., 2004; Bindoff et al., 2019), equivalent to 2.1–2.4 Pg C yr⁻¹ (DeVries et al., 2023). Stand-alone ocean model studies suggest that, at equilibrium, BCP-mediated processes result in baseline atmospheric CO₂ concentrations being 150–200 ppm lower than they would be otherwise in their absence (Maier-Reimer et al., 1996; Falkowski et al., 2000). A fully interactive Earth System Model (ESM) refines this estimate, showing a 163 ppm increase in pre-industrial CO₂ when marine ecosystems are removed (Tjiputra et al., 2025). The BCP therefore represents a substantial ecosystem service, valued at USD 0.8–1.1 trillion per year globally (Berzaghi et al., 2025). Understanding how the biological processes that underpin the BCP, and its efficiency, will respond to climate change remains an important and active area of research (Jin et al., 2020).

Through air–sea gas exchange, atmospheric CO₂ dissolves into the surface ocean, increasing the dissolved inorganic carbon (DIC) pool. In the ocean, this DIC is absorbed and sequestered through two mechanisms, the solubility pump and BCP. The solubility pump is driven by the combination of thermodynamic and physical processes that slowly transport DIC from the surface ocean into the interior via mode and deep water formation. It is the dominant mechanism for the uptake and transport of anthropogenic CO₂ from the surface ocean to below the mixed layer (Bindoff et al., 2019). The BCP further facilitates this vertical transport of carbon through the export of particulate organic carbon (POC), synthesised via phytoplankton net primary productivity (NPP) in the sunlit surface ocean. The subsequent sinking and export (C_exp) of POC out of the euphotic zone and into the deep ocean (>1000 m) contributes to the long-term removal of carbon on timescales of hundreds to thousands of years (Falkowski et al., 1998, 2000; Buesseler et al., 2007; Visser, 2025). Phytoplankton photosynthesis lowers the partial pressure of CO₂ in the surface ocean while export production sustains the vertical gradients in DIC necessary to promote the absorption of CO₂ from the atmosphere. Therefore, the BCP emerges as a critical, climate-relevant pathway in the ocean's role as a long-term carbon sink (Falkowski et al., 2000; Sigman and Boyle, 2000).

In the upper ocean, most of the produced POC is rapidly remineralised, with only a small fraction (∼20 %) being exported below 100 m (Falkowski et al., 2000; Laws et al., 2000; Ducklow et al., 2001; Henson et al., 2012; Siegel et al., 2014; Nowicki et al., 2022). The efficiency of the BCP, commonly denoted as the e-ratio ($={\mathrm{C}}_{\text{exp}}/\text{NPP}$), principally determines how much fixed carbon is exported from the surface ocean to depth, indirectly influencing the long-term sequestration of atmospheric CO₂. A consensus emerging from multi-model intercomparison studies is that global warming will lead to higher sea surface temperatures (SSTs), increased stratification, and reduced oxygen solubility by the end of the 21st century (Steinacher et al., 2010; Bopp et al., 2013; Bindoff et al., 2019; Kwiatkowski et al., 2020). In response to these changing physical environmental conditions, marine ecosystems and phytoplankton community composition (Bopp et al., 2005; Basu and Mackey, 2018; Tréguer et al., 2018; Benedetti et al., 2021; Henson et al., 2021; Frémont et al., 2022; Fisher et al., 2025) are expected to be impacted, which will have consequences for ocean productivity and export production.

Capturing the complexity of the BCP remains a major challenge in ESMs (Doney et al. 2024). Across successive generations of Coupled Model Intercomparison Project(s) (CMIP), there is substantial intermodel variability in projected 21st century NPP and C_exp (Bopp et al., 2013; Laufkötter et al., 2015, 2016; Kwiatkowski et al., 2020; Tagliabue et al., 2021; Wilson et al., 2022), reflecting large differences in how models represent key biogeochemical processes and marine ecosystem complexity (Séférian et al., 2020; Henson et al., 2022). Under the high-emissions SSP5-8.5 scenario, CMIP6 models project that export production will change by −41 % to +1.8 % by the end of the 21st century relative to the pre-industrial period (Henson et al., 2022) and relative to present-day (1995–2014) conditions, global NPP is projected to decline by 1.76 %, with a model standard deviation of 9.01 % (Tagliabue et al., 2021). Furthermore, future declines in NPP could likely be underestimated within the CMIP6 ensemble, considering their contemporary trends relative to remote-sensing (Ryan-Keogh et al., 2025).

Within the CMIP ensemble, a wide range of biogeochemical models is represented (e.g. BEC – Moore et al., 2001, 2004; BFM – Vichi et al., 2007; PISCES – Aumont et al., 2015; REcoM2 – Hauck et al., 2013; TOPAZ – Dunne et al., 2012), each differing in its complexity and parameterisation of key biogeochemical processes. Although model intercomparison projects have provided valuable insights into the global BCP, it remains unclear how differences in these biogeochemical parameterisations, and complexity thereof, contribute to the divergence among CMIP models in their projections of the BCP over the 21st century. Most model intercomparison studies evaluate different architectures against common criteria to identify model shortcomings and ensemble-wide patterns (Fennel et al., 2022). However, pinpointing sources of divergence and the role of individual parameterisations is challenging within this context. An alternative is to compare different parameterisations within the same biogeochemical model under identical physical forcing, enabling variability arising solely from biogeochemical differences to be isolated.

This study employs the PISCES ocean biogeochemical model (Aumont et al., 2015) in a suite of distinct configurations under high emissions RCP8.5 scenario to isolate how differences in biogeochemical parameterisations and complexity drive variability in the BCP over the 21st century. The methodology details the various PISCES configurations and their nuances while the results are presented in two parts. In the first, we compare present-day PISCES outputs with remote-sensing data to evaluate whether variations in biogeochemical parameterisations and complexity lead to discernible differences in NPP and export production at global and biome scales. The second examines how these variables evolve between present/reference (1986–2005) and future (2091–2100) projections, highlighting how differences in top-down processes propagate through NPP and C_exp. Ultimately, the study identifies key biogeochemical parameterisations that strongly shape ecosystem dynamics, driving differences in NPP and export production, and contributing to intramodel variability in projections of the BCP.

2 Methodology

2.1 PISCES configurations

Table 1 presents a summary of the five PISCES configurations used in this study, which differ primarily in their parameterisation of phytoplankton growth processes. The standard PISCES (PST) model of Aumont et al. (2015) served as the foundation for all configurations. The PST model simulates the lower trophic levels of marine ecosystems and resolves the biogeochemical cycles of carbon and the main nutrients of P, N, Si and Fe. The model has four living compartments, two phytoplankton functional types (PFTs; diatoms and nanophytoplankton) and two zooplankton size classes (microzooplankton and mesozooplankton). Particulate detritus, produced through phytoplankton–zooplankton–nutrient interactions, is partitioned into small and large particles. Both size classes arise from multiple processes, including zooplankton grazing, phytoplankton and zooplankton mortality, and aggregation. However, their dominant sources differ as small particles are primarily linked to nanophytoplankton–microzooplankton interactions, whereas large particles are mainly associated with those of diatom-mesozooplankton. The PST model employs a hybrid of Monod-quota formulations in representing the cycling of major nutrients. Phytoplankton C : N : P stoichiometry is fixed, with N and P limitations governed by a Monod parameterisation. As a result, growth rates follow a rectangular hyperbolic function of the limiting external nutrient concentration (Monod, 1949; Flynn, 2003). In contrast, Fe limitation is represented with a quota-based approach (Aumont et al., 2015), where growth depends on the internal nutrient quota of a cell, which itself is regulated by the external nutrient concentration (Droop, 1968).

Table 1 List of PISCES configurations.

The PST and low-Fe PISCES (PSF) configurations are those described in Tagliabue et al. (2020). The PSF configuration modified PST by adjusting biological and maximum iron quotas for phytoplankton growth to better reflect low-iron conditions and their implications under climate change. Iron availability plays a key role in regulating phytoplankton growth and global NPP (Tagliabue et al., 2017); however, ESMs often show large uncertainties in representing the ocean iron cycle and its effect on productivity (Tagliabue et al., 2016). Unlike PST and PSF, the PISCES-QUOTA (P5Z) model allows for variable C : N : P stoichiometry (Kwiatkowski et al., 2018), with growth exclusively governed by quota-based formulations (Table 1). P5Z resolves 39 prognostic variables, compared to 24 in PST, and importantly includes three phytoplankton groups: picophytoplankton, nanophytoplankton, and diatoms.

Wrightson et al. (2022) further developed the model by introducing an explicit diazotroph PFT (Crocosphaera and Trichodesmium), along with temperature-dependent nutrient uptake parameterised via elemental use efficiencies (EUEs), resulting in the P6Z configuration. These additions addressed a key limitation in P5Z, which represented nitrogen fixation as an implicit flux of ammonia, preventing mechanistic investigation of diazotroph dynamics. Given that marine nitrogen fixation plays a critical role in shaping the response of NPP to climate change (Wrightson and Tagliabue, 2020; Bopp et al., 2022), explicitly resolving diazotrophs in P6Z allowed responses to environmental drivers, especially rising ocean temperatures, to be investigated. Preliminary work conducted in this study showed little difference in global C_exp or NPP patterns across the various P6Z versions in Wrightson et al. (2022). Therefore, in this study, we used only the configuration that included temperature dependent EUEs for Fe and P for the generic N-fixing PFT (Table 1).

Anugerahanti and Tagliabue (2024) built upon the base P5Z model by incorporating manganese (Mn), following Hawco et al. (2022), which included its role in limiting phytoplankton productivity in the Southern Ocean, where observations show Mn as either a primary or co-limiting micronutrient alongside Fe (Wu et al., 2019; Browning et al., 2021; Hawco et al., 2022). The resulting P5M configuration explicitly simulates coupled Mn and zinc (Zn) cycles (Richon and Tagliabue, 2021), including biological uptake, scavenging, and regeneration for Mn and Zn.

The various PISCES configurations (Table 1) were identically forced with physical outputs from the IPSL-CM5A-LR (https://doi.org/10.1594/WDCC/CMIP5.IPILr8) climate model on the ORCA tripolar grid with 2° horizontal spatial resolution. Two model runs were conducted for each configuration. The first simulation was forced by outputs derived from historical (1851–2005) climate forcing and the second with the high emission Representative Concentration Pathway 8.5 (RCP8.5) scenario (2006–2100). For this study, we conducted our analysis of the BCP using averaged model outputs over two time windows, the “reference” (1986–2005) and “future” (2091–2100).

2.2 Remote-sensing products

The five PISCES configurations were assessed for the reference period against an ensemble of satellite-derived NPP and C_exp products, averaged over 1998–2005. This period was chosen to align with the model simulations while accommodating the earliest availability of satellite observations. Although the averaging periods differ between the model outputs and remote-sensing products, the use of multi-year means reduces the influence of interannual variability and allows for a consistent comparison of large-scale patterns. Model outputs and remote-sensing products were compared globally as well as across RECCAP2 biomes (Fig. S1 in the Supplement) by ocean basin (Fay and McKinley, 2014; Doney et al., 2024).

Several algorithms have been described to estimate primary production based on ocean-colour remote-sensing observations, each differing in complexity and formulation and exhibiting substantial variability across regions such as the subtropical gyres and the Southern Ocean (Ryan-Keogh et al., 2023a; Westberry et al., 2023). Common algorithms include the: Vertically Generalised Production Model (VPGM), Carbon-based Production Model (CbPM), and Carbon Absorption Fluorescence Euphotic-resolving (CAFE) model. This study used six different NPP datasets (Fig. S2), all of which are derived from ESA Ocean Colour Climate Change Initiative (OC-CCI) remotely-sensed chlorophyll concentrations (Sathyendranath et al., 2019). Five were from the regular 25 km gridded NPP datasets at 8 d resolution of Ryan-Keogh et al. (2023a), who considered the Eppley-VGPM (Eppley, 1972), Behrenfeld-VGPM (Behrenfeld and Falkowski, 1997), Behrenfeld-CbPM (Behrenfeld et al., 2005), Westberry-CbPM (Westberry et al., 2008), and the Silsbe-CAFE model (Silsbe et al., 2016). The sixth NPP dataset was the 9 km monthly resolution product from the Biological Pump and Carbon Exchange Processes (BICEP) project (Kulk et al., 2020, 2021) which used a modified version of the base algorithm of Longhurst et al. (1995).

Similar to NPP, a number of algorithms exist to compute C_exp from remotely-sensed observations (Jönsson et al., 2023). Export production can be related to NPP through a simple scaling factor (Eppley and Peterson, 1979) or include multiple relationships, such as mixed-layer depth (MLD), chlorophyll and SST (Jönsson et al., 2023). For simplicity, we restricted our analysis to export algorithms that estimate C_exp from empirical relationships between NPP and SST. We applied five commonly used formulations (Fig. S3 and Table S1 in the Supplement): Laws et al. (2000), Henson et al. (2011), two equations from Laws et al. (2011), and Li and Cassar (2016). These algorithms have been evaluated against in situ observations (e.g. Dunne et al., 2005; Mouw et al., 2016; Bisson et al., 2018), with reported skill scores of $R^{\mathrm{2}}=\mathrm{0.64}-\mathrm{0.70}$ (Li and Cassar, 2016; Jönsson et al., 2023). It is important to note that they estimate export at different depth horizons such as Henson et al. (2011) at 100 m and Li and Cassar (2016) at the climatological MLD, meaning biases in subtropical gyres and polar regions may differ slightly across algorithms. We used SST from the 0.05° Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) product, which is a Level 4 product that combines multi-sensor satellite and in situ data (Good et al., 2020). From the daily fields for 1998–2005, monthly SST averages were computed for OSTIA. Thereafter, C_exp was computed using the five algorithms applied to the monthly SST and six NPP fields (interpolated to the 2° model grid), resulting in 30 individual C_exp estimates. For each algorithm, outputs from the six datasets were then averaged to produce a single representative C_exp field, yielding five ensemble algorithm outputs (Fig. S3).

3 Results

3.1 Model and remote-sensing variability

At both the global (Table 2) and biome scales (Fig. 1), differences in the complexity of phytoplankton processes and parameterisations (Table 1) among the five PISCES configurations leads to substantial intramodel variability in both the reference and future states (Figs. 2 and 3), particularly for NPP. For the reference period (Table 2), the coefficient of variation (CV) for global NPP across the five configurations is ∼20 %, indicating substantial divergence among models.

Figure 1 Model and remote-sensing (RMT) estimates of (a) NPP and (b) C_exp, integrated over each RECCAP2 biome (refer to Fig. S1). Black bars indicate ±1 standard deviation across the remote-sensing ensemble.

The two Monod-quota configurations, PST and PSF, simulate 15–21 Pg C yr⁻¹ less NPP than the purely Quota-based configurations (P5Z, P6Z and P5M). In contrast, the relative range of C_exp is much lower (CV≈5 %), showing export production is far more consistent across configurations and reflecting the fact that no direct adjustments were made to the equations or parameterisations explicitly controlling carbon export across the different configurations. Variability arises as an indirect outcome, emerging primarily from upstream differences in primary productivity patterns, which cascade down to influence export. Collectively, this results in the Monod-quota configurations resolving a more efficient BCP than the purely Quota-based ones. However, despite these intramodel differences, all the PISCES configurations largely fall within the range of remote-sensing estimates for both NPP and C_exp (Table 2). Compared to the CMIP5 ensemble, NPP for the PISCES configurations lies within the range of variability, whereas C_exp exceeds the upper bound.

Table 2 Reference and future projections of NPP, C_exp, and e-ratio for the five PISCES configurations. C_exp is defined as the mass of POC sinking through the 100 m depth horizon, consistent with CMIP conventions (Palevsky and Doney, 2018). The e-ratio (C_exp/NPP) reflects the efficiency of the BCP as governed by particle formation and sinking processes. NPP is integrated over the upper 100 m when calculating the e-ratio. NPP and C_exp are reported as globally integrated values (Pg C yr⁻¹), and the e-ratio as an area-weighted global mean. Relative changes (%) from the reference to future period are shown, alongside CMIP5 estimates from Fu et al. (2016) and remote-sensing (RMT) values from this study's ensemble and Doney et al. (2024).

At the biome scale (Fig. 1), no single configuration consistently reproduces all regions, with intramodel variability in NPP (Fig. 1a) and C_exp (Fig. 1b) varying considerably. A similar state of variability exists for remote-sensing derived estimates. The Southern Ocean shows the strongest intramodel variability (CV≥30 %; Fig. S4). P5Z and P6Z yield the highest NPP and C_exp, while the Monod-quota models give the lowest. P5M stands apart from both groups, with the inclusion of the manganese cycle dampening phytoplankton growth in this region, resulting in lower NPP than the other Quota-based configurations. This reduction in productivity propagates downstream, ultimately leading to lower C_exp.

In the equatorial and subtropical gyres, the Quota-based configurations simulate higher NPP, yet the simpler Monod-quota models produce slightly higher C_exp, indicating that greater primary production does not directly translate into greater export in these regions (see also Fig. 2). The inclusion of a diazotroph PFT in P6Z introduces regional differences relative to P5Z within these biomes, most notably in the South Atlantic and North Pacific subtropical gyres; however, C_exp remains broadly similar between the two configurations (Fig. 1b). The absence of a uniform NPP response to the added diazotroph group reflects basin-specific nutrient limitation regimes, particularly differences in Fe and P availability that constrain diazotroph growth across ocean basins (Wrightson et al., 2022). For PSF, the low-iron biological parameterisation results in slightly higher NPP in the equatorial biomes and subtropical gyres, but marginally lower productivity in the Southern Ocean compared to PST. This occurs because PSF prescribes smaller Fe : C quotas for phytoplankton, enabling growth under lower iron supply. As a result, productivity in low latitude biomes responds more strongly to assumptions about biological iron uptake, where Fe and N limitation play a larger role, whereas productivity in the Southern Ocean is less affected due to its already pervasive iron limitation (Tagliabue et al., 2020).

Figure 2 Top row shows zonally integrated (a) NPP and (b) C_exp, and the area-weighted mean (c) e-ratio for the reference period. The lower row (d–f) shows the respective changes in the future relative to the reference period.

3.2 Future projections of the BCP

From Table 2, all configurations project future global declines in both NPP and C_exp by 2100, with NPP decreasing by 6.30±3.00 % and C_exp by 12.94±1.65 %; however, the e-ratio remains largely unchanged. For NPP and C_exp, the largest absolute and relative declines are observed in the Quota-based configurations but all configurations fall within the variability of the CMIP5 ensemble (Table 2).

Figure 2 expounds on the biome scale assessment in Fig. 1, presenting zonally integrated NPP and C_exp, and area-weighted e-ratios for the reference period, along with their projected future changes. For the reference period (Fig. 2a–c), all configurations reproduce the expected pattern of elevated NPP and C_exp in the equatorial region and at high latitudes, including the Southern Ocean and the subpolar to polar regions of the Northern Hemisphere, with lower values in the subtropical gyres. Differences arise mainly in magnitude, with the Quota-based configurations simulating higher NPP than the simpler Monod-quota ones (Fig. 2a); however, this does not translate into uniformly higher C_exp (Fig. 2b), resulting in a less efficient BCP (Fig. 2c). Intramodel variability in C_exp is less than NPP, with PST and PSF producing slightly higher export production between ∼20° N/S but notably lower C_exp than the Quota-based configurations in the Southern Ocean (40–75° S).

For future projections (Fig. 2d–f), all configurations simulate declining NPP between 30° N/S (Fig. 2d). Beyond this band, the models diverge. Around 45° S, PST and PSF show an increase in NPP of ∼25 %, whereas the Quota-based configurations show only a 4 %–5 % rise. North of 30° N, the Monod-quota configurations exhibit increasing NPP (from $\sim \mathrm{1}\mathit{\%}-\mathrm{20}\mathit{\%}$ by 60° N), while the Quota-based configurations instead show a sharp decline of 15 %–25 % at 30° N, which gradually weakens poleward. All configurations eventually project higher NPP in the northern polar region (Fig. 2d). In contrast, future changes in C_exp (Fig. 2e) follow a broadly consistent pattern across configurations, differing mainly in magnitude rather than trajectory. All configurations show declining export production over most latitudes, with a regional exception in the Southern Ocean, where C_exp increases. This increase occurs farther north (∼45° S) in PST and PSF, reflecting their distinct NPP response, while the Quota-based models show a similar feature centred farther south. North of 30° N, all configurations show declining C_exp, though the strength of the reduction varies Collectively, the e-ratio remains relatively stable across latitudes, with no clear model-specific pattern between 30° N/S (Fig. 2f). All models show a 10 %–20 % decline in the e-ratio near 45° S, an increase of $\sim \mathrm{5}\mathit{\%}-\mathrm{10}\mathit{\%}$ south of 60° S, and a consistent decrease north of 30° N (Fig. 2f).

Intramodel variability in NPP, C_exp, and e-ratio across the reference and future periods reflects both the direct effects of differing phytoplankton parameterisations (Table 1) and the cascading, top-down consequences of these choices. To interpret the latitudinal patterns in Fig. 2, it is necessary to examine how phytoplankton biomass responds within each configuration. Accordingly, Fig. 3 shows the reference and relative future changes in total phytoplankton biomass, grouped into small (SML) and large (LRG) phytoplankton to accommodate differences in PFT structure across configurations.

Figure 3 Top row shows zonally integrated (a) Small (SML) and (b) Large (LRG) phytoplankton biomass, and the area-weighted mean (c) LRG/SML ratio for the reference period. Across the different PISCES configurations, PST and PSF have two PFTs, whereas the Quota-based configurations incorporate three (P5Z and P5M) and four (P6Z). To enable a comparison of phytoplankton biomasses, PFTs were grouped into “small” (nano-, pico- and N-fixers) and “large” (diatoms) categories. The lower row (d–f) shows the respective changes in the future relative to the reference period.

For the reference period (Fig. 3a–c), the spatial patterns and intramodel differences in small (Fig. 3a) and large (Fig. 3b) phytoplankton biomass closely parallel those of NPP (Fig. 2a) and C_exp (Fig. 2b). These patterns are expected and stem from the differing ecological functions and environmental sensitivities of small versus large phytoplankton in the model. Small phytoplankton dominate most of the ocean (LRG/SML<1; Fig. 3c) but contribute little to export because they remain in the surface layer, whereas large phytoplankton form larger, faster-sinking particles that disproportionately drive carbon export (Aumont et al., 2015). A striking feature is the substantially higher biomass of large phytoplankton (Fig. 3b) in P5M relative to the other configurations, particularly in the equatorial region and north of 30° N. Because P5M was specifically tuned to represent Mn and other co-limiting trace-metal effects in the Southern Ocean, its divergence elsewhere likely reflects model-specific tuning that produces unrealistic behaviour outside that region. This may also explain why P5M yields slightly higher global C_exp (Table 2) than the other Quota-based configurations.

Responding to future changes in marine environmental conditions, small and large phytoplankton show global declines of 7.88±2.97 % and 14.24±1.82 %, respectively by 2100. Across configurations, the latitudinal patterns of future changes in small (Fig. 3d) and large (Fig. 3e) phytoplankton biomass broadly follow those of NPP and C_exp (Fig. 2d and e). Intramodel divergence is most pronounced for small phytoplankton (Fig. 3d). South of 30° S, PST and PSF show increasing small phytoplankton biomass toward higher latitudes, whereas the Quota-based configurations simulate the opposite trend. Around ∼40° S this contrast is strongest, marking the clearest separation between the two model frameworks. North of 30° N, the configurations also diverge, with opposing latitudinal trends between 30–60° N. These differences likely stem from the contrasting phytoplankton complexity (i.e. the number of PFTs) and growth formulations across the Monod-quota and Quota-based configurations, which strongly shape how small phytoplankton, and therefore NPP, respond to future environmental change. This top-down interpretation is further supported by the broadly consistent spatial patterns and projected changes in N and P inventories across configurations (Fig. S5).

Large phytoplankton show comparatively little intramodel variability (Fig. 3e), with all configurations exhibiting similar spatial patterns and future changes, differing mainly in magnitude. The largest deviations occur again in P5M outside the Southern Ocean, reflecting its anomalously high large phytoplankton biomass in the reference state (Fig. 3b). Overall, the LRG/SML ratio is projected to decline globally by 6.83±3.33 %, with the largest absolute and relative changes occurring for PST and PSF, especially near ∼45 and 60° S (Fig. 3f). A declining ratio indicates a shift towards proportionally more small phytoplankton at the global scale, although regional exceptions occur, most notably in parts of the Southern Ocean, where large phytoplankton become locally more dominant (Fig. 3f). Nevertheless, the concurrent global declines in both small and large phytoplankton biomass ultimately reduces future NPP and C_exp.

4 Discussion

4.1 Differences in phytoplankton parameterisations driving intramodel variability

Our results show that the more complex Quota-based configurations, which each include some combination of additional PFTs, a more detailed growth formulation, and expanded nutrient cycling, simulate a less efficient BCP than the simpler Monod-quota configurations for both the reference and future periods. For the reference period, the added complexity across the five PISCES configurations produces notable differences at both global (Table 2) and biome (Fig. 1) scales, most clearly in NPP. Intramodel variability in C_exp is smaller and appears mainly as an indirect consequence of differences in NPP rather than from explicit changes to export related parameterisations. Despite these differences, all configurations reproduce broadly similar large-scale patterns in NPP, C_exp, and the distribution of small and large phytoplankton for the reference period, with intramodel differences arising in the magnitude of these respective fields.

The more complex Quota-based configurations simulate higher global NPP than the Monod-quota configurations, but this increase does not translate uniformly into higher C_exp because the additional productivity is carried by small phytoplankton, which enhance upper-ocean biomass but contribute little to export. From Fig. 2b, export production is greater in PST and PSF between ∼20° N/S, whereas outside these latitudes, particularly in the Southern Ocean, C_exp is highest in the Quota-based models. These contrasting, and in some regions contradictory, patterns in NPP and C_exp behaviour for the reference period between the configurations likely arise from top-down controls linked to the greater PFT complexity in the Quota-based configurations and its influence on zooplankton dynamics. These interactions also likely underpin the intramodel differences seen in the future projections of NPP and the BCP.

Between the standard model (Aumont et al., 2015) and PISCES-QUOTA, the two biggest updates are the representation of phytoplankton growth processes to use a quota-formulation and the addition of a dedicated picophytoplankton PFT (Kwiatkowski et al., 2018). Furthermore, the zooplankton feeding parameterisations are also adjusted. In the standard PST and PSF configurations, mesozooplankton feed on nanophytoplankton and diatoms, generating large sinking particles via faecal pellet production. In PISCES-QUOTA, however, picophytoplankton are grazed exclusively by microzooplankton rather than mesozooplankton. From this study, it is not possible to fully disentangle the respective contributions of the quota-based growth formulation and the inclusion of additional PFTs to variability in NPP and C_exp within the Quota-based configurations. The higher NPP simulated by these configurations (Table 2) is partly attributable to their explicit treatment of nitrogen assimilation costs as flexible stoichiometry reduces respiratory losses when ammonium (${\mathrm{NH}}_{\mathrm{4}}^{+}$) is available. This then yields higher effective growth rates than in the Monod-quota configurations, particularly in low latitudes and during summer at high latitudes. However, NPP differences between the Monod-quota and Quota-based configurations remain modest in low latitudes (Fig. 2a), whereas in polar regions, a clearer association emerges between elevated small phytoplankton biomass and enhanced NPP in the Quota-based configurations. Kwiatkowski et al. (2018) concluded that accounting for variable phytoplankton stoichiometry has a limited impact on the global carbon cycle such that the added PFT complexity is likely the primary parameterisation driving the divergence in behaviour between the Monod-quota and Quota-based configurations in our study.

The competition index (CI; Fig. 4a) shows that the single small nanophytoplankton PFT in the Monod-quota configurations outcompetes the nanophytoplankton PFT in the Quota-based configurations across the equatorial and lower-latitude regions, although it is consistently outcompeted by the picophytoplankton group across all regions (Fig. 4b). Splitting the small phytoplankton PFT of the standard model into two distinct PFTs in PISCES-QUOTA allows for greater ecosystem complexity and differentiation across environmental niches. For the open-ocean regions of the lower latitudes (0–30° N/S), picophytoplankton dominate, comprising 50 %–70 % of total phytoplankton biomass for the reference period (Fig. S6). Consequently, relative to the Monod-quota configurations, nanophytoplankton biomass is smaller in the Quota-based configurations, but this reduction is offset by picophytoplankton biomass, leading to similar overall NPP in the region. However, because mesozooplankton feed exclusively on nanophytoplankton and diatoms, the reduced nanophytoplankton biomass in the Quota-based configurations limits their food supply. This leads to lower mesozooplankton grazing rates (Fig. S7d) and a reduced production of large sinking particles. This mechanistically explains why export production is higher in the equatorial region in PST and PSF compared with the Quota-based configurations (Fig. 2b).

Figure 4 The competition index (CI) quantifies where the standard nanophytoplankton PFT (mean of PST and PSF) outperforms the nano- and picophytoplankton PFTs of the Quota-based configurations. CI is computed from the mean realised growth rates (μ) over the upper 100 m and defined as $\text{CI}=({\mathit{\mu }}_{\mathrm{1}}-{\mathit{\mu }}_{\mathrm{2}})/\left(\right|{\mathit{\mu }}_{\mathrm{1}}|+|{\mathit{\mu }}_{\mathrm{2}}\left|\right)$, where μ₁ is the mean Monod-quota nano-PFT and μ₂ is either the nano- or picophytoplankton PFT from the Quota-based configurations. Positive values showcase regions where the nanophytoplankton PFT of the Monod-quota configurations outcompetes the nano- and/or picophytoplankton PFT of the Quota-based models. Panels show CI for the reference (a, b) period and future shifts (c, d).

Outside of 30° N/S, the CI shows a shift in competitive behaviour, with the nanophytoplankton PFT for the Quota-based configurations now outcompeting the single nanophytoplankton PFT of the Monod-quota models (Fig. 4a). Because picophytoplankton already outcompete the Monod-quota nano- PFT everywhere (Fig. 4b), this means that small phytoplankton as a whole can proliferate more effectively in the Quota-based configurations beyond the tropics. This explains the higher small phytoplankton biomass (Fig. 3a) and elevated NPP (Fig. 2a) in these regions. This pattern is strongest in the Southern Ocean, where the Quota-based configurations show substantially higher NPP and export production compared to PST and PSF. The elevated small phytoplankton biomass fuels enhanced micro- and mesozooplankton grazing rates (Fig. S7c and d), enhancing export production in the Southern Ocean (Fig. 2b) and this mechanism is consistent with the findings of Laufkötter et al. (2013).

The mechanistic patterns identified in the reference period provide a critical foundation for understanding why certain regions show stronger divergence than others in the future projections of NPP and the BCP across the PISCES configurations. By 2100, our results show global NPP and C_exp are projected to decline by 6.30±3.00 % and 12.94±1.65 %, respectively, accompanied by a shift towards dominance of small phytoplankton species. These findings align with previous modelling studies (Bindoff et al., 2019), and fall within the variability of the CMIP5 ensemble (Bopp et al., 2013; Fu et al., 2016). The comparable magnitude of variability in NPP and C_exp (Table 2) indicates that differences in parameterisations among the selected PISCES configurations can generate a spread in results of similar order to that found across CMIP5 models (Séférian et al., 2020). While this does not imply that the present ensemble captures the full diversity of CMIP model structural differences, it nevertheless highlights the sensitivity of biogeochemical outputs to relatively subtle differences in the representation of phytoplankton growth processes and ecosystem complexity which do contribute to intermodel variability.

Global future reductions in NPP and C_exp likely reflect reduced nutrient supply to the surface associated with enhanced stratification and increased interior residence times (Dunne, 2023); however, the e-ratio shows little to no significant change for the five PISCES configurations, indicating that BCP efficiency is maintained relative to the reference period. In Fig. 4c and d, the CI shows little change between the reference and future periods, indicating that the relative competitive balance among PFTs remains largely unchanged. Consequently, the Quota-based configurations continue to sustain higher absolute NPP in the future (Table 2) than the Monod-quota configurations. However, they also maintain lower C_exp, reflecting persistent regional differences in phytoplankton-zooplankton interactions (Fig. S7) that stem from the contrasting levels of PFT complexity, particularly within the small phytoplankton community, across the PISCES configurations.

Intramodel variability is greatest between 30–60° N/S for future projections of NPP and C_exp. Small phytoplankton biomass increases from 30° S to a peak near ∼40° S in the Monod-quota configurations before weakening toward 60° S, whereas the Quota-based configurations show the opposite latitudinal pattern, with both model families converging at ∼60° S (Fig. 3d). NPP responds accordingly, though the impact on C_exp is more muted. Under future changes in marine environmental conditions, such as warmer SSTs and enhanced stratification, the simpler Monod-quota models project an increase in small phytoplankton biomass, whereas the Quota-based configurations project a decline in this region. This divergence arises from the added phytoplankton complexity in the Quota-based configurations. By explicitly representing picophytoplankton, P5Z, P6Z, and P5M resolve a broader range of ecological niches. As marine environmental conditions shift, the response of picophytoplankton dominates the overall behaviour of the small phytoplankton community, leading to a decrease in biomass (Fig. 3d). In contrast, the single nanophytoplankton PFT in the Monod-quota models, being more generalised and arguably over-parameterised, responds differently and shows an increase for 30–60° S. A similar pattern emerges in the Northern Hemisphere, where both model families project an overall decline in small phytoplankton biomass but the latitude–dependent trends differ (Fig. 3d). Again, this reflects the distinct niche sensitivities of the picophytoplankton PFT in the Quota-based configurations versus the more simplified single nanophytoplankton PFT present in the Monod-quota configurations.

4.2 Complexity and the persistence of uncertainty

For both the reference and future periods, our results show that differences in biogeochemical parameterisations constitute a major axis of divergence across the five PISCES configurations. Importantly, not all parameterisations exert influence at the same scale. Some, such as the introduction of the picophytoplankton PFT in PISCES-QUOTA, generates substantial shifts in global NPP, which in turn leads to moderate but consistent differences in C_exp, for both the reference and future periods relative to the Monod-quota configurations. Other parameterisations exert more limited effects, such as low-iron parameterisation in PSF, added diazotroph PFT in P6Z, and explicit manganese cycling for P5M, impacting phytoplankton growth rates or expanding phytoplankton complexity in ways that manifest significant differences in NPP and C_exp behaviour at the biome scale rather than globally. Together, our findings highlight that certain parameterisations fundamentally reshape model behaviour, whereas others provide more subtle or regionally confined refinements.

Added complexity in biogeochemical models is intended to improve the representation of marine ecosystem functioning. However, the current uncertainty ranges in remote-sensing estimates of global NPP and C_exp remain extremely large at 46.3 % and 109.8 %, respectively (Doney et al., 2024). As shown in Fig. 1, all five PISCES configurations fall well within the spectrum of remote-sensing estimates at the biome scale such that understanding whether the added complexity is improving representation of the marine system is difficult to discern. Although not the focus of the study, a skill assessment of the different configurations relative to the ensemble mean of the various remote-sensing NPP and C_exp datasets showed increased complexity may improve global NPP magnitude without improving, and possibly worsening, spatial pattern fidelity, while for C_exp, there is little difference in skill across configurations (Fig. S8). This raises an important point, one that is actively discussed within the modelling community, determining the optimal balance of model complexity, and the most effective combination of parameterisations, when utilising biogeochemical models for future projections of the carbon cycle and the BCP.

One may argue that more complex biogeochemical models permit more sources of error and degrees of freedom; however, Flynn (2010) makes the bold statement that models that ignore key biogeochemical processes are inherently dysfunctional. Especially for future projections beyond the extent of observations as the absence of process driven biogeochemical feedback will lead to compensating errors and large future uncertainties. Furthermore, one also needs to consider the computational expenses that accompany greater complexity as PISCES-QUOTA is far more costly to operate than the simpler standard model (Kwiatkowski et al., 2018). These quid pro quos preclude a single “best” configuration, but they illuminate the key considerations when choosing among PISCES configurations of differing complexity. Finally, we tested these different configurations under the very strong climate forcing of the RCP8.5 scenario. Results may be different under weaker levels of global warming currently expected or strongly mitigated scenarios.

5 Conclusion

In this study, we used five distinct PISCES configurations of varying complexity that differ in their number of PFTs, growth formulations, imposed nutrient limitations, and the representation of nutrient cycling. We showed that these parameterisation differences alone can generate substantial divergence and intramodel variability under the high emissions RCP8.5 scenario in projections of the BCP, with NPP responding more sensitively than C_exp. Yet, despite these contrasts, all configurations produce for the reference period NPP and C_exp magnitudes that fall within the broad range of remote-sensing estimates, making it difficult to assess whether added complexity unequivocally improves model realism. Moreover, the spread in future NPP and C_exp projections across the PISCES configurations is comparable to that seen in the CMIP5 ensemble, suggesting that relatively subtle differences in phytoplankton processes and their parameterisation may underpin a substantial fraction of the intermodel variability in ESM projections.

The introduction of a picophytoplankton PFT in PISCES-QUOTA, and its presence in the other Quota-based configurations, emerges as one of the most influential parameterisation choices, producing higher global NPP yet slightly lower C_exp, and driving opposing future NPP responses and latitudinal trends in C_exp between 30–60° N/S compared with the Monod-quota configurations. Other parameterisations, such as the low-iron scheme in PSF, the added diazotroph PFT in P6Z, and explicit manganese cycling in P5M, exert more modest, regionally confined effects, influencing NPP and C_exp primarily at biome scales rather than driving large-scale divergence in model behaviour. Finding that differing parameterisations produce contrasting model outputs is unsurprising, but the nature of our study allowed us to identify which parameterisations generate substantial intramodel variability and which exert only minimal influence. The discrete number of configurations inevitably constrains our conclusions to parameterisations centred on phytoplankton processes. Nevertheless, within the PISCES and broader ESM modelling community, this study represents a novel contribution, highlighting how different configurations of a biogeochemical model behave in projecting contemporary and future states of NPP and C_exp.