The BFM was coded under the conventions of the open-source, Fortran-based Framework for Aquatic Biogeochemical Models (Bruggeman and Bolding, 2014) as described in Álvarez et al. (2023) and now this biogeochemical model, called FABM-BFM hereinafter, is integrated in the distributed version of FABM (Bruggeman et al., 2024). FABM enables complex BGC models to be developed as sets of stand-alone, process-specific modules. These are combined at runtime through coupling links to create a custom-made biogeochemical/ecosystem model. One feature of the modularity of FABM-BFM is that the number of plankton types and their coupling links can be defined at runtime without changing or recompiling the code, by adjusting the fabm.yaml configuration file. The whole structure of the ecosystem used in this work (number of plankton types, size structure, allometric relationships, prey-predator relationships) was defined with the characteristics described in the following sections. The script that generates the fabm.yaml is available in the running repository (see Code Availability Section).

To demonstrate the applicability and how FABM-BFM can be used to develop diversity-capable models, we present the application to the size-structured plankton community model that resolves 60 phytoplankton and 37 zooplankton size classes (BFM98) in a one-dimensional pelagic ocean setting representing an open-ocean ecosystem in the NW Mediterranean Sea, the Bouée pour l’acquisition de Séries Optiques à Long Terme (BOUSSOLE) observatory. The site is visited monthly by oceanographic cruises to collect optical and BGC variables (Golbol et al., 2000). Seawater is generally collected at 12 discrete depths within the 0–400 m water column (5, 10, 20, 30, 40, 50, 60, 70, 80, 150, 200, 400 m). Sampling is performed with a SeaBird SBE911+ CTD/rosette system equipped with pressure (Digiquartz Paroscientific), temperature (SBE3+) and conductivity (SBE4) sensors and 12-L Niskin bottles, from which water samples are collected for subsequent absorption measurements and high-performance liquid chromatography (HPLC) analyses. The HPLC procedure used to measure phytoplankton pigment concentrations is fully described in Ras et al. (2008). The fractions of 3 phytoplankton size classes were calculated following the diagnostic pigments method (DiCicco et al., 2017). Satellite measurements were extracted from the CMEMS daily products of phytoplankton absorption and light extinction coefficient (OCEANCOLOUR_MED_BGC_L3_MY_009_143) averaging a 20x20km box around the BOUSSOLE site (E.U. Copernicus Marine Service Information, 2022).

Forced with climatologically averaged atmospheric input and initialization, the model settles into a repeated seasonal cycle after approximately 2–3 years of integration. The CLIM simulation does not lose diversity when the biomass of each plankton type is relaxed to their initial concentration with a time constant of 5 years ( Supplementary Material 2 ), which in practice means that all plankton types are always present in the system and the net inflow of biomass to the system is less than 1 ng C L^-1 d^-1. Therefore, from here thereinafter, a plankton type was considered present in the community if the biomass is larger than 0.01 mg C m^-3. Despite the fact that BFM98 contains 386 state variables compared to the 55 of BFM9, the integration and initialization times are less than three times longer ( Supplementary Material 1 ).

The annual evolution of real-date bulk biogeochemical properties simulated by BFM9 and BFM98 were compared to observations collected at the DYFAMED and BOUSSOLE observatories. Chl-a concentration (mg Chl-a m^-3) at surface peaks in spring and decreases in summer comparably to in situ measurements ( Figure 2a ). At the depth of the Chl-a maximum (DCM) the model tends to underestimate Chl-a at the beginning of the summer and overestimate it during summer and autumn. Integrated phytoplankton ( Figure 2b ) and bacterial biomass ( Figure 2c ) in the first 130 m of the water column show a clear alternance between phytoplankton and bacterial dominated communities. One remarkable difference between BFM9 and BFM98 is the larger concentrations of phytoplankton in the latter ( Figure 2b ), concurrent with less bacteria ( Figure 2c ) and less heterotrophic nanoflagellates ( Figure 2d ). All zooplankton integrated and mean biomasses tend to peak after spring in relatively good agreement with the scarce observational data ( Figures 2e, f ), but BFM98 tends to sustain larger biomasses of micro- and meso-zooplankton throughout the year.

Dissolved organic carbon concentrations simulated by the model tend to be slightly smaller than those observed but both model configurations are able to capture the accumulation of DOC during summer and autumn in surface (20 m) and the relative lower concentrations at 90m ( Figure 2g ). The ratio between particulate organic carbon to particulate nitrogen ratio observed at the site tends to be smaller than 6 during spring to become larger in summer, whereas BFM98 (B+P+NAP+Z⁽⁶⁾+Z⁽⁵⁾) tends to be more conservative around the Redfield 6 value and BFM9 be smaller ( Figure 2h ). Total annual net primary production integrated in the first 90 m simulated is 114.7 and 118.3 gC m^-2 by BFM9 and BFM98, respectively ( Figure 2i ). In both cases underestimated when compared to short-term incubations done in situ but comparable to estimates in the NW Mediterranean Sea (Siokou-Frangou et al., 2010).

The spectrally-resolved version of BFM resolves light transmission in the water column as it is absorbed and scattered by water, phytoplankton, CDOM and detritus. Figure 3a shows the vertical distribution (0 to 200 m) of total Chl-a concentration along an average year. The simulated Chl-a accumulates in spring in the first 30 m to later remain confined to deeper layers from 40 m at the beginning of the summer to around 55 m in autumn ( Figure 3a ), in a pattern observed at the site (Antoine et al., 2006) and in general in the open waters of the NW Mediterranean Sea (Lavigne et al., 2015). The inherent optical properties simulated in the proposed configuration were consistent with previous lower complexity versions of the model (Álvarez et al., 2023) with a_PH(450) following the distribution of Chl-a ( Figure 3b ) and a_CDOM(450) forming a sub-superficial maximum in summer ( Figure 3c ). The contribution of NAP to total non-water absorption [a_NAP(450)] is slightly larger during spring given a large presence of detritus (not shown). The annual evolution of simulated a_PH(450) ( Figure 3d ) and light extinction coefficient (k_D(490), Figure 3e ) was compared to observations collected remotely in the north west Mediterranean Sea, showing a relative good agreement. BFM98 tends to simulate slightly lower k_D(490) out of the spring bloom compared to BFM9.

The simulated vertical distribution (0 to 200 m) and temporal variation (January to December) of phytoplankton size classes according to Chl-a concentrations ( Figure 4 right panels) is consistent in shape and magnitude with observations derived from HPLC data ( Figure 4 left panels), which clearly reveal the dominance of nanoplankton in the site and the importance of diatoms and picoplankton in the bloom and in the summer DCM, respectively (see Álvarez et al. (2023) for the same analysis of PSCs at BOUSSOLE with BFM9). A visual comparison of the dominance of different phytoplankton size classes at the DCM throughout the year ( Figure 4 bottom panels) shows how the diatoms in the model seem to peak later than observed and nanoflagellates tend to dominate the early bloom whereas picoplankton dominates later in the year.

The number of phytoplankton and zooplankton types present in the community (functional richness) across different years indicates seasonal variations and inter-annual trends. On the y-axis, Figure 5 shows the functional richness of both the phytoplankton and zooplankton community. Each line represents the simulation of a different year from 2000 to 2022, color-coded accordingly. For phytoplankton, while most changes in richness at surface are driven by the abrupt increase in biomass of the spring bloom ( Figure 5a ), when the DCM is formed richness remains low ( Figure 5b ) despite relatively high values of biomass ( Figure 2a ). Richness of zooplankton types is less variable; at any given point in time only 14 to 20 types are present in the community both at surface ( Figure 5c ) and at the DCM ( Figure 5d ).

The CLIM simulations, where BFM98 was integrated forward with climatologically averaged atmospheric forcing during 100 repeated years, showed that after the first 2–3 years, no significant drift in the functional richness of phytoplankton ( Figure 6a ) nor zooplankton (not shown) was observed. The TREND simulation was identical to CLIM, except that it included a constant decrease in nutrient concentrations from 2005 to 2099 of -13.1% in phosphate and -19.8% in nitrate. Phytoplankton richness in TREND ( Figure 6b ) showed smaller seasonal amplitude by the end of the hundred years. To explore such differences between CLIM and TREND, we employed community turnover. Community turnover computes the difference in plankton types in two timepoints as the total number of types gained (NG) and/or lost (NL) with respect to the total number of types observed in both timepoints (NT). The month-to-month turnover (Righetti et al., 2019) computes the average turnover ((NG+NL)/NT) between all consecutive months of the year, therefore mainly capturing the seasonal differences in community composition. The seasonal turnover of the climatological run (CLIM) stabilizes after approximately one decade around 0.17 ( Figure 6c ), but it shows some variability around this value. The seasonal turnover of the TREND run is smaller and reaches 0.10 at the end of the 100 years, corroborating the smaller seasonal differences in phytoplankton community composition with decreasing nutrient availability, and, in practice, indicating that the model diversity is able to change in response to environmental forcing.

The ability to simulate realistic plankton diversity beyond the current state-of-the-art models that rely on functional types and aggregated taxa, alongside biogeochemistry can reveal previously unrecognized dynamics in plankton community interactions and their cascading effects on biogeochemical cycles. By explicitly representing a diverse plankton community, our model facilitates future research into the interplay between biogeochemical fluxes and biodiversity. Changes in plankton diversity may have profound implications for ecosystem functioning, as they can influence primary production and nutrient remineralization pathways given differential nutrient uptake strategies among coexisting phytoplankton groups (Kwon et al., 2022); alter the efficiency of the biological carbon pump, since size and other particle properties influence particle sinking and fluxes (Cael et al., 2021a); and influence and the energy transfer across trophic levels though shifts in predator-prey interactions and trophic dynamics (Serra-Pompei et al., 2022).

Our diversity-capable FABM-BFM model incorporates species selection through migration, and thanks to the species selection approach implemented formally through migration (Bruggeman and Kooijman, 2007), effectively representing eco-evolutionary processes that influence community dynamics. By defining a sufficiently complex and broad trait space, the model allows for adaptive shifts in community structure in response to environmental conditions, implicitly capturing how marine communities alter their functional structure, relocate, or adapt in response to environmental shifts. Investigating diversity-sustaining mechanisms is an active area of research and a diversity-capable model as this presented could be used as a laboratory to test those, including the role of trophic interactions (Edwards and Yool, 2000; Vallina et al., 2014), transport (Junger et al., 2023) or even explicit adaptation (Sauterey et al., 2015; Le Gland et al., 2021) could be further explored. The explicit representation of evolutionary dynamics is beyond the scope of this work, but could complement this approach in future studies, given that the role that evolution plays in marine microbial ecology and carbon cycling is intriguing (Ward et al., 2019). Whether it is implicit or explicitly, to incorporate the processes of adaptation into biogeochemical models is a pressing need as there is increasing observational evidence for significant regional regime shifts at the microbial level particularly in response to climate change (e.g Cael et al., 2021b). Exploring the evolutionary aspects of plankton diversity could offer important insights into how these communities may adapt to changing environmental conditions in the coming decades.

One of the fundamental challenges that diverse-capable ecosystem models face is the trade-off between increasing complexity and the resulting increase in uncertainty. As models incorporate more plankton types to capture the intricate dynamics of marine ecosystems, they can improve the representation of diverse species interactions, biogeochemical cycles, and ecosystem functions. However, it also introduces greater uncertainty, as more parameters need to be estimated and many interactions remain poorly understood or data-limited. Balancing this trade-off requires careful calibration and validation. A diversity capable module developed within a coupled optical-biogeochemical model permits the whole system to be constrained by multiple observational datasets, including those biogeochemical and optical, and furthermore incorporating novel sources of observational data, such as high-throughput plankton imaging (Pierella Karlusich et al., 2022) or molecular data (Levine et al., 2025).

One key challenge in modeling biodiversity through discrete plankton types is scalability. Like other trait-based models (Dutkiewicz et al., 2020), our model captures only a fraction of biodiversity, highlighting a fundamental limitation of our approach: the trait variability included in the model dictates its ability to simulate biodiversity responses to environmental changes. In the BFM98 configuration, the phytoplankton community composition is dynamic, but future model development should incorporate additional traits (such as temperature optima) to better represent ecosystem shifts under multifaceted environmental change. Such further refinement in resolving species-level interactions and fine-scale ecological processes could enhance the predictive power and accuracy of the proposed model in simulating future changes under varying climatic conditions.

Under the light of those challenges, the modular and scalable nature of FABM-BFM presents several advantages. Traditional models with discrete plankton types become computationally expensive as the number of state variables increases. However, the flexibility of FABM-BFM allows for dynamic ecosystem configurations at runtime, enabling more efficient computational scaling. The 1D versions of the FABM-BFM are relatively fast, making it feasible to conduct extensive analyses that require numerous simulations, such as sensitivity testing of model assumptions or data assimilation procedures. This also makes century-long experiments highly achievable and facilitates the transition to three-dimensional configurations. Additionally, the possibility to perform multi-model simulations is especially valuable for exploring biodiversity-related questions (Jones and Cheung, 2015). Such multi-model ensemble approach (facilitated by FABM) is an essential aspect to perform scenario simulations and associate uncertainties in model projections, that become more robust with respect to the potential differences among the specific ecological models.

By allowing the representation of multiple traits in a broad trait space, FABM-BFM can simulate biodiversity under climatic projections, and explore how predicted changes in diversity could impact ocean biogeochemistry, ocean-color and ecosystem stability or resilience. A key advantage of the proposed framework in this regard is that it is built on the BFM, which underpins the operational biogeochemical products provided by the Copernicus Marine Service (Coppini et al., 2023). Incorporating biodiversity into the BFM framework enhances its potential for future operational applications, particularly in the context of climate change.

Despite the advancements presented in this study, several limitations remain that future work should address. First, the model does not yet incorporate all relevant ecological traits, such as thermal norms (see Supplementary Material 1 ) that determine species’ temperature-dependent performance. Including these traits would improve the ability of the model to simulate species distributions and community responses to climate variability. Second, while our approach allows for the representation of several hundred plankton types, it remains a simplified functional-type framework. Furthermore, the model has not yet been systematically compared to more detailed species-level models or real-world plankton community observations. Assessing how well the current biodiversity representation captures observed ecological processes and community dynamics is a necessary step for further validation and refinement. Finally, the absence of lateral transport in the 1D configuration constrains the influence of connectivity, necessitating prescribed migration fluxes to approximate transport effects. A key question for future work is whether the transport processes in a fully three-dimensional implementation would be sufficient to maintain diversity, or if additional ecological processes must be incorporated to sustain realistic plankton community structures over time, such as more complex predator-prey dynamics or automatic differentiation of traits.