Chromophoric dissolved organic matter dynamics revealed through the optimization of an optical–biogeochemical model in the northwestern Mediterranean Sea

Chromophoric dissolved organic matter (CDOM) significantly contributes to the non-water absorption budget in the Mediterranean Sea. The absorption coefficient of CDOM, a_CDOM(λ), is measurable in situ and can be retrieved remotely, although ocean-colour algorithms do not distinguish it from the absorption of detritus. These observations can be used as indicators for the concentration of other relevant biogeochemical variables in the ocean, e.g. dissolved organic carbon. However, our ability to model the biogeochemical processes that determine CDOM concentrations is still limited. Here we propose a novel parameterization of the CDOM cycle that accounts for the interplay between the light- and nutrient-dependent dynamics of local CDOM production and degradation, as well as its vertical transport. The parameterization is included in a one-dimensional (1D) configuration of the Biogeochemical Flux Model (BFM), which is here coupled to the General Ocean Turbulence Model (GOTM) through the Framework for Aquatic Biogeochemical Models (FABM). Here the BFM is augmented with a bio-optical component that resolves spectrally the underwater light transmission. We run this new GOTM-(FABM)-BFM configuration to simulate the seasonal a_CDOM(λ) cycle at the deep-water site of the Bouée pour l'acquisition de Séries Optiques à Long Terme (BOUSSOLE) project in the northwestern Mediterranean Sea. Our results show that accounting for both nutrient and light dependence of CDOM production improves the simulation of the seasonal and vertical dynamics of a_CDOM(λ), including a subsurface maximum that forms in spring and progressively intensifies in summer. Furthermore, the model consistently reproduces the higher-than-average concentrations of CDOM per unit chlorophyll concentration observed at BOUSSOLE. The configuration, outputs, and sensitivity analyses from this 1D model application will be instrumental for future applications of BFM to the entire Mediterranean Sea in a three-dimensional configuration.

1 Introduction

The inventory and composition of dissolved organic matter (DOM) is of utmost importance to marine ecosystems as it is the energy base and carbon source for heterotrophic life in the ocean (Azam et al., 1983). Furthermore, as the largest pool of reduced carbon in the oceans, DOM plays an important role in the global carbon cycle (Legendre et al., 2015). Although most of the dissolved exudates that form the DOM are non-light-absorbing (Mühlenbruch et al., 2018), part of the pool absorbs light mainly in the ultraviolet (UV) and blue spectral range of the electromagnetic radiation. This portion of the DOM is referred to as chromophoric dissolved organic matter (CDOM, Jerlov, 1951; Sieburth and Jensen, 1968). CDOM absorbs light depending on both the CDOM concentration and its mass-specific spectral absorption coefficient. a_CDOM(λ) is a measurable proxy for fundamental processes regulating CDOM dynamics. The high absorption of CDOM in the blue part of the electromagnetic spectrum affects the water leaving radiance, a radiometric quantity that can be retrieved by satellites, and interferes with estimates of chlorophyll a concentration (Chl a) from ocean-colour observations (Siegel et al., 2013). This is particularly important for the Mediterranean Sea because its surface waters (0–20 m) appear greener than the global ocean, although its phytoplankton concentration is typically lower (Chl a<0.2 mg m⁻³) (Morel and Gentili, 2009b; Claustre et al., 2002; Volpe et al., 2007). Besides the presence of CDOM, there are other possible factors influencing this optical behaviour, including the particular pigment ratios in the Mediterranean phytoplankton community (Organelli et al., 2011), the abundance of small coccolithophores (Gitelson et al., 1996), and the influence of Saharan dust (Claustre et al., 2002). However, the high contribution of CDOM was found to be a major factor since the CDOM concentration in the basin is about twice that in the Atlantic at the same latitudes (Morel and Gentili, 2009b).

a_CDOM(λ) can be measured in situ at selected locations and can be retrieved on a global scale from remote-sensing platforms, although remotely it is not distinguishable from the absorption of detritus. Remote-sensing platforms provide radiometric observations, which are used to compute a combined value for the absorption coefficient of CDOM and of detritus, a_DG(λ) (Werdell et al., 2018). On the other hand, spectrally resolved biogeochemical (BGC) models provide a_CDOM(λ) estimates by simulating CDOM concentrations and prescribing the optical properties of the pool. Such BGC–optical models could advance the understanding of the dynamics shaping a_CDOM(λ). However, they rely on the accurate representation of both the dynamics of CDOM cycling and its optical properties. CDOM cycling is essentially controlled by in situ biological production (Romera-Castillo et al., 2010), terrestrial and atmospheric sources and physical transport, microbial consumption (Nelson and Gauglitz, 2016; Legendre et al., 2015; Stedmon and Markager, 2005), and deep ocean circulation and/or vertical mixing (Coble, 2007), and it is photoreactive and efficiently destroyed in the upper layers of the water column by solar radiation (Mopper and Kieber, 2000). Most marine ecosystem models that explicitly resolve CDOM follow Bissett et al. (1999) (e.g. Xiu and Chai, 2014; Dutkiewicz et al., 2015), in which CDOM is assumed to have a source that is a fixed fraction of the dissolved primary production (dpp) of phytoplankton, is remineralized slowly, and is bleached under high-light conditions. These formulations seem to work well in open waters (Dutkiewicz et al., 2015) but not in the northwestern (NW) Mediterranean Sea, where CDOM is not a constant proportion of dpp, even when co-varying with phytoplankton Chl a (Lazzari et al., 2021b). The aim of this work is to decipher the processes that influence CDOM dynamics in the NW Mediterranean Sea by advancing and using a BGC–optical model in synergy with in situ observations.

A major challenge in modelling CDOM concentrations in the Mediterranean Sea concerns the origin of the CDOM pool in superficial waters – i.e. how large is the proportion of CDOM transported from allochthonous sources compared to CDOM produced locally by primary and secondary production (Santinelli et al., 2002; Copin-Montégut and Avril, 1993)? In the basin, the higher-than-average CDOM concentrations seem to be sustained by allochthonous sources as fluxes of DOM depositions from the atmosphere are 2–5 times larger in the Mediterranean Sea than those estimated for the global ocean, which explains the abundance of humic-like substances (Santinelli, 2015; Galletti et al., 2019). Therefore, accurate modelling of the balance between CDOM transport (atmospheric input, advection, and/or vertical replenishment from deeper stocks) and destruction at the surface (photolysis and bacterial consumption) is crucial to accurately represent the role of CDOM of allochthonous origin in shaping the dynamics of a_CDOM(λ) in the NW Mediterranean Sea.

As for locally produced CDOM, the challenge is that the mechanisms for CDOM production and destruction involve processes that are decoupled from phytoplankton dpp, such as photolysis, zooplankton sloppy feeding, and bacterial production and consumption. In the open waters of the NW Mediterranean Sea, far from terrestrial inputs, photobleaching and biological production and consumption are suggested to be important factors determining the seasonality of a_CDOM(λ) (Pérez et al., 2016; Xing et al., 2014; Organelli et al., 2014). a_CDOM(λ) and Chl a concentration roughly covary at the surface, but the causative mechanisms are different. While a_CDOM(λ) dynamics at the surface are mainly explained by vertical transport of CDOM from depth and by photochemical destruction, Chl a concentration at the surface is mainly driven by photoacclimation of phytoplankton and nutrient limitation. There is also a close coupling between subsurface a_CDOM(λ) and phytoplankton Chl a via microbial activities or interactions in the planktonic food web (Pérez et al., 2016; Xing et al., 2014; Organelli et al., 2014), although the relative contribution to CDOM production by phytoplankton and by bacteria is still controversial (Romera-Castillo et al., 2010). Although it is generally assumed that phytoplankton in open oceans are not a direct source of CDOM but rather a source of labile organic matter that is microbially transformed and subsequently produces CDOM (Nelson et al., 1998; DeGrandpre et al., 1996), observations in the NW Mediterranean Sea suggest direct CDOM production by phytoplankton (Oubelkheir et al., 2007, 2005; Xing et al., 2014). Therefore, to incorporate CDOM in an ecosystem model of the NW Mediterranean Sea, one must consider the composition of DOM produced by phytoplankton in addition to the cycling of dissolved matter and the microbial loop in the area.

We approached this challenge by using the Biogeochemical Flux Model (BFM, Vichi et al. 2007) in a one-dimensional (1D) configuration. The BFM was coupled to the General Ocean Turbulence Model (GOTM, Burchard et al., 1999), a 1D water column turbulent kinetic energy model, by using the Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014). In such a configuration, the model accounts for the vertically differentiated processes of photobleaching and transport from deep inventories. Local CDOM production includes the excretion of phytoplankton, which depends on light and nutrient availability, the bacterial production and consumption, and zooplankton activities. The version of BFM used here resolves the spectral light transmission underwater (Terzić et al., 2021; Lazzari et al., 2021a) and simulates the inherent optical properties (IOPs) of CDOM, phytoplankton, and organic detritus in the water column. The model output can be compared with a range of optical observations made in the field and remotely.

With this model configuration, we simulated a data-rich monitoring site in the NW Mediterranean. The Bouée pour l'acquisition de Séries Optiques à Long Terme (BOUSSOLE) observatory is located in the Ligurian Sea, about 32 nautical miles off the French coast (Antoine et al., 2006). In this area, the water depth varies between 2350 and 2500 m, and the currents are extremely slow. This peculiarity is due to its location near the centre of the cyclonic circulation that characterizes the Ligurian Sea, and it makes reasonable the approximation of a 1D configuration. The BOUSSOLE observatory consists of a mooring where a buoy, specifically designed for collecting radiometric and bio-optical quantities, collects optical data at a high temporal resolution. Oceanographic cruises visit the site monthly to collect a complementary dataset of algal pigment concentrations and IOPs, which include light absorption coefficients of phytoplankton, detritus, and CDOM (Antoine et al., 2006). Numerous studies combining optics and biochemistry have exploited these two complementary sets of observations collected at the BOUSSOLE site, e.g. fluorescence to infer Chl a (e.g. Bayle et al., 2015; Mignot et al., 2011) and phytoplankton community composition (Sauzède et al., 2015), to infer the size structure of the phytoplankton community from measured light spectra (Organelli et al., 2013), and to estimate community production from particle backscattering coefficients (Barbieux et al., 2022; Barnes and Antoine, 2014). One-dimensional models developed or employed at the BOUSSOLE site have focused separately on physics and biochemistry (e.g. D'Ortenzio et al., 2008; Ulses et al., 2016) or optics (Blum et al., 2012), but to our knowledge, no spectrally resolved BGC model has been used previously to simulate the site.

Methods that merge models and observations can make significant progress in reducing model uncertainty (e.g. Teruzzi et al., 2021; Ciavatta et al., 2018), and data-intensive processes such as model calibration and validation are of utmost importance to make model applications reliable (Arhonditsis and Brett, 2004). A major challenge when working with complex coupled hydrodynamic–BGC numerical models is the need to make an appropriate choice of model parameter values that ensure optimal performance in terms of reproducing observational data. Genetic algorithms are commonly used by BGC modellers to solve nonlinear optimization problems such as parameter estimation (e.g. Kriest et al., 2017; Kuhn and Fennel, 2019; Falls et al., 2022; Wang et al., 2020; Mattern and Edwards, 2017). In this work, we calibrated the model by using the Parallel Sensitivity Analysis and Calibration utility (ParSAC, Bruggeman and Bolding, 2020), which implements the Differential Evolution (DE) algorithm (Storn and Price, 1997). The DE algorithm is a stochastic, population-based, and direct-search parallel method which operates through genetic operations, namely mutation, crossover, and selection. Its objective is to eventually guide the population to the most likely values of model parameters based on the agreement of model output with observations.

With these three elements, (i.e. the 1D hydrodynamic–BGC–optical model GOTM-(FABM)-BFM, the observations collected at BOUSSOLE, and the optimization tool), here, we propose and optimize a novel formulation for the simulation of CDOM at the BOUSSOLE site. The specific objective of the present study is to improve the representation of the variability of a_CDOM(λ) within a coupled BGC–optical model of the ocean by (i) linking biodegradable CDOM produced by phytoplankton to both nutrient availability and light intensity and (ii) improving the representation of the sources of allochthonous CDOM. We optimized optical and BGC parameters to improve the representation of the observed variables at the BOUSSOLE site and investigated the model simulation of the unobserved variables (e.g. DOC, bacterial concentration) and bio-optical relationships between a_CDOM(λ) and Chl a and DOC. Based on the optimized model configuration, we conducted two numerical experiments to gain new insights into (i)  the BGC processes that determine the seasonal dynamics of CDOM at the surface and at the depth of the Chl a maximum (DCM) and (ii) the role of local production versus allochthonous inputs in the total pool of CDOM.

2 Methods

The model tested in the present work is a one-dimensional (depth) configuration of the BFM, extended here with a bio-optical component, as described in Sect. 2.1. The observational data, which are explained in detail in Sect. 2.2, come from the BOUSSOLE project. A test case of the BOUSSOLE site in the NW Mediterranean from the surface to 2400 m is used, and the setup of the site is described in Sect. 2.3. The optimization strategy is described in Sect. 2.4, which also outlines the details of the experiments.

2.1 Model description

The modelling framework consisted of four elements: (i) a hydrodynamic model that simulates the vertical mixing of BGC variables, (ii) a coupling layer that connects the hydrodynamic model to the BGC model, (iii) a BGC model that simulates the sink and source terms of BGC variables, and (iv) a bio-optical model that resolves the underwater field of spectral light based on the presence of optical constituents in the water column.

2.1.1 The hydrodynamic model GOTM and the framework FABM

The General Ocean Turbulence Model (GOTM, Burchard et al., 1999) is a one-dimensional water column built around an extensive library of turbulence closure models. GOTM is very configurable and simulates the vertical structure of the water column – notably saline, thermal and turbulence dynamics. Being a 1D model, horizontal gradients at the application sites need to be either negligible or parameterized. A key advantage of a 1D model is that it is feasible to make long simulations with high vertical resolution as the computational resources required are much smaller compared to full 3D models. The open-source, Fortran-based Framework for Aquatic Biogeochemical Models (FABM, Bruggeman and Bolding, 2014) provides a coupling layer that enables the flexible coupling of ecosystem processes into the GOTM. 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 BGC–ecosystem model. At each simulation time step, the BGC equations are applied to each layer, and the rates of sink and source terms are calculated at the current time and in the current space using local variables (e.g. light, temperature, and concentrations of state variables) provided by the GOTM. The rates pass to the hydrodynamic model via FABM, which handles numerical integration of the BGC processes and the transport of BGC substances (e.g. nutrients, dissolved and particulate organic matter) between the layers. Updated states are then passed back to the BGC model via FABM.

The implementation of the BFM in FABM comprises 54 state variables. These include representations of dissolved inorganic carbon (O^(DIC)); inorganic forms of nitrogen ($N^{\left({\mathrm{NO}}_{\mathrm{3}}\right)})$; ammonium ($N^{\left({\mathrm{NH}}_{\mathrm{4}}\right)})$; phosphorus ($N^{\left({\mathrm{PO}}_{\mathrm{4}}\right)})$ and silicate (N^(Si)); four phytoplankton types (P^(DIATOM), P^(NANO), P^(PICO), and P^(DINO)); heterotrophic bacteria (B) and four grazers (Z^(HETEN), Z^(MICRO), Z^(OMNI), and Z^(CARNI)); particulate organic matter (R^(DET)), three pools of dissolved organic matter differentiated by reactivity into labile (R⁽¹⁾), semi-labile (R⁽²⁾), and semi-refractory (R⁽³⁾); and CDOM differentiated by the same reactivities (X⁽¹⁾, X⁽²⁾ and X⁽³⁾). In the following, subscripts are appended to the symbols for the components to indicate the elemental constituent for which the state variable stands, including carbon (C), nitrogen (N), phosphorus (P), Chl a (only in phytoplankton), and silica (Si, only in P^(DIATOM)). In the spectrally resolved version of BFM used in this work, the transmission of light in the water column is resolved in 33 wavebands centred from 250 to 3700 nm, with 25 nm spectral resolution within the visible range (Terzić et al., 2021; Lazzari et al., 2021a). Vertically resolved irradiance is absorbed and scattered by water, phytoplankton, CDOM, and detritus. Although the code implementation involved a redesign of the BFM code into a FABM-compliant modular structure, the core of the overall conceptual model of the spectrally resolved BFM remains intact compared to previous applications (e.g. Lazzari et al., 2021b; Álvarez et al., 2022) and is described below. The sources for the code of GOTM, FABM, and the spectrally resolved BFM adapted to work under the FABM convention can be found in the “Code availability” section.

2.1.2 The BGC model BFM and the bio-optical component

A complete record of the partial differential equations for each state variable that conform to BFM can be found in several publications (Vichi et al., 2007a, b) along with the details of the spectrally resolved version (Terzić et al., 2021; Lazzari et al., 2021a). The following sections describe the main dynamics of the cycling of phytoplankton functional types (PFTs), DOC and CDOM, the optical treatment of constituents (PFTs, CDOM, and detritus), and the solution for light transmission. Following Vichi et al. (2007), the dynamical equations are written in “rate of change” form, where the right-hand side contains the terms representing significant processes for each living or non-living component. For each process, the state variable subject to change appears in front of a vertical bar, and the superscript and subscript after the bar are the three-letter acronym of the process represented and the state variable involved in the process, respectively.

Cycling of BGC variables: PFTs, DOC, and CDOM

Primary producers are divided into four types that roughly represent the functional spectrum of phytoplankton in marine systems: diatoms (P^(DIATOM)), nanoflagellates (P^(NANO)), picophytoplankton (P^(PICO)), and dinoflagellates (P^(DINO)). The variation in carbon concentration (mg C m⁻³ d⁻¹), being i of each of the PFTs, is calculated as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{d}{P}_{\text{C}}^{\left(i\right)}}{\mathrm{d}t}}=\left.{\displaystyle \frac{\mathrm{d}{P}_{\text{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\text{gpp}}{O^{\left(\text{DIC}\right)}}}-\left.{\displaystyle \frac{\mathrm{d}{P}_{\text{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\text{dpp}}{R_{\text{C}}^{\left(\mathrm{2}\right)}}}-\left.{\displaystyle \frac{\mathrm{d}{P}_{\text{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\text{rsp}}{O^{\left(\mathrm{DIC}\right)}}}\\ & -\sum _{j=\text{HETEN,MICRO,OMNI, CARNI}}\left.{\displaystyle \frac{\mathrm{d}{P}_{\text{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\text{grz}}{Z_{\text{C}}^{\left(j\right)}}};\\ & i=\text{DIATOM,NANO,PICO,DINO}.\end{array}\end{array}$$

Here, we briefly describe spectrally resolved gross (gpp) and dissolved primary production (dpp), whereas the exact formulations for the respiration (rsp) and grazing (grz) processes can be found in the BFM manual (Vichi et al., 2020). Gross primary production of each phytoplankton species (mg C m⁻³ d⁻¹) is computed as follows:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{gpp}}{O^{\left(\mathrm{DIC}\right)}}}=r_{{\mathrm{0}}_{\mathrm{P}}}^{\left(i\right)}\times f{t}_{\mathrm{P}}^{\left(i\right)}\left(T\right)\times f{n}_{\mathrm{P}}^{\left(i\right)}\left(\mathrm{nutrients}\right)\\ & \times \left(\mathrm{1}-\exp \left({\displaystyle \frac{-{\mathit{\varphi }}^{\left(i\right)}\times {\mathit{\theta }}^{\left(i\right)}\times {\int }_{\mathrm{387.5}}^{\mathrm{800}}{a}_{\mathrm{PS}}^{\ast \left(i\right)}\left(\mathit{\lambda }\right)\times E_{\mathrm{0}}\left(\mathit{\lambda }\right)d\mathit{\lambda }}{r_{{\mathrm{0}}_{\mathrm{P}}}^{\left(i\right)}}}\right)\right)\\ & \times {\mathrm{P}}_{\mathrm{C}}^{\left(i\right)},\end{array}\end{array}$$

where $r_{{\mathrm{0}}_{\text{P}}}$ is the maximum productivity rate (d⁻¹), regulated by dimensionless temperature- (ft_P(T)) and nutrient-dependent (fn_P(nutrients)) factors. Growth-limiting nutrients are $N^{\left({\mathrm{PO}}_{\mathrm{4}}\right)}$, $N^{\left({\mathrm{NO}}_{\mathrm{3}}\right)}$, and N^(Si) for P^(DIATOM) and $N^{\left({\mathrm{PO}}_{\mathrm{4}}\right)}$ and $N^{\left({\mathrm{NO}}_{\mathrm{3}}\right)}$ for the rest of the phytoplankton types. Light regulates gpp through an increasing, saturating function of the number of photons absorbed per Chl a unit, the photochemical efficiency (ϕ, mg C µE⁻¹) and the Chl a ($P_{\text{Chl}a}^{\left(i\right)})$-to-carbon ($P_{\text{C}}^{\left(i\right)})$ ratio (θ, mg Chl a  mg C⁻¹). Total absorbed irradiance is the integral of the Chl a-specific absorption spectrum of photosynthetic pigments [$a_{\text{PS}}^{\ast }\left(\mathit{\lambda }\right)$] times E₀(λ) (see Sect. 2.1.2.3 and Eq. 26) between 387.5 and 800 nm (µmol quanta mg Chl a⁻¹ d⁻¹). The photoacclimation term used in the calculation of $P_{\text{Chl}a}^{\left(i\right)}$ is resolved spectrally in the same way, and the exact formulation can be found in Álvarez et al. (2022).

DOC in the BFM is characterized by reactivity into labile ($R_{\text{C}}^{\left(\mathrm{1}\right)})$, semi-labile ($R_{\text{C}}^{\left(\mathrm{2}\right)})$, and semi-refractory ($R_{\text{C}}^{\left(\mathrm{3}\right)})$, and, accordingly, CDOM is represented by the same three reactivity levels ($X_{\text{C}}^{\left(\mathrm{1}\right)}$, $X_{\text{C}}^{\left(\mathrm{2}\right)}$, and $X_{\text{C}}^{\left(\mathrm{3}\right)})$. CDOM is calculated as a fraction of the DOC produced by phytoplankton that only produce $R_{\text{C}}^{\left(\mathrm{2}\right)}$, zooplankton that only produce $R_{\text{C}}^{\left(\mathrm{1}\right)}$, and bacteria that only produce $R_{\text{C}}^{\left(\mathrm{3}\right)}$ (green lines in Fig. 1, all in mg C m⁻³ d⁻¹).

With respect to the semi-labile DOC produced by phytoplankton ($R_{\text{C}}^{\left(\mathrm{2}\right)})$, the BFM considers two components: (i) a constant percentage of carbon fixed by phytoplankton is lost from cells each day, along with dissolved organic phosphorus and nitrogen (DOP and DON, respectively), and (ii) photosynthetic overflow under conditions of nutrient limitation leads to the loss of photosynthesized carbon, which cannot be assimilated into biomass due to the absence of other nutrients. Following Thornton's (2014) definitions, the first process is defined as the passive leakage of molecules through the cell membrane, and the second process of exudation is defined as the loss of excess carbon due to changes in environmental conditions such as nutrient availability. Thus, the total amount of phytoplankton dpp in the BFM is a dynamic term that changes in response to environmental conditions and is calculated as follows:

$$\begin{array}{}\text{(3)}& \left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{dpp}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}={\mathit{\beta }}_{\mathrm{P}}^{\left(i\right)}\times \left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{gpp}}{O^{\left(\mathrm{DIC}\right)}}}+\left(G_{\mathrm{P}}-G_{\mathrm{P}}^{\mathrm{balance}}\right),\end{array}$$

where the first addend in the right-hand side of Eq. (3) represents the passive leakage of molecules across the cell membrane and is considered to be a fraction (β_P) of the gpp, and the second addend represents the difference between the carbon assimilated into biomass (G_P) and the total photosynthesized carbon in the case of intracellular nutrient deficiency ($G_{\text{P}}^{\text{balance}})$. For each PFT, G_P(i) is calculated as gpp minus dpp and rsp, and $G_{\text{P}\left(i\right)}^{\text{balance}}$ is calculated as follows:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& G_{\mathrm{P}\left(i\right)}^{\mathrm{balance}}=min\left(\right.G_{\mathrm{P}\left(i\right)},{\displaystyle \frac{\mathrm{1}}{p_{\mathrm{P}}^{\min }}}\times \left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{P}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{N^{\left({\mathrm{PO}}_{\mathrm{4}}\right)}}},\\ & {\displaystyle \frac{\mathrm{1}}{n_{\mathrm{P}}^{\min }}}\times \sum _{k={\mathrm{NO}}_{\mathrm{3}},{\mathrm{NH}}_{\mathrm{4}}}\left.{\displaystyle \frac{\mathrm{d}{P}_{\text{N}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{N^{\left(k\right)}}}\left.\right),\end{array}\end{array}$$

where $P_{\text{P}}^{\text{min}}$ and $n_{\text{P}}^{\text{min}}$ are the minimum phosphorous and nitrogen quota, respectively, and the two nitrogen sources are $N^{\left({\mathrm{NO}}_{\mathrm{3}}\right)}$ and $N^{\left({\mathrm{NH}}_{\mathrm{4}}\right)}$. The exact formulations for the processes of phosphorous and nitrogen compounds uptake (upt) can be consulted in Vichi et al. (2020).

To investigate what fraction of these two $R_{\text{C}}^{\left(\mathrm{2}\right)}$ fluxes is chromophoric and whether this proportion varies with environmental conditions, we included light and nutrient dependencies on phytoplankton CDOM production. Assuming that modelled exudation is dominated by non-chromophoric carbohydrates, under decreasing conditions of nutrient availability where exudation accounts for a larger proportion of dpp, the simulated proportion of chromophoric material in total dpp decreases. On the other hand, with a modelled leakage that contains more chromophoric material when cells have a higher intracellular pigment concentration, the simulated chromophoric fraction of the total dpp increases with decreasing light availability. Thus, whereas the DOC that leaked constantly from phytoplankton cells is considered to have an amount that is chromophoric and proportional to cell pigmentation ($f_{\text{P}}^{X\mathrm{2}}\times \mathit{\theta }/{\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}})$, DOC exuded under intra-cellular nutrient shortage is considered to be non-absorbing, and therefore it is not partitioned into CDOM. The coloured fraction of the total DOC of phytoplankton origin ($f_{R\mathrm{2}}^{X\mathrm{2}})$ is the sum of all coloured material produced by PFTs with respect to total DOC production. The coloured fraction produced by each phytoplankton type i ($f_{\text{P}\left(i\right)}^{X\mathrm{2}})$ depends on the relative proportions of DOC that originated through leakage or exudation and on the intracellular concentration of Chl a as follows:

$$\begin{array}{}\text{(5)}& f_{\mathrm{P}\left(i\right)}^{X\mathrm{2}}=f_{\mathrm{P}\left(i\right)}^{\max X\mathrm{2}}\times {\displaystyle \frac{{\mathit{\theta }}^{\left(i\right)}}{{\mathit{\theta }}_{\mathrm{Chl}}^{\mathrm{0}\left(i\right)}}}\times {\displaystyle \frac{{\mathit{\beta }}_{\mathrm{P}}^{\left(i\right)}\times \left.\frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}\right|\frac{\mathrm{gpp}}{O^{\left(\mathrm{DIC}\right)}}}{\left.\frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}\right|\frac{\mathrm{dpp}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}}.\end{array}$$

The non-absorbing part of this flux (1-$f_{R\mathrm{2}}^{X\mathrm{2}})$ is directed to the state variable $R_{\text{C}}^{\left(\mathrm{2}\right)}$ and consumed by bacteria (upt):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& \mathrm{d}{R}_{\mathrm{C}}^{\left(\mathrm{2}\right)}/\mathrm{d}t=\sum _{i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}\left(\left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{dpp}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}\times \left(\mathrm{1}-f_{\mathrm{P}\left(i\right)}^{X\mathrm{2}}\right)\right)\\ & -\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}.\end{array}\end{array}$$

CDOM is directed to the state variable $X_{\text{C}}^{\left(\mathrm{2}\right)}$, consumed by bacteria (upt), and photobleached (deg):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{2}\right)}/\mathrm{d}t=\sum _{i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{4}}\left(\left.{\displaystyle \frac{\mathrm{d}{P}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{dpp}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}\times f_{\mathrm{P}\left(i\right)}^{X\mathrm{2}}\right)\\ & -\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{R_{\mathrm{C}}^{\left(\mathrm{2}\right)}}}-\left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{2}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{deg}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}.\end{array}\end{array}$$

The formulations for the processes of bacterial consumption (upt) can be consulted in Vichi et al. (2020).

Labile DOC ($R_{\text{C}}^{\left(\mathrm{1}\right)})$ is produced by the excretion of microzooplankton (Z^(MICRO)) and nano-heterotrophs (Z^(HETEN)) (exc) and therefore represents the sources of DOC associated with the zooplankton-mediated mortality of phytoplankton and bacteria, and it is consumed quickly by bacteria. The released fraction of carbon by mesozooplankton excretion (both Z^(OMNI) and Z^(CARNI)), on the other hand, is assumed to have no dissolved products and is therefore directed to the state variable $R_{\text{C}}^{\text{(DET)}}$.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& \mathrm{d}{R}_{\mathrm{C}}^{\left(\mathrm{1}\right)}/\mathrm{d}t=\sum _{j=\mathrm{HETEN},\mathrm{MICRO}}\left(\left.{\displaystyle \frac{\mathrm{d}{Z}_{\mathrm{C}}^{\left(j\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{exc}}{R_{\mathrm{C}}^{\left(\mathrm{1}\right)}}}\times \left(\mathrm{1}-f_{\mathrm{Z}\left(j\right)}^{X\mathrm{1}}\right)\right)\\ & -\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{R_{\mathrm{C}}^{\left(\mathrm{1}\right)}}}\end{array}\end{array}$$

Labile CDOM ($X_{\text{C}}^{\left(\mathrm{1}\right)})$ is explicitly resolved as a constant fraction of exc determined by the parameter $f_{Z\text{(MICRO)}}^{X\mathrm{1}}$ for microzooplankton and $f_{Z\text{(HETEN)}}^{X\mathrm{1}}$ for nano-heterotrophs:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& \mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{1}\right)}/\mathrm{d}t=\sum _{j=\mathrm{HETEN},\mathrm{MICRO}}\left(\left.{\displaystyle \frac{\mathrm{d}{Z}_{\mathrm{C}}^{\left(j\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{exc}}{R_{\mathrm{C}}^{\left(\mathrm{1}\right)}}}\times f_{\mathrm{Z}\left(j\right)}^{X\mathrm{1}}\right)\\ & -\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{upt}}{R_{\mathrm{C}}^{\left(\mathrm{1}\right)}}}-\left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{1}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{deg}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}.\end{array}\end{array}$$

The formulations for the processes of zooplankton excretion (exc) can be consulted in Vichi et al. (2020).

The production of DOC by phytoplankton and as a by-product of zooplankton activities is an important driver of secondary production by heterotrophic prokaryotes. Labile and semi-labile DOC, both non-absorbing ($R_{\text{C}}^{\left(\mathrm{1}\right)}$ and $R_{\text{C}}^{\left(\mathrm{2}\right)})$ and coloured ($X_{\text{C}}^{\left(\mathrm{1}\right)}$ and $X_{\text{C}}^{\left(\mathrm{2}\right)})$ are consumed (upt) by bacteria and excreted (exc) in the form of semi-refractory DOC ($R_{\text{C}}^{\left(\mathrm{3}\right)})$, resulting in bacterial DOC recycling that contributes to the sequestration of organic carbon:

$$\begin{array}{}\text{(10)}& \mathrm{d}{R}_{\mathrm{C}}^{\left(\mathrm{3}\right)}/\mathrm{d}t=\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{exc}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}\times \left(\mathrm{1}-f_B^{X\mathrm{3}}\right)-\left.{\displaystyle \frac{\mathrm{d}{R}_{\mathrm{C}}^{\left(\mathrm{3}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{rem}}{O^{\left(\mathrm{DIC}\right)}}}.\end{array}$$

Part of this recycled DOC is coloured, and therefore $X_{\text{C}}^{\left(\mathrm{3}\right)}$ is explicitly resolved as a constant fraction of $R_{\text{C}}^{\left(\mathrm{3}\right)}$ production determined by the parameter $f_B^{X\mathrm{3}}$:

$$\begin{array}{}\text{(11)}& \mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{3}\right)}/\mathrm{d}t=\left.{\displaystyle \frac{\mathrm{d}{B}_{\mathrm{C}}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{exc}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}\times f_B^{X\mathrm{3}}-\left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{3}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{deg}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}-\left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{3}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{rem}}{O^{\left(\mathrm{DIC}\right)}}}.\end{array}$$

Neither $R_{\text{C}}^{\left(\mathrm{3}\right)}$ nor $X_{\text{C}}^{\left(\mathrm{3}\right)}$ are explicitly consumed by bacteria, but both are remineralized (rem) at a very slow temperature-controlled ($Q_{{\mathrm{10}}_{\mathrm{R}}})$ rate (r_R, d⁻¹):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \left.{\displaystyle \frac{\mathrm{d}{R}_{\mathrm{C}}^{\left(\mathrm{3}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{rem}}{O^{\left(\mathrm{DIC}\right)}}}=Q_{{\mathrm{10}}_{\mathrm{R}}}^{\frac{T-\mathrm{10}}{\mathrm{10}}}\times r_{\mathrm{R}}\times R_{\mathrm{C}}^{\left(\mathrm{3}\right)};\\ & \left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(\mathrm{3}\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{rem}}{O^{\left(\mathrm{DIC}\right)}}}=Q_{{\mathrm{10}}_{\mathrm{R}}}^{\frac{T-\mathrm{10}}{\mathrm{10}}}\times r_{\mathrm{R}}\times X_{\mathrm{C}}^{\left(\mathrm{3}\right)}.\end{array}\end{array}$$

Unlike DOC, all three pools of CDOM are photobleached (deg) at a rate that reaches a maximum defined by the parameter b_X(i) (d⁻¹) when PAR is above the threshold defined by the parameter I_X(i) (µmol quanta m⁻² d⁻¹) and that decreases linearly at lower PAR values (Dutkiewicz et al., 2015):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \left.{\displaystyle \frac{\mathrm{d}{X}_{\mathrm{C}}^{\left(i\right)}}{\mathrm{d}t}}\right|{\displaystyle \frac{\mathrm{deg}}{R_{\mathrm{C}}^{\left(\mathrm{3}\right)}}}=b_{X\left(i\right)}\times min\left(\mathrm{PAR}/I_{X\left(i\right)},\mathrm{1}\right)\times X_{\mathrm{C}}^{\left(i\right)};\\ & i=\mathrm{1},\mathrm{2},\mathrm{3}.\end{array}\end{array}$$

Figure 1 Schematic representation of the hypothesized interactions regulating the concentrations of CDOM and DOC in the surface layers of the open waters of the NW Mediterranean Sea. The boxes show the BFM state variables, and the yellow triangles show the three state variables that represent three reactivities of CDOM. Green arrows indicate fluxes of DOC and CDOM production (excretion$R_{\text{C}}^{\mathrm{1}}$, excretion$R_{\text{C}}^{\mathrm{3}}$, exudation and leakage; mg C m⁻³ d⁻¹), orange arrows indicate fluxes of DOC and/or CDOM degradation through direct bacterial consumption (uptake$R_{\text{C}}^{\mathrm{1}}$, uptake$X_{\text{C}}^{\mathrm{1}}$, uptake$R_{\text{C}}^{\mathrm{2}}$, uptake$X_{\text{C}}^{\mathrm{2}}$; mg C m⁻³ d⁻¹) or indirect remineralization (rem$R_{\text{C}}^{\mathrm{3}}$, rem$X_{\text{C}}^{\mathrm{3}}$; mg C m⁻³ d⁻¹), and black arrows indicate photobleaching processes (deg$X_{\text{C}}^{\mathrm{1}}$, deg$X_{\text{C}}^{\mathrm{2}}$, deg$X_{\text{C}}^{\mathrm{3}}$; mg C m⁻³ d⁻¹). The names in italics ($f_{\text{Z}}^{X\mathrm{1}}$, $f_{R\mathrm{2}}^{X\mathrm{2}}$, and $f_B^{X\mathrm{3}}$; dimensionless) represent the fractions that divide the flux between CDOM and non-coloured DOC. The subscript i appended to each living component and detritus indicates the elemental constituents among carbon (C), nitrogen (N), and phosphorus (P), with Chl a only being used for phytoplankton and silica (Si) only being used for diatoms.

Bio-optical properties of CDOM, detritus, and phytoplankton

The optical constituents of the spectrally resolved version of BFM are the four PFTs, detrital particles, and the three forms of CDOM corresponding to their resistance to bacterial metabolic activity. CDOM absorbs light and has a negligible contribution to scattering; therefore, CDOM absorption at 450 nm, a_CDOM(450) (m⁻¹), is first calculated as the product of CDOM biomass ($X_{\text{C}}^{\left(i\right)}$, mg C m⁻³) and the mass-specific absorption coefficients at 450 nm ($a_{X\left(i\right)}^{\mathrm{450}}$, m² mg C⁻¹). CDOM absorption as a function of wavelength (λ), a_CDOM(λ) is modelled with an exponential function, decreasing with increasing λ:

$$\begin{array}{}\text{(14)}& a_{\text{CDOM}}\left(\mathit{\lambda }\right)=\sum _{i=\mathrm{1},\mathrm{2},\mathrm{3}}{a}_{X\left(i\right)}^{\mathrm{450}}\times \exp \left[-S_{X\left(i\right)}^{\mathrm{350}-\mathrm{500}}\times \left(\mathit{\lambda }-\mathrm{450}\right)\right]\times X_{\text{C}}^{\left(i\right)},\end{array}$$

where $S_{X\left(i\right)}^{\mathrm{350}-\mathrm{500}}$ is the spectral slope of the CDOM absorption coefficients between 350 and 500 nm.

Detritus absorbs and scatters light. Detritus absorption at 440 nm (a_NAP(440), m⁻¹) is calculated as the product of the carbon component of detritus ($R_{\text{C}}^{\text{(DET)}}$, mg C m⁻³) and their mass-specific absorption coefficients at 440 nm ($a_{\text{R}}^{\mathrm{440}}$, m² mg C⁻¹), and then the absorption spectrum of detritus is modelled as an exponential function which decreases with increasing λ:

$$\begin{array}{}\text{(15)}& a_{\mathrm{NAP}}\left(\mathit{\lambda }\right)=a_{\mathrm{R}}^{\mathrm{440}}\times \exp \left[-S_{\mathrm{R}}^{\mathrm{350}-\mathrm{500}}\times \left(\mathit{\lambda }-\mathrm{440}\right)\right]\times R_{\mathrm{C}}^{\left(\mathrm{DET}\right)},\end{array}$$

where $S_{\text{R}}^{\mathrm{350}-\mathrm{500}}$ is the spectral slope of the detritus absorption coefficients, which we have set to 0.013 nm⁻¹ (Table 1). Detrital scattering at 550 nm (b_NAP(550), m⁻¹) is calculated as the product of $R_{\text{C}}^{\text{(DET)}}$ and the mass-specific scattering coefficient at 550 nm ($b_{\text{R}}^{\mathrm{550}}$, m² m C⁻¹), and then the scattering spectrum of detritus, b_NAP(λ), is computed as an exponential function of λ:

$$\begin{array}{}\text{(16)}& b_{\mathrm{NAP}}\left(\mathit{\lambda }\right)=b_{\text{R}}^{\mathrm{550}}\times (\mathrm{550}/\mathit{\lambda }{)}^{e_{\text{R}}}\times R_{\mathrm{C}}^{\left(\text{DET}\right)}\end{array}$$

where e_R is an exponent set to 0.5 (Table 1).

To describe the optical properties of phytoplankton, we used Chl a-specific absorption spectra ($a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right))$ of different phytoplankton taxa growing in culture under different light, nutrient supply, and temperature conditions. We considered 177 spectra, digitized from the literature and provided as supplementary material in Álvarez et al. (2022). Each $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ was decomposed into the contribution of Chl a; photosynthetic carotenoids; photoprotective carotenoids (PPCs); and a fourth pigment group that was made up of phycoerythrin for Synechococcus species, Chl b for green algae taxa and Prochlorococcus species, and Chl c for the other taxa. To obtain the relative pigment concentrations for each pigment group we fitted a multiple-linear-regression model to the weight-specific absorption spectra for each pigment group obtained from Bidigare et al. (1990) in order to obtain the lowest sum of residuals between the predicted and the measured $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ (Hickman et al., 2010). The Chl a-specific absorption spectra for photosynthetic pigments only ($a_{\text{PS}}^{\ast }\left(\mathit{\lambda }\right))$ were reconstructed using the obtained relative concentrations excluding PPCs (Babin et al., 1996). From this collection of spectra, $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ and $a_{\text{PS}}^{\ast }\left(\mathit{\lambda }\right)$ were averaged for picophytoplankton, nanophytoplankton, diatoms, and dinophyta taxa. The resulting four $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ were used to calculate the total absorption by phytoplankton (a_PH(λ), m⁻¹) as the sum of the products of PFT Chl a ($P_{\text{Chl}a}^{\left(i\right)}$, mg Chl a m⁻³) and the Chl a-specific absorption coefficients (m² mg Chl a⁻¹):

$$\begin{array}{}\text{(17)}& a_{\mathrm{PH}}\left(\mathit{\lambda }\right)=\sum _{i=\mathrm{DIATOM},\mathrm{NANO},\mathrm{PICO},\mathrm{DINO}}{a}_{\mathrm{PH}}^{\ast \left(i\right)}\left(\mathit{\lambda }\right)\times P_{\mathrm{Chl}a}^{\left(i\right)}.\end{array}$$

Scattering coefficients of phytoplankton were assumed to be functions of phytoplankton biomass ($P_{\text{C}}^{\left(i\right)}$, mg C m⁻³) and the C-specific spectra taken from Dutkiewicz et al. (2015). Scattering coefficients by all phytoplankton (b_PH(λ), m⁻¹) were calculated as the sum of the product of phytoplankton biomass and the C-specific scattering spectrum for each of the four PFTs (m² mg C⁻¹):

$$\begin{array}{}\text{(18)}& b_{\mathrm{PH}}\left(\mathit{\lambda }\right)=\sum _{i=\mathrm{DIATOM},\mathrm{NANO},\mathrm{PICO},\mathrm{DINO}}{b}_{\mathrm{PH}}^{\ast \left(i\right)}\left(\mathit{\lambda }\right)\times P_{\mathrm{C}}^{\left(i\right)}.\end{array}$$

The phytoplankton absorption cross-sections for all pigments (${\stackrel{\mathrm{‾}}{a}}_{\text{PH}}^{\ast }$, m² mg Chl a⁻¹), photosynthetic pigments (${\stackrel{\mathrm{‾}}{a}}_{\text{PS}}^{\ast }$, m² mg Chl a⁻¹), and the scattering cross-section (${\stackrel{\mathrm{‾}}{b}}_{\text{PH}}^{\ast }$, m² mg Chl a⁻¹) were computed as the average of $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$, $a_{\text{PS}}^{\ast }\left(\mathit{\lambda }\right)$, and $b_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ between 387.5 and 800 nm, respectively. These aggregated values were used to perturb the magnitude of their respective spectra without modifying the spectral shape. The computation of ${\stackrel{\mathrm{‾}}{a}}_{\text{PS}}^{\ast }$ and ${\stackrel{\mathrm{‾}}{b}}_{\text{PH}}^{\ast }$ is equivalent to Eq. (19), substituting $a_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$ with $a_{\text{PS}}^{\ast }\left(\mathit{\lambda }\right)$ and $b_{\text{PH}}^{\ast }\left(\mathit{\lambda }\right)$, respectively:

$$\begin{array}{}\text{(19)}& {\stackrel{\mathrm{‾}}{a}}_{\mathrm{PH}}^{\ast }=\underset{\mathrm{387.5}}{\overset{\mathrm{800}}{\int }}{a}_{\mathrm{PH}}^{\ast }\left(\mathit{\lambda }\right)d\mathit{\lambda }/\left(\mathrm{800}-\mathrm{387.5}\right).\end{array}$$

Total absorption, scattering, and backscattering (a(λ), b(λ), and bb(λ), respectively) depend on the additive contribution of seawater and the BGC constituents along the water column. Total a(λ) (m⁻¹) is calculated from the absorption by water (a_W(λ), m⁻¹, Pope and Fry (1997)), CDOM, detritus, and phytoplankton:

$$\begin{array}{}\text{(20)}& a\left(\mathit{\lambda }\right)=a_{\mathrm{W}}\left(\mathit{\lambda }\right)+a_{\mathrm{CDOM}}\left(\mathit{\lambda }\right)+a_{\mathrm{NAP}}\left(\mathit{\lambda }\right)+a_{\mathrm{PH}}\left(\mathit{\lambda }\right).\end{array}$$

Total b(λ) (m⁻¹) is calculated as the sum of the scattering of water (b_W(λ), m⁻¹, Smith and Baker (1981)), detritus, and phytoplankton:

$$\begin{array}{}\text{(21)}& b\left(\mathit{\lambda }\right)=b_{\mathrm{W}}\left(\mathit{\lambda }\right)+b_{\mathrm{NAP}}\left(\mathit{\lambda }\right)+b_{\mathrm{PH}}\left(\mathit{\lambda }\right).\end{array}$$

Total bb(λ) (m⁻¹) is computed as the sum of each constituent's scattering multiplied by constant and λ-independent backscattering-to-scattering ratios (bbr_W, bbr_R, and $bb{r}_{\text{PH}}^{\left(i\right)})$ that were 0.5 for water (Morel, 1974), 0.005 for detritus (Gallegos et al., 2011), 0.002 for P^(DIATOM) (Dutkiewicz et al., 2015), 0.0071 for P^(NANO) (Gregg and Rousseaux, 2017), 0.0039 for P^(PICO) (Gregg and Rousseaux, 2017; Dutkiewicz et al., 2015), and 0.003 for P^(DINO) (Dutkiewicz et al., 2015) (Table 1):

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& bb\left(\mathit{\lambda }\right)=b_{\mathrm{W}}\left(\mathit{\lambda }\right)\times bb{r}_{\mathrm{W}}+b_{\mathrm{NAP}}\left(\mathit{\lambda }\right)\times bb{r}_{\mathrm{R}}\\ & +\sum _{i=\mathrm{DIATOM},\mathrm{NANO},\mathrm{PICO},\mathrm{DINO}}{b}_{\mathrm{PH}}^{\ast \left(i\right)}\left(\mathit{\lambda }\right)\times {\mathrm{P}}_{\mathrm{C}}^{\left(i\right)}\times bb{r}_{\mathrm{PH}}^{\left(i\right)}.\end{array}\end{array}$$

Spectrally resolved light transmission and computation of E₀(λ,z)

Incoming spectral irradiance at the top of the ocean was obtained from the OASIM model interfaced with the European Centre for Medium-Range Weather Forecast (ECMWF) atmospheric model (Lazzari et al., 2021a). Input data contained two downward streams just below the surface of the ocean: direct [E_d(λ, 0⁻)] and diffuse [$E_{\text{s}}(\mathit{\lambda },{\mathrm{0}}^{-})$] (both in W m⁻² and averaged at 33 wavebands). Irradiance along the depth of the water column (z) was described with three state variables (all in W m⁻²): downwelling direct [E_d(λ,z)] and diffuse [E_s(λ,z)] components and an upward diffuse component [E_u(λ,z)]. The vertically resolved three-stream propagation was resolved with the following system of differential equations (Aas, 1987; Ackleson et al., 1994; Gregg, 2002):

$$\begin{array}{}\text{(23)}& {\displaystyle }\mathrm{d}{E}_{\mathrm{d}}\left(\mathit{\lambda },z\right)/\mathrm{d}Z=-{\displaystyle \frac{a\left(\mathit{\lambda },z\right)+b\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{d}}}}\times E_{\mathrm{d}}\left(\mathit{\lambda },z\right),\text{(24)}& {\displaystyle }\begin{array}{rl}& \mathrm{d}{E}_{\mathrm{s}}\left(\mathit{\lambda },z\right)/\mathrm{d}Z=-{\displaystyle \frac{a\left(\mathit{\lambda },z\right)+r_{\mathrm{s}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{s}}}}\\ & \times E_{\mathrm{s}}\left(\mathit{\lambda },z\right)+{\displaystyle \frac{r_{\mathrm{u}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{u}}}}\times E_{\mathrm{u}}\left(\mathit{\lambda },z\right)\\ & +{\displaystyle \frac{b\left(\mathit{\lambda },z\right)-r_{\mathrm{d}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{d}}}}\times E_{\mathrm{d}}\left(\mathit{\lambda },z\right)\end{array},\text{(25)}& {\displaystyle }\begin{array}{rl}& -\mathrm{d}{E}_{\mathrm{u}}\left(\mathit{\lambda },z\right)/\mathrm{d}Z=-{\displaystyle \frac{a\left(\mathit{\lambda },z\right)+r_{\mathrm{u}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{u}}}}\\ & \times E_{\mathrm{u}}\left(\mathit{\lambda },z\right)+{\displaystyle \frac{r_{\mathrm{s}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{s}}}}\times E_{\mathrm{s}}\left(\mathit{\lambda },z\right)\\ & +{\displaystyle \frac{r_{\mathrm{d}}\times bb\left(\mathit{\lambda },z\right)}{{\mathit{\upsilon }}_{\mathrm{d}}}}\times E_{\mathrm{d}}\left(\mathit{\lambda },z\right),\end{array}\end{array}$$

where a(λ,z) is the total absorption (Eq. 20), b(λ,z) is the total scattering (Eq. 21), and bb(λ,z) is the total backward scattering (Eq. 22) at depth z (all in m⁻¹). r_s, r_u and r_d are the effective scattering coefficients, and υ_d, υ_s, and υ_u are the average cosines of the three streams, which are constant for diffuse radiance but vary with solar zenith angle for direct radiance. For details on the derivation of the solution, see Terzić et al. (2021).

At the centre of the depth layer, E_d(λ,z), E_s(λ,z), and E_u(λ,z) were averaged between the top and the bottom of the layer, and total scalar irradiance, E₀(λ,z), was computed by scaling the three streams by their inverse average direction cosines. E₀(λ,z) was converted from irradiance values (in units of W m⁻²) to photon flux (given in µmol quanta m⁻² d⁻¹) by multiplying by λ in metres, dividing by the product of the Avogadro's constant (N_A), Planck's constant (h), and speed of light (c) and converting mol quanta to µMol quanta and seconds to days:

$$\begin{array}{}\text{(26)}& \begin{array}{rl}& E_{\mathrm{0}}\left(\mathit{\lambda },z\right)=\left[\right.E_{\mathrm{d}}\left(\mathit{\lambda },z\right)/\stackrel{\mathrm{‾}}{{\mathit{\upsilon }}_{\mathrm{d}}}+E_{\mathrm{s}}\left(\mathit{\lambda },z\right)/\stackrel{\mathrm{‾}}{{\mathit{\upsilon }}_{\mathrm{s}}}\\ & +E_{\mathrm{u}}\left(\mathit{\lambda },z\right)/\stackrel{\mathrm{‾}}{{\mathit{\upsilon }}_{\mathrm{u}}}\left.\right]\times {\displaystyle \frac{\mathit{\lambda }}{h\cdot c\cdot N_{\mathrm{A}}}}\\ & \times {\displaystyle \frac{\mathrm{1}\mathrm{m}}{{\mathrm{10}}^{\mathrm{9}}\mathrm{nm}}}\times {\displaystyle \frac{{\mathrm{10}}^{\mathrm{6}}\mathrm{\mu }\mathrm{mol}\mathrm{quanta}}{\mathrm{1}\mathrm{mol}\mathrm{quanta}}}\times {\displaystyle \frac{\mathrm{86}\mathrm{400}\mathrm{s}}{\mathrm{1}\mathrm{d}}}.\end{array}\end{array}$$

E₀(λ,z) was the light available for phytoplankton growth (Eq. 2) at depth z, and its integral value from 387.5 and 800 nm constituted the photosynthetically available radiation (PAR) that was used as the light input for CDOM photobleaching (Sect. 2.1.2.1, Eq. 13).

2.2 Observations

The Mediterranean Sea is a semi-enclosed marginal sea. The Ligurian Sea, in the NW Mediterranean Sea, is characterized by oligotrophic to moderately eutrophic waters. The BOUSSOLE and the DYFAMED (Dynamique des Flux Atmosphériques en Méditerranée) monitoring sites, both in the Ligurian Sea, provided bio-optical and physical–chemical measurements, respectively (Sect. 2.2.1). In addition, we considered ocean-colour data collected by satellites (Sect. 2.2.2) (Fig. 2).

2.2.1 In situ data at BOUSSOLE and DYFAMED sites

The BOUSSOLE project includes a permanent optical mooring deployed at 43^∘22^′ N, 7^∘54^′ E, where the depth is 2440 m, collecting optical data at a high temporal resolution, and oceanographic cruises that visit the site monthly to collect optical and BGC variables (Golbol et al., 2000). Only data from monthly cruises were used in the present work. 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 CDOM absorption measurements are described in Organelli et al. (2014) and references therein. We only remind here that the spectral slope of a_CDOM(λ) was determined between 350 and 500 nm ($S_{\text{CDOM}}^{\mathrm{350}-\mathrm{500}})$, although shorter UV λ intervals have also been used to report spectral slopes in the Mediterranean Sea (275–295 nm, Catalá et al., 2018). Particulate light absorption spectra, a_P(λ), were measured using the quantitative filter pad technique (QFT) (Mitchell, 1990; Mitchell et al., 2000), and the protocol is described in Antoine et al. (2006). a_P(λ) was decomposed into phytoplankton and detrital absorption coefficients, a_PH(λ) and a_NAP(λ), respectively, using the numerical decomposition technique of Bricaud and Stramski (1990). The HPLC procedure used to measure phytoplankton pigment concentrations is fully described in Ras et al. (2008). The fractions of picophytoplankton, nanophytoplankton, and microphytoplankton were calculated following the diagnostic pigments method (Uitz et al., 2006).

The Dynamique des Flux Atmosphériques en Méditerranée (DYFAMED) time series station (Marty, 2002) is located at 43.25^∘ N and 7.54^∘ E, where data have been collected since 1991. Temperature (T) and salinity (S) data were collected within the 0–2400 m water column with the same CTD/rosette configuration used at BOUSSOLE with nutrients sampled at discrete depths. T and S data from 1994 to 2014 were downloaded from the OceanSITES Global Data Assembly Centre, and nitrate, nitrite, and phosphate concentrations were downloaded from Coppola et al. (2021).

2.2.2 Data from satellites

Satellite-retrieved IOPs and Chl a were obtained from the OCEANCOLOUR_MED_BGC_L3_MY_009_143 product (EU Copernicus Marine Service Information, 2022). The product provides multi-sensor R_RS(λ) spectra (SeaWiFS, MODIS-AQUA, NOAA20-VIIRS, NPP-VIIRS, and Sentinel3A-OLCI) merged using the state-of-the-art multi-sensor merging algorithms and shifted to the standard SeaWiFS wavelengths (412, 443, 490, 555, and 670 nm). IOPs (a_PH(443) and a_DG(443)) are derived from R_RS(λ) via the Quasi-Analytical Algorithm version 6 (QAAv6) model (Lee et al., 2002). Chl a concentration is determined using the Mediterranean regional algorithm, an updated version of MedOC4 (Volpe et al., 2019) and AD4 (Berthon and Zibordi, 2004), merged according to the water type present (Mélin and Vantrepotte, 2015; Volpe et al., 2019). For all variables, to obtain the time series from January 2011 to December 2014 at the study site, we extracted a 20×20 km box around the BOUSSOLE site (from 43.27 to 43.47^∘ N in latitude and from 7.77 to 8.03^∘ E in longitude) and averaged it when available observations covered at least 20 % of the area.

Figure 2 Map showing the location of the BOUSSOLE mooring site (red cross) and the DYFAMED site (black cross) in the Ligurian Sea (NW Mediterranean Sea). The box surrounding the sites indicates the 20×20 km area considered for averaging data from the daily 1 km Copernicus Marine Service product to obtain the time series of satellite data at the study site.

2.3 Study region and setup of the 1D model for the BOUSSOLE site

The configuration of the one-dimensional spectrally resolved version of the BFM used to simulate the BOUSSOLE site has the same parameter values as the 3D configuration of the BFM for the Mediterranean Sea, except for those reported in Table 1. Some updates comprised parameters related to the growth of nanophytoplankton (P^(NANO)), in particular those related to the uptake of phosphorous (a₁, $P_{\text{P}}^{\text{min}}$, and $P_{\text{P}}^{\text{opt}})$. This update was set to reproducing the high contribution of nanophytoplankton to the total Chl a content at the BOUSSOLE site (Antoine et al., 2006) (see Fig. 4). We also updated the parameters for the production by phytoplankton and the consumption by bacteria of semi-labile DOC (β_P and $v_B^{R\mathrm{2}})$ to simulate accurately the seasonality of DOC concentrations in the superficial waters of the site (Avril, 2002) (see Fig. 5). The parameter values related to the partition between dissolved and particulate excretion in Z^(MICRO) and Z^(HETEN) and the mortality of Z^(CARNI) and Z^(OMNI) are equal to the up-to-date values in Álvarez et al. (2022) together with all the optical parameters. All the other parameter values are equal to those shown in Lazzari et al. (2012).

The 1D configuration of GOTM-(FABM)-BFM for the BOUSSOLE site consists of (i) a water column of 2448 m subdivided into 196 vertical levels with a higher resolution at the surface and bottom, which were generated using the iGOTM web tool (https://igotm.bolding-bruggeman.com/, last access: 12 October 2021); (ii) atmospheric forcing that included hourly 10 m winds, 2 m air temperature and humidity, sea level pressure, precipitation, and cloud cover from ECMWF ERA5; atmospheric pCO₂ that was set to a constant value of 400 ppm; and spectral plane downwelling irradiance at 15 min frequency from the OASIM model that has been validated for the BOUSSOLE site (Lazzari et al., 2021a); and (iii) initial conditions for T, S, and BGC variables.

Observed vertical profiles of T and S at the DYFAMED site from 2009 to 2014 were edited to match the GOTM profile format and were used as initial conditions and constrained by a monthly nudging in order to reproduce the intensity and timing of the mixing as closely as possible to observations. Monthly vertical profiles of BGC variables (nitrate, phosphate, silicate, O₂, DIC, and alkalinity) were obtained from the Copernicus reanalysis (Cossarini et al., 2021), linearly interpolated onto the model vertical grid, and used for initialization and constrained by a yearly nudging. All other variables were set to constant values throughout the water column. The initial values of $R_{\text{C}}^{\left(\mathrm{1}\right)}$, $R_{\text{C}}^{\left(\mathrm{2}\right)}$, and $R_{\text{C}}^{\left(\mathrm{3}\right)}$ were 80, 400, and 480 mg C m⁻³, respectively, consistent with DOC concentrations at depth in the western Mediterranean Sea (Catalá et al., 2018; Santinelli, 2015; Galletti et al., 2019). The initial values of $X_{\text{C}}^{\left(\mathrm{1}\right)}$, $X_{\text{C}}^{\left(\mathrm{2}\right)}$, and $X_{\text{C}}^{\left(\mathrm{3}\right)}$ were 1, 3, and 1.25 mg C m⁻³, respectively. Given that below 200 m depth the BFM simulates negligible concentrations of $X_{\text{C}}^{\left(\mathrm{1}\right)}$ and $X_{\text{C}}^{\left(\mathrm{2}\right)}$, the initial concentration of $X_{\text{C}}^{\left(\mathrm{3}\right)}$ was chosen to match a_CDOM(450) values observed at the site below 100 m (Organelli et al., 2014) and a_CDOM(250) and a_CDOM(325) values observed in the western Mediterranean Sea below 200 m (Catalá et al., 2018). All PFTs were initialized to 8 mg C m⁻³ and 0.16 mg Chl a m⁻³. Bacteria and zooplankton were initialized to 1 mg C m⁻³.

A 6-year hindcast (January 2009 to December 2014) was conducted, covering the period for which all the necessary forcing data were available. The first 2 years were considered as model spin-up, and the next 4 years were considered for the optimization. The source for the setup can be found in the “Code availability” section.

Table 1 Values of the BFM optical and BGC parameters that were not optimized. The parameters with the source “this study” were manually chosen in order to obtain a good approximation of model results to findings by Avril (2002) and Antoine et al. (2006) with regard to the concentration of DOC at DYFAMED and PFT Chl a at BOUSSOLE, respectively.

2.4 Optimization strategy and simulations for hypothesis testing

The final step in completing the GOTM-(FABM)-BFM configuration at the BOUSSOLE site was to integrate the model and observations to drive the simulated output as closely as possible to the observations. Being computationally light, the one-dimensional BFM configuration has the advantage of allowing many simulations in a reasonable time, which provided us with a powerful tool for optimizing a relatively large number of parameters (Sect. 2.4.1) through a genetic algorithm that minimizes the misfit between the model and the observations (Sect. 2.4.2) and for investigating processes based on the proposed new parameterizations (Sect. 2.4.3).

2.4.1 Parameters optimized

Observations used for optimization were all collected at the BOUSSOLE site at a monthly temporal resolution and roughly at the same number of discrete depths and included pico-, nano-, and micro-phytoplankton Chl a (pico Chl a, nano Chl a, and micro Chl a, respectively); a_PH(450); a_NAP(450); a_CDOM(450); and $S_{\text{CDOM}}^{\mathrm{350}-\mathrm{500}}$. The potential number of independent parameters included in the optimization problem is limited by the observations available and by the fact that optimized parameters must appear in some process independently. A total of 25 optical and BGC parameters were optimized (Table A1 in Appendix A), and all of them were given a positive range of values for initialization and evolution without conditions or trade-offs for the consistency between pairs of parameters. The rest of the parameters not included in the optimization kept constant values (Table 1). The only observations related to CDOM were a_CDOM(λ), which depend on both the mass concentration and the mass-specific absorption coefficients across the electromagnetic spectrum. We assumed the absorption coefficients of the three reactivities to be constant at 450 nm ($a_{X\mathrm{1}}^{\mathrm{450}}$, $a_{X\mathrm{2}}^{\mathrm{450}}$, and $a_{X\mathrm{3}}^{\mathrm{450}})$, and we optimized parameters related to the production of CDOM ($f_{Z^{\text{(MICRO)}}}^{X\mathrm{1}}$, $f_{Z^{\text{(HETEN)}}}^{X\mathrm{1}}$, $f_{P^{\text{(DIATOM)}}}^{X\mathrm{2}}$, $f_{P^{\text{(NANO)}}}^{X\mathrm{2}}$, $f_{P^{\text{(PICO)}}}^{X\mathrm{2}}$, $f_{P^{\text{(DINO)}}}^{X\mathrm{2}}$, and $f_B^{X\mathrm{3}})$ and those related to the photobleaching of $X_{\text{C}}^{\left(\mathrm{3}\right)}$ (b_X(3) and I_X(3)). The spectral slope of the CDOM pool has been associated with aromaticity and the average molecular weight of the CDOM compounds (Blough and Green, 1995) but appears to be less linked to CDOM concentration. Therefore, the spectral slopes between 350 and 500 nm for the three reactivities ($S_{X\left(\mathrm{1}\right)}^{\mathrm{350}-\mathrm{500}}$, $S_{X\left(\mathrm{2}\right)}^{\mathrm{350}-\mathrm{500}}$, and $S_{X\left(\mathrm{3}\right)}^{\mathrm{350}-\mathrm{500}})$ were optimized to ensure the proper simulation of a_CDOM(λ) in wavebands other than 450 nm. For phytoplankton, independent observations regarding mass concentrations (Chl a) and optical properties (a_PH(λ)) were available; therefore, both types of parameters were optimized. We optimized the absorption cross-sections (${\stackrel{\mathrm{‾}}{a}}_{\text{PH}}^{\ast })$ and the maximum Chl a : C quota (${\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}})$ as they appear independently in light attenuation (Eq. 17) and in photoacclimation (Eq. 19 in Álvarez et al., 2022) and in CDOM production (Eq. 5), respectively. The absorption cross-sections of photosynthetic pigments (${\stackrel{\mathrm{‾}}{a}}_{\text{PS}}^{\ast })$ and the photochemical efficiencies (${\mathit{\varphi }}_{\text{C}}^{\mathrm{0}})$ of the PFTs appear only in photosynthesis (Eq. 2). We opted to optimize only one of these two interdependent parameters and therefore kept constant ${\stackrel{\mathrm{‾}}{a}}_{\text{PS}}^{\ast }$ and optimized ${\mathit{\varphi }}_{\text{C}}^{\mathrm{0}}$. The choice responds to the fact that ${\stackrel{\mathrm{‾}}{a}}_{\text{PS}}^{\ast }$ can be inferred from the observed a_PH(λ) of multiple phytoplankton taxa, and ${\mathit{\varphi }}_{\text{C}}^{\mathrm{0}}$, instead, is not well documented in the literature, being more difficult to measure. For detritus, we optimized the reference absorption coefficients at 440 nm ($a_{\text{R}}^{\mathrm{440}})$ and kept constant all the BGC parameters that alter the mass concentrations. Given that a_NAP(λ) observations were available, optical parameters related to the spectral shape of absorption by detritus ($S_{\text{R}}^{\mathrm{350}-\mathrm{500}})$ could have been included in the optimization, but we decided to maintain this parameter as constant given the small contribution of a_NAP(λ) to the total non-water absorption as compared to CDOM and phytoplankton. All parameters related to scattering and backscattering of both phytoplankton (${\stackrel{\mathrm{‾}}{b}}_{\text{PH}}^{\ast }$ and bbr_PH) and detritus ($b_{\text{R}}^{\mathrm{550}}$, e_R, and bbr_R) were not included in the optimization.

2.4.2 Parameter optimization method

To perform the estimation of model parameter values using observational data, we used the Parallel Sensitivity Analysis and Calibration utility (ParSAC, Bruggeman and Bolding, 2020). ParSAC applies the Differential Evolution (DE) algorithm (Storn and Price, 1997), which is inspired by the rules of natural selection. Different individuals have different fitness values in terms of how well they simulate observations. A higher fitness value indicates a combination of parameters that reproduce better observations and therefore must persist in the next generation. ParSAC formulates such a fitness (i.e. the probability that the candidate parameter values are the true parameter set representing reality) as a multi-objective function calculated as the log-transformed likelihood between the outcome of the model and the observations provided, assuming the observed values are log-normally distributed, with the median being equal to the model prediction. For any variable observed from a total number of M observed variables, there will be N_m number of pairs consisting of observations O_n and corresponding simulation-based estimates P_n (model results at the temporal and spatial observation locations). The sum of squares of the residuals for each variable observed m (SSQ_m) is computed as follows:

$$\begin{array}{}\text{(27)}& {\text{SSQ}}_m=\sum _{n=\mathrm{1}}^{N_m}{\left[P_n-O_n\right]}^{\mathrm{2}}.\end{array}$$

For M variables observed, all differences between the model and observations are combined as follows:

$$\begin{array}{}\text{(28)}& \text{lnLikelihood}=\sum _{m=\mathrm{1}}^M-N_m\times \ln \left({\text{SSQ}}_m\right).\end{array}$$

Each individual in a given population of size S=288 is a 25-dimensional target vector that represents a candidate solution to the problem. The parameter values in the target vectors in the first generation are sampled uniformly between minimum and maximum values, which are listed in Table A1 (Appendix A). The DE algorithm creates new generations of individuals by applying cycles of mutation, crossover, and selection operating on the target vectors. Mutation means that, once all simulations in a generation are completed, the DE algorithm creates a new set of S mutant vectors. In this implementation, each mutant vector is created by randomly selecting three targets (x_r1, x_r2, and x_r3) and applying a mutation operator that consists of $x_{\mathrm{r}\mathrm{1}}+\mathrm{0.5}\cdot (x_{\mathrm{r}\mathrm{2}}-x_{\mathrm{r}\mathrm{3}})$. Crossover means to increase the diversity of the parameter vectors; crossover is done to generate the final trial set of individuals. Here, crossover is performed between each target vector and its corresponding mutant vector, retaining each element of the target vector in the population with a probability of 0.9, otherwise introducing the corresponding element from the mutant vector. Selection, which is the final phase, is where the lnLikelihood of the trial vector is compared to the lnLikelihood of the corresponding target vector. The trial vector replaces the target if the lnLikelihood value obtained from the trial is higher than or equal to the lnLikelihood obtained from the target vector. After the selection procedure, the size of the population is again S. The cycle of mutation, crossover, and selection is carried out iteratively, and eventually, the DE algorithm provides an estimate of the optimal parameter values that minimizes the misfit between the model output and the observations (i.e. maximizes the sum of likelihoods in the population). Genetic algorithms are particularly robust in overcoming local minima (Storn and Price, 1997). As a parallel search technique that runs several vectors simultaneously, superior parameter configurations can help other vectors escape local minima and therefore avoid misconvergence. After passing through each generation, ParSAC checks whether the end condition is met, which, in this work, was a maximum number of 180 generations, after which the range of lnLikelihoods in the population was less than 0.005. The optimal solution for the 25 parameter values was extracted as the mean of the parameter distribution in generation 180. To determine the success achieved through optimization, we calculated the correlation coefficient (R) and the root mean squared error (RMSE) for each observation compared to the model output for each generation.

2.4.3 Simulations and hypothesis-testing experiments

The model configuration described in Sect. 2.3 and the optimized parameter values obtained as described in Sect. 2.4.2 represent the optimized model configuration. We used a series of hypothesis-testing experiments to investigate the mechanisms that determine the vertical and temporal patterns of CDOM distribution at the BOUSSOLE site. Those experiments consisted of three more simulations in which we changed some of the assumptions of the optimized configuration. All simulations are summarized in Table 2, all were run for 6 years with 2 years of spin-up, and the results were analysed for the last 4 years (January 2011 to December 2014).

In the first experiment, EXP-Biology, we investigated the role of nutrient and light limitations in controlling the biogenic in situ production of biodegradable CDOM and therefore in generating the observed a_CDOM(λ) values. This experiment consisted of comparing the results of the optimized configuration with two additional simulations. These additional simulations are (1) the constant leakage ratio, in which we considered the proportion $X_{\text{C}}^{\left(\mathrm{2}\right)}$: $R_{\text{C}}^{\left(\mathrm{2}\right)}$ in phytoplankton dpp to be a constant fraction of leakage and therefore removed the dependence of $f_{R\mathrm{2}}^{X\mathrm{2}}$ on $\mathit{\theta }/{\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}}$, and (2) the constant dpp ratio, in which we considered the fraction $X_{\text{C}}^{\left(\mathrm{2}\right)}$: $R_{\text{C}}^{\left(\mathrm{2}\right)}$ in dpp to be constant and therefore removed the dependencies of $f_{R\mathrm{2}}^{X\mathrm{2}}$ on $\mathit{\theta }/{\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}}$ and on leakage/dpp.

In the second experiment, EXP-Physics, we investigated the role of allochthonous CDOM in maintaining CDOM concentrations in the surface ocean. This experiment consisted of comparing the results of the optimized configuration with an additional simulation of intense bleaching, in which we prescribed the same photobleaching parameters for $X_{\text{C}}^{\left(\mathrm{3}\right)}$ (b_X(3) and I_X(3)) as we did for the other two CDOM pools and quadrupled the CDOM : DOC fractions in all biotic sources ($f_{Z^{\text{(MICRO)}}}^{X\mathrm{1}}$, $f_{Z^{\text{(HETEN)}}}^{X\mathrm{1}}$, $f_{P^{\text{(DIATOM)}}}^{X\mathrm{2}}$, $f_{P^{\text{(NANO)}}}^{X\mathrm{2}}$, $f_{P^{\text{(PICO)}}}^{X\mathrm{2}}$, $f_{P^{\text{(DINO)}}}^{X\mathrm{2}}$, and $f_B^{X\mathrm{3}})$. By comparing the optimized configuration with the intense-bleaching simulation, we evaluated whether an increased, locally produced CDOM could compensate for a reduced allochthonous input and the extent to which the optimized simulation was able to represent the role of $X_{\text{C}}^{\left(\mathrm{3}\right)}$, and thus of any source of allochthonous CDOM, in generating the observed a_CDOM(λ) values.

Table 2 List of simulations performed to test the findings provided by the optimized framework. The list reports all the formulations and/or parameter values which distinguish the test simulations from the optimized model configuration.

3 Results

3.1 Skill of the optimized simulation

In this section we evaluate how the optimized model configuration reproduces seasonal cycles and vertical gradients in terms of biological and optical variables. We compared the output of the optimized simulation with in situ observed data of mixing and nutrients (Sect. 3.1.1), BGC variables (Sect. 3.1.2), and IOPs (Sect. 3.1.3) and with satellite-retrieved data not used in the optimization (Sect. 3.1.4). The optimized parameter values and the quantitative model data skill metrics (R and RMSE) achieved throughout optimization can be found in Appendix A.

3.1.1  Mixing and nutrients

The temporal evolution of a climatological year (2011–2014) of the temperature in the water column is presented for the observations (Fig. 3a) and for the optimized model run (Fig. 3b). Overall, there is good agreement in terms of both the timing of the stratification and the mixing events. The seasonal cycles of nitrate (Fig. 3c vs. 3d) and phosphate (Fig. 3e vs. 3f) are well captured, along with the depth of the nitracline, which ranges from 45 m at the onset of the stratification to 60 m at the end of the year, with NO₃ and PO₄ being underestimated by approximately 25 % in the deeper layer.

Figure 3 Annually averaged vertical profiles as a function of time of (a–b) temperature, (c–d) nitrate, and (e–f) phosphate, obtained from observations at the DYFAMED site (a, c, e) and from the optimized simulation (b, d, f). The white line in the temperature panels represents the mixed-layer depth (mld) obtained in (a) with a threshold of 0.2 ^∘C for temperature (de Boyer Montégut et al., 2004) and in (b) with a threshold of $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻² for turbulent kinetic energy. The line in (c) and (d) represents the nitracline (nitrate concentration equals 2 µM).

3.1.2 Chl a, DOC, and CDOM

Observations show maximum Chl a concentrations of more than 1 mg m⁻³ in the first 30 m of the water column in March–April and the formation of a DCM in April–May that deepens down to around 50 m and gradually shoals after October (Fig. 4a). The model succeeds in capturing the important spatial features of the site's mean total Chl a (TChl a) field, such as the formation of a DCM in April–May that reaches 50 m depth in September–October (Fig. 4e). However, the simulated maximum Chl a concentrations are about 45 % smaller than the observations (Fig. 4e vs. 4a), and, therefore, the optimized simulation generally underestimates the vertically integrated Chl a. The subdivision of Chl a into size fractions follows the main features observed at BOUSSOLE, which include a larger contribution of phytoplankton larger than 2 µm to the superficial bloom in March–April (Fig. 4g and h) and the noticeable contribution of nanophytoplankton to the TChl a at the site (Antoine et al., 2006) (Fig. 4c).

Figure 4 Annually averaged vertical profiles as a function of time of log-transformed total Chl a – (a) and (e), picophytoplankton Chl a – (b) and (f), nanophytoplankton Chl a – (c) and (g), and microphytoplankton Chl a – (d) and (h) obtained from observations at the BOUSSOLE site (left panels) and from the optimized simulation (right panels). The thick line in (e) represents the depth of the DCM.

Simulated DOC and CDOM showed clear and very different seasonal dynamics (Fig. 5). After a period of winter mixing, when the concentrations of both variables were homogeneously distributed in the column, with DOC values between 38–40 µM (Fig. 5a) and CDOM concentrations at about 1.2 mg C m⁻³(Fig. 5b), superficial concentrations of both variables started to increase in March. While CDOM concentrations reached a maximum in May and decreased rapidly at the surface afterwards, DOC concentration had a maximum in June and gradually decreased until the following winter. CDOM formed a sub-superficial maximum in May–June that deepened down to 50 m in October–November; this feature was absent in DOC.

3.1.3 Inherent optical properties (IOPs)

The optimized model configuration captured important features of the mean site a_PH(450) field, which mainly follows the distribution of TChl a (Fig. 6a). After the general vertical homogeneity established in winter, the maximum values (0.05 m⁻¹) are found at the surface in spring; in summer, the maximum is located around 40–50 m, and the lowest values (0.005 m⁻¹) are found at the surface and below 75 m (Fig. 6b). Observed a_NAP(450) showed very low values (0.01 m⁻¹) (Fig. 6c) comparable to those simulated (Fig. 6d). Observed a_CDOM(450) showed clear seasonal dynamics with rather vertically homogeneous values (0.018–0.02 m⁻¹) until April and a subsurface maximum from April to September at 30–50 m depth. After September, a_CDOM(450) tended to be more homogeneous from the surface to 130 m while still being relatively high compared to winter (Fig. 6e). In the optimized simulation, a_CDOM(450) follows the distribution found for the observations and was therefore homogeneously distributed during winter (0.02 m⁻¹) and increased to values up to 0.04 m⁻¹ at the surface in April. Thereafter, the maximum values were confined to the depth of the yellow-substance maximum (YSM) formed at 30 m depth in June and deepened throughout summer down to 50 m in December (Fig. 6f). Our depth profile data showed that the simulated depth of the YSM (Fig. 6f) was only slightly shallower than that of the DCM (Fig. 6e), which is consistent with reports from Oubelkheir et al. (2007) and Xing et al. (2014) in the NW Mediterranean Sea.

3.1.4 Validation with satellite products

Satellite-retrieved Chl a, a_PH(443), and a_DG(443) from January 2011 to December 2014 at the BOUSSOLE site are compared with simulated values, averaged over the first optical depth of the water column (Fig. 7). In the 4-year time series considered, the simulated first optical depth ranged from 4.8 to 28 m. The simulated Chl a values are slightly underestimated with respect to the satellite value (bias of −0.028 mg m⁻³), especially in the period from June to October. If the mixing periods in February and March, when simulated TChl a concentrations drop to almost zero, are excluded from the comparison, Pearson's correlation coefficient is 0.55. The simulated a_PH(450), on the other hand, has a correlation coefficient with a_PH(443) of 0.69 with no bias, suggesting that the optimized values of ${\stackrel{\mathrm{‾}}{a}}_{\text{PH}}^{\ast }$ compensated for the underestimation of Chl a. Simulated a_DG(450) has a correlation coefficient with a_DG(443) of 0.40 with some degree of underestimation (bias of −0.0031 m⁻¹). The optimized simulation could not capture increments in a_DG(443) during November–December, probably because these maxima were not recorded by the BOUSSOLE in situ measurements at a monthly resolution.

Figure 5 Annually averaged vertical profiles as a function of time of (a) total dissolved organic carbon (DOC) and (b) total CDOM in the optimized simulation. The double line in (b) shows the position of the 50 and 60 µmol quanta m⁻² s⁻¹ of PAR that represent the PAR threshold for maximum bleaching of the biodegradable (b_X(1) and b_X(2)) and the semi-refractory CDOM (b_X(3)), respectively.

3.2 Dynamics of CDOM in relation to other BGC products

In this section, we analyse the results of the optimized simulation to investigate the following: (i) the vertical and seasonal dynamics of a_CDOM(λ) in relation to phytoplankton, bacteria, and DOC (Sect. 3.2.1); (ii) the mutual proportions of CDOM, phytoplankton, and detritus in the absorption budget (Sect. 3.2.2); and (iii) the ability to reproduce observed bio-optical relationships between a_CDOM(λ) and other BGC variables (Sect. 3.2.3).

3.2.1 Seasonal and vertical dynamics of a_CDOM(450)

The average 4-year depth profiles of TChl a, a_CDOM(450), bacterial carbon (B_C), and dissolved organic carbon (DOC) simulated for the BOUSSOLE site (depth <200 m) and their seasonal patterns are shown in Figs. 8 and 9. During winter (January to March), when deep convective water mixing occurred, TChl a was high (0.2–0.3 mg m³) and homogeneously distributed within the 0–75 m layer, while it rapidly decreased to values close to 0.02 mg m⁻³ below 80 m. DOC and a_CDOM(450) show vertical homogeneity in winter, with values around 40 µM and 0.02 m⁻¹, respectively (Fig. 8, winter). A sub-superficial maximum of TChl a formed in spring between 30–40 m, a a_CDOM(450) maximum formed at 30 m depth, whereas the maximum values of DOC and B_C remained within the first 10 m. There, DOC reached the highest seasonal concentration, with 66 µM (Fig. 8, spring). The sub-superficial maximum of TChl a deepened in summer to 45  m. A reduction of surface a_CDOM(450) values was observed from spring to summer, and the subsurface maximum was reinforced and deepened to 40 m. A B_C subsurface maximum was also formed slightly above the a_CDOM(450) maximum and at least 10 m above the DCM (Fig. 8, summer). Chl a-depleted surface waters and a weaker DCM persisted in fall. A maximum of a_CDOM(450) was also observed in fall but in slightly deeper waters, closer to the DCM (Fig. 8 Fall).

In terms of temporal dynamics, TChl a reached maximum concentrations (0.7 mg m⁻³) at the end of winter and gradually decreased to values of 0.05 mg m⁻³ in summer because of the depletion of nutrients due to the strong stratification of the water column (Fig. 9a). At the same time, a DCM formed at approximately 40 m depth, and its TChl a content increased throughout summer (Fig. 9b). The increase of a_CDOM(450) started at the same time as the TChl a increase but persisted much longer, when TChl a was much lower and B_C remained high (Fig. 9a). In the DCM, the maximum a_CDOM(450) values were at a similar level of maxima from March to June (Fig. 9b). While in the DCM, a_CDOM(450) and DOC followed the same temporal pattern (Fig. 9b), at the surface, the maximum of a_CDOM(450) was reached 2 months earlier than that of the DOC, which continues to accumulate during most of the spring, whereas a_CDOM(450) decreases as a result of photobleaching (Fig. 9a). The similar temporal patterns of a_CDOM(450) and B_C at the surface (Fig. 9a) and especially at the DCM (Fig. 9b) indicate a possible involvement of bacteria in CDOM production, although in our optimized simulation, semi-refractory CDOM is mainly of allochthonous origin (see Sect. 3.3.2 for a more detailed analysis).

Figure 6 Annually averaged vertical profiles as a function of time of (a–b) a_PH(450), (c–d) a_NAP(450), and (e–f) a_CDOM(450), obtained from observations at the BOUSSOLE site (a, c, e) and from the optimized simulation (b, d, f). The thick line in (f) represents the depth of the yellow-substance maximum (z_YSM) (depth of maximum a_CDOM(400) above the 1028.88 isopycnal).

3.2.2 Light absorption budget

The relative contributions of CDOM, phytoplankton, and detritus to the total non-water light absorption at 450 nm are shown for the optimized simulation at the BOUSSOLE site in Fig. 10. During the 4-year simulation, the contribution of CDOM to the light absorption budget in the entire water column down to 100 m varied between 12 % and 99 %, with an average value of 64 %. At the surface (Fig. 10a), detritus was the minor contributor to the absorption budget (1 % to 10 %). Phytoplankton was the major absorbing substance only during the bloom in March and April, whereas a_CDOM(450) dominated the rest of the year and contributed more than 50 % to the light absorption. This dominance even held in summer and fall when a_CDOM(450) was minimal as a result of photobleaching (Fig. 6f). In accordance with the results reported by Organelli et al. (2014), the surface waters of BOUSSOLE can therefore be classified as CDOM type for most of the year, except for the algal bloom in March and April. In the DCM (Fig. 10b), on the other hand, with the exception of a few short mixing periods, a_PH(450) always contributed the most to absorption but never exceeded 60 %. Detritus was still a minor component (1 % to 16 %), and a_CDOM(450) varied between a minimum of 30% in May and a maximum of 60 %–70 % in winter. Note that the depth of the YSM was consistently about 5 m shallower than the DCM (Fig. 8 in spring, summer and fall). Even though the maximum values are not exactly at the DCM, a_CDOM(450) is an important contributor to total absorption in this layer as the configuration of the optimized model is able to reproduce the formation of a CDOM subsurface maximum in spring that progressively intensifies throughout the summer, closely connected to the DCM.

3.2.3 Bio-optical relationships with a_CDOM(λ)

The overall relationship between TChl a and a_CDOM(λ), shown in Fig. 11a, reflects the covariation of these variables over vertical and seasonal scales. In fact, the optimized simulation showed minimum a_CDOM(450) values of 0.018 m⁻¹ (see also Fig. 6f) with very low TChl a concentrations, in agreement with data observed at the BOUSSOLE site. a_CDOM(450) also tended to increase with TChl a concentrations, but the variability was large, with a_CDOM(450) ranging from 0.025 to 0.07 m⁻¹ for TChl a values equal to 1 mg m⁻³ (Fig. 11a). Furthermore, in agreement with Organelli et al. (2014), data shown in Fig. 11a were generally above the relationships proposed by Morel and Gentili (2009b) and Bricaud et al. (2010), derived from a global bio-optical model in the first case and from in situ measurements collected in the southeast Pacific Ocean in the second one. The main reason is the higher-than-average contribution of CDOM for a given Chl a content in the Mediterranean Sea in contrast to other oceans. Both TChl a and a_CDOM(450) were used as observations in the optimization; therefore, the emergence of such a relationship was expected.

The optimized model configuration also allowed the emergence of bio-optical relationships with unobserved variables, such as DOC. The simulated a_CDOM(250) as a function of DOC concentration showed a fairly linear relationship (Fig. 11b). The same relationships found by Catalá et al. (2018) and Galletti et al. (2019) (grey areas and lines in Fig. 11b) were based on in situ measurements collected over the entire Mediterranean Sea in April–May of different years. Our modelled data for spring show a good agreement with their samples taken in the upper part of the water column (with a_CDOM(250) values larger than 0.9 m⁻¹), while our modelled data for fall–winter agree better with samples taken below 200 m (a_CDOM(250) values smaller than 1 m⁻¹); this was a period when the water column was more homogeneous and when the values of both DOC (Fig. 5a) and a_CDOM(λ) (Fig. 6f) in the superficial part were more similar to those in the deep part of the water column. Note that Catalá et al. (2018) and Galletti et al. (2019) reported a_CDOM(λ) values at 254 nm, while modelled data report a_CDOM(λ) values in the 250 nm waveband that range from 187.5 to 312.5 nm.

Figure 7 Time series (January 2011 to December 2014) at the BOUSSOLE site of TChl a (mg m⁻³), a_PH(450), and a_DG(450) averaged over the first optical depth (calculated as $z_{\text{EU}}/\mathrm{4.6}$). The optimized simulation (in blue) is compared to the same satellite-retrieved products at 443 nm (in black) reporting the correlation coefficient (R), the modelling efficiency (MEF), the root mean squared error (RMSE), and the average error (bias), all calculated following Stow et al. (2009). Comparable variables measured at the BOUSSOLE site are shown for reference (grey dots).

Figure 8 Seasonal vertical profiles of phytoplankton TChl a, a_CDOM(450), bacterial biomass, and DOC. Solid lines show the median value for the season of 4 years (2011 to 2014), and the coloured areas show the percentiles in 10 % steps from darker to lighter.

3.3 Sources of CDOM

In this final section of results, using a series of hypothesis-testing experiments (Table 2), we investigated the mechanisms that determine the vertical and temporal patterns of CDOM distribution at the BOUSSOLE site. We step through two aspects: the biogenic in situ production of biodegradable CDOM (EXP-Biology, Sect. 3.3.1) and the role of CDOM of allochthonous origin in maintaining CDOM concentrations in the surface ocean (EXP-Physics, Sect. 3.3.2).

3.3.1 In situ production of CDOM and synchronization with phytoplankton Chl a

EXP-Biology focuses on the physiological processes affecting CDOM production by phytoplankton, examining which is the coloured fraction of the total DOC of phytoplankton origin (Eq. 5). In the optimized simulation, the parameters for the production of $X_{\text{C}}^{\left(\mathrm{2}\right)}$ by phytoplankton ranged from 3.2 % to 8.9 % ($f_{P\text{(DIATOM)}}^{X\mathrm{2}}=\mathrm{8}\pm \mathrm{0.9}$, $f_{P\text{(NANO)}}^{X\mathrm{2}}=\mathrm{4.3}\pm \mathrm{1.6}$, $f_{P\text{(PICO)}}^{X\mathrm{2}}=\mathrm{3.9}\pm \mathrm{0.7}$), except for P^(DINO), whose value did not converge (Appendix A). These parameters, regulated by $\mathit{\theta }/{\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}}$, represent which fraction of the leakage is coloured, whereas DOC from exudation is assumed to be non-coloured. Therefore, the fraction $X_{\text{C}}^{\left(\mathrm{2}\right)}$ : $R_{\text{C}}^{\left(\mathrm{2}\right)}$ in the optimized simulation was variable in time and depth for two reasons: (1) it depends on the relative proportions of DOC from leakage or exudation, and (2) it depends on the θ to ${\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}}$ ratio of phytoplankton. Two additional simulations used a constant fraction of CDOM in the leakage flux (constant leakage ratio, Table 2) and a constant fraction of CDOM in the total dpp flux (constant dpp ratio, Table 2), setting these CDOM fractions to 6.2 % and 1.8,%, respectively, which are the average values obtained from the optimized simulation. The constant percentage values in the constant leakage ratio and the constant dpp ratio were chosen to homogenize the values of $f_{R\mathrm{2}}^{X\mathrm{2}}$ at the beginning of the year between the formulations and to highlight the vertical and temporal differences in $f_{R\mathrm{2}}^{X\mathrm{2}}$ between the three formulations.

Figure 9 Simulated annual evolution of phytoplankton TChl a, a_CDOM(450), bacterial biomass (B_C), and DOC at (a) the surface (0–9 m) and at (b) the depth of DCM. The solid lines show the average of the year from 2011 to 2014, and the coloured areas show the interannual variability.

Figure 12a shows the percentage of CDOM flux with respect to the total DOC flux ($f_{R\mathrm{2}}^{X\mathrm{2}})$ at the surface (0–9 m) and at the depth of the DCM in the three simulations of EXP-Biology. In the constant-leakage-ratio simulation, when the sub-surface Chl a maximum develops in March–April, the proportions of CDOM production differ between the surface (1 %) and the DCM (3 %) as nutrient availability begins to decrease at the surface compared to the DCM. This difference remained fairly stable until the next mixing, except for an increase to 4 % in the DCM at the end of November. In the optimized simulation, the fraction of CDOM at the surface and at the DCM also differs from March onwards, but the percentages are lower in spring and summer (0.5 % and 2 %, respectively) and gradually increase in fall and winter as light intensity decreases and phytoplankton photoacclimate accordingly. In terms of $X_{\text{C}}^{\left(\mathrm{2}\right)}$ production flux, this means that, whereas the constant-dpp-ratio simulation is the one producing more $X_{\text{C}}^{\left(\mathrm{2}\right)}$ at the surface and less $X_{\text{C}}^{\left(\mathrm{2}\right)}$ at the DCM (Fig. 12b vs. 12c), the constant-leakage-ratio simulation increases substantially CDOM production in the DCM compared to the other two (Fig. 12c), and the optimized simulation substantially decreases CDOM production at the surface (Fig. 12b). With respect to the observable variable a_CDOM(450), the optimized formulation is able to reproduce the dynamics of a_CDOM(450) at the surface during the transition from the winter productive phase to the summer better than the other two (Fig. 12d), while the constant-leakage-ratio simulation captures the dynamics of the DCM in summer quite similarly to the optimized one, suggesting that nutrient regulation predominantly determines the dynamics in the DCM (Fig. 12e).

Figure 10 Annual cycle of the relative contributions of CDOM, phytoplankton, and detritus to total non-water light absorption at 450 nm (a) in the first 0–9 m and (b) at the DCM of the BOUSSOLE site. The solid lines show the annual average from 2011 to 2014, and the coloured areas show the interannual variability.

3.3.2 The role of allochthonous CDOM and/or mediation by bacteria

In the optimized simulation, the deep inventory of $X_{\text{C}}^{\left(\mathrm{3}\right)}$ replenishes the surface during mixing and is slowly photobleached (b_X(3)=0.00119 d⁻¹ and I_X(3)=50 µE m⁻² d⁻¹, Appendix A). The optimized value of b_X(3) was much lower than the prescribed value for the labile and semi-labile pools (0.167 d⁻¹). In EXP-Physics (Fig. 13), we investigated whether maintaining the prescribed bleaching rate while increasing local production could have an equivalent effect. The optimized model configuration simulated a moderate increase in $X_{\text{C}}^{\left(\mathrm{3}\right)}$ bleaching rates towards the surface compared to the simulation of intense bleaching, where the bleaching parameters for $X_{\text{C}}^{\left(\mathrm{3}\right)}$ were set to be equal to the other two pools (Fig. 13a, blue line vs. magenta line). In the optimized simulation, the lower $X_{\text{C}}^{\left(\mathrm{3}\right)}$ degradation contributed to the minimum values of a_CDOM(450) at the surface being above 0.018–0.02 m⁻¹ throughout the year, with $X_{\text{C}}^{\left(\mathrm{3}\right)}$ contributing more than 80 % of the total pool in fall–winter and around 50 % in spring (Fig. 13b, blue line). In the intense-bleaching simulation, on the other hand, where $X_{\text{C}}^{\left(\mathrm{3}\right)}$ was bleached more intensively in the illuminated part of the water column and biota-produced CDOM was quadrupled, the minimum values of fall and winter were equally maintained by $X_{\text{C}}^{\left(\mathrm{2}\right)}$ and $X_{\text{C}}^{\left(\mathrm{3}\right)}$, but during the productive season, $X_{\text{C}}^{\left(\mathrm{1}\right)}$ and $X_{\text{C}}^{\left(\mathrm{2}\right)}$ were largely overestimated (Fig. 13b, magenta line).

The intense-bleaching simulation reproduced consistent average values of a_CDOM(450), though with stronger seasonal differences than the optimized simulation at depths close to the DCM (Fig. 13c). However, below the DCM, the contribution of allochthonous $X_{\text{C}}^{\left(\mathrm{3}\right)}$ is smaller in the intense-bleaching simulation than in the optimized simulation, and the biotically produced CDOM cannot maintain the observed values of a_CDOM(450) (Fig. 13d). These results indicate that the optimization procedure increased $X_{\text{C}}^{\left(\mathrm{3}\right)}$ concentration by reducing bleaching because even a very high (4-fold) increase in local production cannot compensate. Therefore, only allochthonous sources could balance bleaching and keep a_CDOM(450) larger than 0.02 m⁻¹ throughout the year in the upper layer. We cannot rule out the possibility that this background a_CDOM(450) is sustained in reality by bacterial production. Optimization decreased b_X(3) instead of increasing $f_B^{X\mathrm{3}}$ (2.3±1.3 %), probably because increasing $f_B^{X\mathrm{3}}$ cannot maintain the high a_CDOM(450) values observed below the DCM since bacteria are hardly present in this layer in our simulations (Fig. 8). In any case, in the optimized simulation, the above-average a_CDOM(450) values are maintained by $X_{\text{C}}^{\left(\mathrm{3}\right)}$, either by allochthonous $X_{\text{C}}^{\left(\mathrm{3}\right)}$ with a lower photobleaching rate than that produced by phytoplankton or by $X_{\text{C}}^{\left(\mathrm{3}\right)}$ produced by bacteria. Below 150 m, the only contributor to the total pool of CDOM is $X_{\text{C}}^{\left(\mathrm{3}\right)}$ in both simulations. With an initialization value for $X_{\text{C}}^{\left(\mathrm{3}\right)}$ of 1.25 mg m⁻³ and a remineralization rate of $\mathrm{3}\times {\mathrm{10}}^{-\mathrm{5}}$ d⁻¹ (Legendre et al., 2015), the values of a_CDOM(450) were consistent with those observed in BOUSSOLE (Fig. 13e) and reported in Organelli et al. (2014).

Figure 11 Variations of a_CDOM(λ) (a) at 450 nm (a_CDOM(450)) as a function of TChl a and (b) at 250 nm (a_CDOM(250)) as a function of total DOC, all simulated at the surface (0–9 m) of the BOUSSOLE site from 2011 to 2014 (dots coloured by season). The in situ TChl a and a_CDOM(450) values collected at BOUSSOLE from 2011 to 2016 are displayed in the background of (a) as grey dots. The regression lines from Catalá et al. (2018) and Galletti et al. (2019) and the ranges of their reported values in the Mediterranean Sea are shown in the background of (b).

4 Discussion

4.1 Optimization of a hydrodynamic–BGC model through a 1D fast system

Some of the uncertainties in coupled hydrodynamic–ecosystem models can be attributed to incomplete knowledge and representation of the fundamental processes being simulated (Frölicher et al., 2016) and poorly constrained parameter values (Murphy et al., 2004), especially for the BGC processes. The search for optimal parameter values by manual trial-and-error calibration or through brute-force approaches can be computationally unaffordable. Furthermore, testing parameters on an individual basis can lead to important information being missed, given parameter interactions and non-linearity of complex ecosystem models. Here we successfully applied an optimization procedure that integrates model output and observations to constrain a relatively large number of parameters in a one-dimensional configuration of the BFM. The constrained parameters contributed to improving the simulation of BGC (Chl a) and optical variables (IOPs) at the data-rich BOUSSOLE site. The use of the many observations available at this site led us to improving the simulation of bio-optical relationships that were either observed (a_CDOM(450) as a function of TChl a) or not observed directly at the site (DOC as a function of a_CDOM(250)) but reported by other studies in the Mediterranean Sea (Catalá et al., 2018; Galletti et al., 2019).

Figure 12 EXP-Biology: in situ production of CDOM by phytoplankton in the optimized simulation (blue), in the simulation of constant leakage ratio where the fraction $f_{R\mathrm{2}}^{X\mathrm{2}}$ was a constant 6.2 % of leakage (orange), and in the simulation of constant dpp ratio where the fraction $f_{R\mathrm{2}}^{X\mathrm{2}}$ was a constant 1.8 % of the total dpp (purple). (a) Annual cycle of $f_{R\mathrm{2}}^{X\mathrm{2}}$ at the surface (dashed lines) and at the depth of the DCM (solid lines) in the three simulations. Annual cycle of $X_{\text{C}}^{\left(\mathrm{2}\right)}$ production flux (b) averaged over the first 0–9 m and (c) in the DCM. Annual cycle of a_CDOM(450) averaged in (d) the first 0–9 m and (e) between 40–60 m; the dots show a_CDOM(450) values measured at the BOUSSOLE site.

Our rigorous optimization had some limits; in particular, it could not constrain all the parameters we had investigated. This was the case with the parameters for dinoflagellates (P^(DINO)), which could not be constrained because this group represents only a small part of the microphytoplankton size fraction compared to diatoms (P^(DIATOM)). This result indicates that the subset of observations we used may not contain enough information, which is a well-known limitation of formal optimization (Ward et al., 2010). It also suggests that the results of the optimization may vary if we use disaggregated observations for diatoms and dinoflagellates. The choice of the parameter list is also crucial. A common strategy is to optimize the parameters to which model output is more sensitive. Since we were interested in the processes that shape CDOM cycling, we chose to include in the optimization as many parameters implied in CDOM cycling as possible. Nevertheless, our results show that the optimization trajectories were, to some extent, driven by the parameters to which the model output was more sensitive. This was the case for the maximum Chl a : C quota (${\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}})$ of P^(NANO) and P^(PICO) (Fig. A2 in Appendix A), which drove the improvement in nano- and pico-phytoplankton Chl a (Fig. A1 in Appendix A) rapidly and intensively. All the parameters we investigated are involved in different processes. But interdependence and unidentified correlations between parameters can cause issues of poor identifiability; i.e. different combinations of parameters might yield the same result. Moreover, the optimized parameter values could compensate for processes that are either missing or poorly formulated in the model (e.g. the origin of $X_{\text{C}}^{\left(\mathrm{3}\right)}$ in surface waters, discussed below). These two aspects imply limitations of the values obtained in this work.

Figure 13 EXP-Physics: role of $X_{\mathrm{C}}^{\left(\mathrm{3}\right)}$ in the optimized simulation (blue) and in the simulation of intense bleaching, where the bleaching rate of $X_{\mathrm{C}}^{\left(\mathrm{3}\right)}$ was set to 0.167 d⁻¹ and where $f_{Z\text{(MICRO)}}^{X\mathrm{1}}$, $f_{Z\text{(HETEN)}}^{X\mathrm{1}}$, $f_B^{X\mathrm{3}}$, and $f_{\mathrm{P}\left(\mathrm{i}\right)}^{X\mathrm{2}}$ were multiplied by 4 with respect to the optimized values (pink). (a) Degradation rate of $X_{\mathrm{C}}^{\left(\mathrm{3}\right)}$ (deG_X3 + rem_X3) as a function of depth and a_CDOM(450) annual cycle averaged in depth layers between (b) 0–9 m, (c) 40–60 m, (d) 80–100 m, and (e) 150–400 m. The dots in (b–e) show a_CDOM(450) observations at the BOUSSOLE site. The inset pie charts show the composition of the CDOM inventory in terms of reactivities ($X_{\mathrm{C}}^{\left(\mathrm{3}\right)}$ in brown, $X_{\mathrm{C}}^{\left(\mathrm{2}\right)}$ in orange, and $X_{\mathrm{C}}^{\left(\mathrm{1}\right)}$ in white) at selected time points in the two simulations, intense bleaching (B) and optimized (O).

The data sets used in this work include fewer properties than the BGC model, and observations themselves are subject to errors. As a result, an excellent agreement between model and observation does not guarantee that the model's representation of unobserved properties is correct (Fennel et al., 2022). In our results, the optimized model output matched relatively well the observations provided (Figs. 4 and 5) and the same variables remotely retrieved (Fig. 7) while still being more consistent with established knowledge of processes at the BOUSSOLE site (Figs. 6 and 11) not considered in the optimization. As the results of the optimization may vary depending on how many observation types are used (Wang et al., 2020), an increasing breadth of observables will allow us to continually evaluate our model and predictions. High-frequency data collected by the BOUSSOLE project (every 15 min during daylight and at hourly intervals at night), such as Chl a fluorescence, bb(λ) of total particles, or remote-sensing reflectance above sea water, are valuable additional sources of information for further analysis. We have not tested the influence of the sampling frequency on the optimization performance. Using monthly mean profiles has been reported to be sufficient to contribute to the optimization of vertical Chl a structure and to optimize biological parameters (Shu et al., 2022; Bisson et al., 2021). But assimilating high-frequency data in a 1D model has been proven to be beneficial for inferring phytoplankton dynamics when compared to assimilating observations that are more limited in time (Kaufman et al., 2018).

CDOM-cycling models are useful only if they are flexible enough (i.e. portable) to provide accurate predictions for a wide range of conditions occurring in marine ecosystems. In that case, they can be integrated into larger 3D models of BGC cycles. The possibility of applying parameter optimization techniques in computationally efficient 1D models before using the resulting parameters in the 3D version seems advantageous. Furthermore, parameters optimized with surface data only (i.e. satellite retrieved) may present inconsistencies when simulating underwater biological fields (Wang et al., 2020). This is particularly problematic in oligotrophic regions where the deep Chl a maximum (DCM) is relatively deep and poorly correlated with surface satellite observations (Cullen, 2015; Fennel and Boss, 2003). Optimizing parameter values at intensively sampled 1D sites can overcome this problem and improve the simulation of both vertical and seasonal variability. Further analysis will be required to determine whether the optimal parameter sets from the 1D version should be applied directly in the 3D model. This decision should balance the risk of a lower agreement between model and data in 3D than in the 1D counterpart, with the opportunity to preserve key features presented in this analysis. These features include seasonal dynamics and vertical gradients and the higher-than-average CDOM light absorption coefficients (Sect. 3.1 and 3.2). Since the approximation of a 1D configuration is reasonable at the BOUSSOLE site, we do not expect that our proposed parameter values are compensating for the absence of lateral transport. However, considering that biogeochemical optimization can be influenced by specific physical forcing and possible biases in circulation dynamics (Pasquier et al., 2023), further analysis is needed to upscale the present 1D results to the 3D Mediterranean domain, where lateral contributions and vertical mixing are different and can interact differently with the optimized biogeochemical processes.

4.2 Hypothesis testing in a 1D virtual laboratory: CDOM-cycling numerical experiments

Our results provided new insights into the in situ production of semi-labile CDOM by phytoplankton. We simulated a clear seasonal and vertical structure of a_CDOM(λ). This is in better agreement with observations if the proportion of coloured material in the dpp is modelled to depend on nutrient and light constraints (Fig. 12). We have proposed a simplified parameterization to account for the effects of light and nutrients on CDOM production. Nutrient limitation in relation to phytoplankton CDOM production acts through the relative proportion between leakage and dpp and seems to be the main cause of the difference in CDOM production between surface and DCM. This could reflect the fact that, when growth is limited due to low inorganic nutrient availability, the photosynthetic excess produces mainly non-coloured carbohydrates as severe N and P limitation favours the liberation of carbohydrates that do not contain these elements (Myklestad, 1995), or, simply, both higher CDOM production rates and nutrient availability occur simultaneously in the DCM (Bushaw et al., 1996; Stedmon et al., 2007; Zhang et al., 2011a). Either way, the CDOM production dependency on nutrient availability helped to reproduce the observation that increased CDOM concentrations are found at the depth of the Chl a maximum (Coble et al. 1998).

Light limitation in relation to phytoplankton CDOM production acts through the elemental ratio of Chl a : C relative to the maximum value and appears to be the main cause of the seasonality of CDOM at the surface, although it contributes to amplifying the differences between surface and DCM in CDOM production. As with nutrients, this does not necessarily represent a causal mechanism but a general pattern of high-light, low-CDOM waters. Passive leakage by diffusion through the cell membrane is controlled by cell permeability, the surface-area-to-volume ratio of the cell, and the intracellular concentration of the compound in question (Bjørrisen, 1988). Thus, our results could indeed be representative of a lower concentration of coloured material in the dpp of low-pigmented cells. However, given the inverse covariation between photosynthetic and photoprotective pigments, surface communities exposed to high light intensities that have a low θ may contain a higher proportion of photoprotective pigments. This, in principle, could contribute to the leakage of coloured molecules (e.g. carotenoids). More likely than lower CDOM production under high light, our results could be representative of increased CDOM destruction and a larger or different impact of photobleaching than assumed in our simulations. Photobleaching was likely the main cause of CDOM destruction in the surface layer from spring to summer at the BOUSSOLE site (Organelli et al., 2014), as generally observed for other oceanic waters (Coble, 2007; Bracchini et al., 2010). However, it has been found that the photochemical lability of CDOM is not necessarily related to the highest solar irradiance but depends on the accumulated light exposure (Yamashita et al., 2013; Swan et al., 2012). This may explain why the modulation of CDOM production by the photoacclimation of phytoplankton (through θ) helped to indirectly improve the effect of light intensity on CDOM concentrations.

Besides biological production in situ, our results evidence the relevance of CDOM of allochthonous origin in keeping a_CDOM(450) above 0.018–0.02 m⁻¹ throughout the year in the upper layer. It cannot be excluded that this background value of a_CDOM(450) is maintained by processes in the real system that are poorly represented in the model (bacterial production) or not represented at all (other sources of allochthonous CDOM such as atmospheric deposition or advection). Nevertheless, the higher-than-average a_CDOM(450) values at the surface are supported by semi-refractory CDOM (Fig. 13b), consistent with the role of atmospheric deposition in providing higher-than-average CDOM concentrations in the basin (Galletti et al., 2019). In the water column, we lack observations on the reactivity of the CDOM pool we simulated. However, our simulations mimic the results reported by Bracchini et al. (2010) in that the refractory fraction of the CDOM dominates the water below the DCM (Fig. 13d). We cannot know whether such refractory CDOM is in fact newly produced by bacteria or whether it is old CDOM vertically transported from the deep reservoir.

4.3 Conclusions and future developments

Adding optical components to BGC models greatly increases the amount of data that can be used for model validation and calibration (IOCCG, 2020). In this work, we have shown how the integration of models and observations through an optimization method provides insights into the dynamics of light-absorbing BGC products such as CDOM. Given its strong absorption of UV and blue light, the accurate simulation of CDOM cycling and its main dynamics has major implications. The potential of using an observable proxy, such as a_CDOM(λ), for the retrieval of BGC products (e.g. DOC) can be significant, e.g. for carbon cycle studies and budgeting. Furthermore, the accurate simulation of the CDOM absorption band in the blue part of the spectrum helps in deciphering the mechanisms underlying the optical properties of the basin, especially in coastal and shelf waters and semi-enclosed seas such as the Mediterranean Sea. Future model developments may aim to simulate the fraction of CDOM that is able to fluoresce, i.e. fluorescent DOM, and therefore extend the points of convergence between observations and model results. Such an opportunity would be advantaged by fitting DOM fluorescence sensors on board BGC-Argo floats.

The results of the in situ production and photobleaching experiments contributed to our understanding of the main drivers of the cycling of CDOM compounds in the 1D configuration of the BFM model. These results highlight the usefulness of the 1D fast system as a virtual laboratory for hypothesis testing regarding CDOM cycling, even when proposing relatively simple parameterizations to include the effect of light and nutrients in the production of chromophoric DOM. By testing such parameterizations in the 1D configuration, we have demonstrated that they both help to represent the vertical and seasonal dynamics of CDOM at the site. They can also guide the development of further mechanistic models aimed at finding out what drives CDOM production by phytoplankton, even though observations on the relative importance of the different processes are still limited.

Appendix A: Calibrated parameters and insights from optimization trajectories

Twenty-five BGC and optical parameters, gathered in Table A1, were constrained (i.e. we found optimal parameter estimates with relatively small uncertainties) against observations of pico-, nano-, and microphytoplankton Chl a concentrations, a_PH(450), a_NAP(450), a_CDOM(450), and $S_{\text{CDOM}}^{\mathrm{350}-\mathrm{500}}$ collected at the BOUSSOLE site from 2011 to 2014. The variation in correlation coefficients (R) between the model and observations, as well as the root mean squared errors (RMSE) for all variables included in the multi-objective function (i.e. lnLikelihood), are shown in Fig. A1. All variables showed increased R and decreased RMSE compared to the observations along the optimization. The difference in lnLikelihood relative to the initial one reached 0.01 at generation 115 and 0.005 at generation 144; so, although the procedure was continued, all R's stop improving (vary by less than 0.02) after about 150 generations (nanophytoplankton Chl a was the last).

Interestingly, the improvements in the variables (RMSE reduction and R increase) and the convergences of the parameters (Fig. A2) were not homogeneous during the calibration process. Instead, the results show fast- and late-moving parameters. During the first 20 generations the main improvement was in the pico- and nano-Chl a (Fig. A1a, b), which corresponded with early constrained values of ${\mathit{\theta }}_{\text{Chl}}^{\mathrm{0}}$ in P^(PICO) and P^(NANO) (Fig. A2h, m). Not all independent parameters could be well constrained. Several parameters did not converge after 180 generations. First, there was the spectral slope of the labile CDOM ($S_{X\left(\mathrm{1}\right)}^{\mathrm{350}-\mathrm{500}})$. This is explained by the fact that, in our model configuration, labile CDOM is consumed as soon as it is produced; it does not accumulate, and its spectral slope has, in practice, a very limited influence on the spectral slope of the whole CDOM pool ($S_{\text{CDOM}}^{\mathrm{350}-\mathrm{500}})$. Second, none of the parameters of dinoflagellates (P^(DINO)) converged. This is because we included an aggregate value of microphytoplanktonic Chl a in the multi-objective function (micro-Chl a), and given the preponderance of diatoms (P^(DIATOM)) in this size class, varying P^(DIATOM) parameters caused all changes in microphytoplankton concentration, and P^(DINO) parameters had very limited influence. Two replicates of the optimization experiment were run in order to examine the robustness of the optimization procedure. Each replicate used different randomly generated seeds and therefore started from a different distribution of parameter values in the initial population. Despite slightly different optimization trajectories of the parameter values in the first 100 generations, the parameters that converged (excluding those of P^(DINO) and $S_{X\left(\mathrm{1}\right)}^{\mathrm{350}-\mathrm{500}})$ did so to values that were between 87 %–113 % of the value reached in the main optimization experiment (Fig. A3).

Table A1 List of parameters considered for optimization, with an indication of the numerical range given at the beginning of the optimization and the mean value and range obtained at the end (generation 180).

Figure A1 Optimization trajectories of the fit metrics: correlation coefficient (R, grey, left axis) and root mean square error (RMSE, blue, right axis). The thick lines show the median of the distribution of the metric, and the coloured areas show the percentiles (5 %–95 %). The vertical arrows show the generation at which the range of model vs. observed R is less than 0.02.

Figure A2 Optimization trajectories of parameter values: phytoplankton parameters are shown in shades of green, parameters related to CDOM are shown in shades of orange, and the only detritus parameter is shown in grey. The thick lines show the median of the distribution of parameter values, and the coloured areas show the percentiles (5 %–95 %) (see Table A1 for the exact values and ranges at convergence).

Figure A3 Trajectories of parameter median values in three optimization replicates: (black) optimization experiment shown in Fig. A2 and used in the main text simulations, (orange) replicate 1 (165 generations), and (purple) replicate 2 (164 generations). Note that the scale of the y axis is zoomed in with respect to Fig. A2 to highlight differences among the median trajectories.