Merging bio-optical data from Biogeochemical-Argo floats and models in marine biogeochemistry

New autonomous robotic platforms for observing the ocean, i.e. Biogeochemical-Argo (BGC-Argo) floats, have drastically increased the number of vertical profiles of irradiance, photosynthetically available radiation (PAR), and algal chlorophyll concentrations around the globe independent of the season. Such data may therefore be a fruitful resource to improve performances of numerical models for marine biogeochemistry. Here we present a work that integrates 1314 vertical profiles of PAR acquired by 31 BGC-Argo floats operated in the Mediterranean Sea between 2012 and 2016 into a one-dimensional model to simulate the vertical and temporal variability of algal chlorophyll concentrations. The model was initially forced with PAR measurements to assess its skill when using quality-controlled light profiles, and subsequently with a number of alternative bio-optical models to analyse the model capability when light observations are not available. Model outputs were evaluated against co-located chlorophyll profiles measured by BGC-Argo floats. Results highlight that the data-driven model is able to reproduce the spatial and temporal variability of deep chlorophyll maxima depth observed at a number of Mediterranean sites well. Further, we illustrate the key role of PAR and vertical mixing in shaping the vertical dynamics of primary producers in the Mediterranean Sea. The comparison of alternative bio-optical models identifies the best simple one to be used, and suggests that model simulations benefit from considering the diel cycle.

1 Introduction

In most biogeochemical models the description of optics is generally (over)simplified. The integration of more complex optical models, where inherent and apparent optical properties (IOPs and AOPs respectively) are already included as model state variables (Fujii et al., 2007), therefore constitutes one of the necessary improvements. The research community is emphasizing the importance of merging different methods in order to improve the skill of numerical models, such as the assimilation of remote-sensing data or the use of in situ data for both initialization and validation purposes. Until recently, the use of the latter was especially critical due the scarcity of observations; however the emergence of autonomous robotic platforms such Biogeochemical-Argo floats (hereafter BGC-Argo) helped reduce the gap in bio-optical measurements acquired around the globe, regardless of the season.

The introduction of BGC-Argo floats has led to a drastic increase in the number of radiometric measurements in the Mediterranean Sea, such as downward planar irradiance and photosynthetically available radiation (PAR), for which specifically developed quality control procedures and refined sensor calibration (Organelli et al., 2016a, 2017a) have made their use widespread (Organelli et al., 2017b; Wojtasiewicz et al., 2018; Gerbi et al., 2016; Leymarie et al., 2018). BGC-Argo can therefore be an important source of high vertical spatial and temporal resolution data that can be integrated in the calibration and tuning of bio-optical numerical models for understanding marine biogeochemistry processes.

No studies have so far tried to assimilate radiometric quantities into numerical models to improve the simulation of chlorophyll dynamics in the Mediterranean Sea and to investigate the causes of the vertical, spatial, and temporal variability of zonal gradients. Assimilating radiometric data could prove to be more robust than chlorophyll assimilation as a result of a more accurate uncertainty characterization of optical measurements (Dowd et al., 2014; Organelli et al., 2016a) compared to other biogeochemical variables, such as fluorescence-derived chlorophyll.

Specific studies are required to demonstrate to what extent the assimilation of radiometric data can improve the model skill in simulating key biogeochemical variables (e.g. chlorophyll, nutrients, primary productivity). In this paper we develop a one-dimensional (1-D) model that assimilates PAR profiles acquired by BGC-Argo floats in order to replicate the vertical and temporal dynamics of phytoplankton chlorophyll concentrations. As a first modelling attempt, the exploration is carried out with a “voxel” approach, where light and mixing conditions were replicated from data available from floats. We analyse and validate model performances through a comparison of model outputs with the high number of co-located vertical profiles of chlorophyll concentrations (derived from fluorescence) measured by BGC-Argo floats. In particular, such analysis allows us to study some of the drivers modulating the deep chlorophyll maximum (DCM) depth in stratified conditions. Subsequently, we test different mixing and bio-optical models that simulate downward irradiance and evaluate their skills in order to estimate how well they perform compared to in situ measurements of PAR. The paper is organized as follows: in the Methods section, the Mediterranean Sea BGC-Argo float network and the model configurations are presented. In the Results and discussion section, we analyse the 1-D biogeochemical simulations and their sensitivity according to the objectives of the work. General remarks are illustrated in the Conclusions section.

2 Methods

2.1 BGC-Argo float data

The Mediterranean Sea BGC-Argo array operating in the period 2012–2016 (Fig. 1) was composed of 31 floats that acquired 1314 vertical profiles of temperature (T) and salinity (S), chlorophyll a concentration (Chl, mg m⁻³), derived from fluorescence measurements between 0 and 1000 m (see Organelli et al., 2017b; Roesler et al., 2017), and radiometric quantities, such as downward planar irradiance (E_d, µW cm⁻² nm⁻¹), at three different wavelengths (λ=380, 412, and 490 nm) and photosynthetically available radiation (PAR, µmol photons m⁻² s⁻¹) integrated between 400 and 700 nm (Kirk, 1994). Radiometric measurements were obtained in the upper 250 m, with vertical resolution of 1 m between 10 and 250 m and 0.20 m between 0 and 10 m. All profiles were acquired around local noon.

Figure 1 Spatial distribution of BGC-Argo float profiles superimposed on subbasin division used in the Mediterranean Copernicus Marine Environment Monitoring Service (CMEMS) system.

The quality control (QC) procedure of radiometric profiles was specifically designed to identify and remove the dark signal, atmospheric clouds, and wave focusing at the surface (Organelli et al., 2016a). Note that the operational definition of PAR used by the BGC-Argo community takes into consideration the planar irradiance E_d rather than the scalar one E_o, therefore differing from its theoretical definition and leading to an underestimation of E_o values by 30 % or more (Mobley et al., 2010). The scalar values of PAR were thus derived according to Baird et al. (2016), although the correction related to the irradiance scattering was neglected due to the lack of information on IOPs (see Sect. S1 in the Supplement).

Vertical profiles of Chl concentration were quality-controlled according to the procedure of the international BGC-Argo programme that removes spikes and corrects for non-zero deep values and non-photochemical quenching at the surface (Schmechtig et al., 2016; Organelli et al., 2017b). Due to a factory calibration bias for WET Labs ECO series Chl fluorometers, Chl concentrations were corrected by a factor of 0.5 (Roesler et al., 2017; Organelli et al., 2017a, b; Barbieux et al., 2018).

All the data used in this study are freely available and compiled into the database published by Organelli et al. (2017b). Seven variables (T, S, Chl, E_d(380), E_d(412), E_d(490), PAR) were vertically interpolated to a resolution of 1 m in the upper 400 m. Finally, we partitioned the profiles geographically into 13 (out of 16) subbasins (Fig. 1), with the majority of profiles located in the northwestern Mediterranean (NWM, 332 profiles), followed by northern Ionian (ION3, 172 profiles) and southern Tyrrhenian (TYR2, 162 profiles) seas. No data were available for the southwestern Ionian (ION1) and the eastern Levantine (LEV4) seas and only one profile was present in the northern Adriatic Sea (ADR1), as well as in the western Levantine Sea (LEV1). The WMO code specification for each BGC-Argo float (along with their operational periods) is provided in the Sect. S2 in the Supplement.

2.2 1-D biogeochemical model

Biogeochemical processes have been simulated according to the voxel approach (“volume element with biological content and processes”; Kohlmeier and Ebenhöh, 2009), discretized along the vertical direction in order to resolve vertical irradiance attenuation and nutrient gradients. Each voxel replicated light and mixing conditions according to the trajectory and measurements of the corresponding BGC-Argo float, thus simulating a pseudo-Lagrangian experiment. No exchanges of mass between voxel and the surrounding field have been considered, which implies smaller mass exchanges due to horizontal diffusion and baroclinic components of the (upper ocean) advection field compared to vertical processes and biogeochemical dynamics. Conversely, a voxel exchanges heat with the atmosphere and receives light in accordance with its moving position. Such an approach, similar to the one adopted by Kohlmeier and Ebenhöh (2009), has already been successfully applied by Mignot et al. (2018) in order to analyse BGC-Argo floats in the North Atlantic.

Furthermore, it is assumed that major biogeochemical transformations can be described by the Biogeochemical Flux Model (BFM) parametrizations, properly driven by a bio-optical model, which has been validated by contrasting model results and experimental data. The model is formulated through a system of partial differential equations:

$$\begin{array}{ll}{\displaystyle }{\partial }_t{C}_i(z,t)& {\displaystyle }={\partial }_z\left[D_{\mathrm{v}}\right(z,t\left){\partial }_z{C}_i\right(z,t\left)\right]+v_{\mathrm{sink},i}{\partial }_z{C}_i(z,t)\\ \text{(1)}& {\displaystyle }& {\displaystyle }+{\mathrm{BFM}}_i(T,S,\mathrm{PAR},\stackrel{\mathrm{‾}}{C}(z,t\left)\right),\end{array}$$

where C_i is the ith biogeochemical tracer simulated (i=1, 50), D_v is the vertical eddy diffusivity derived with the vertical mixing model described in Sect. 2.2.1, v_sink is the sinking velocity, and BFM_i is the reaction term corresponding to the tracer C_i. T, S, and PAR are the data measured by BGC-Argo floats. Since the surfacing of BGC-Argo floats is programmed at around local noon, the variability related to diurnal variation in solar irradiance is taken into consideration according to Kirk (1994).

Dutkiewicz et al. (2015)

Table 1 Model configurations considered in the present work. All simulations include diurnal variability except the two cases with continuous light (CL1 and CL2), which use 24 h averaged irradiance values.

Table 2 Parameters derived for optical models using BGC-Argo float data. For each version of OPT2 the regression is performed in the depth range indicated by z_max. Data points are averaged for layers of 15 m thickness.

The biogeochemical model BFM (Vichi et al., 2013) is a biomass-based numerical model that simulates the biogeochemical fluxes of carbon, phosphorus, nitrogen, silicon, and oxygen, characterizing the lower trophic level (producers, consumers, and recyclers) of the marine ecosystem. Its application is based on the coupled transport–biogeochemical model OGSTM-BFM (Lazzari et al., 2012, 2016). It includes four phytoplankton functional types (diatoms, nanoflagellates, picophytoplankton, and dinoflagellates), carnivorous and omnivorous mesozooplankton, bacteria, heterotrophic nanoflagellates, and microzooplankton. Each variable is described in terms of internal carbon, phosphorus, and nitrogen concentrations. Particulate and dissolved organic matter are also included, with the latter partitioned in labile, semi-labile, and semi-refractory phases (for parameters' specifications see Sect. S3 in the Supplement). The present study is focused mainly on Chl, reserving for future analysis (according to data availability and optical model complexity) a study of plankton functional type (PFT) resource competition dynamics (Ryabov and Blasius, 2011, 2014).

In particular, the Mediterranean Sea Monitoring and Forecasting Centre (Med-MFC) has operatively produced analyses, forecasts, and reanalyses of a series of biogeochemical state variables (e.g. Chl, nutrients, pCO₂) for Copernicus Marine Environment Monitoring Services (CMEMS) since 2015 using the MedBFM model (Lazzari et al., 2010, 2012, 2016), which embeds the OGSTM-BFM and assimilates surface Chl from satellite observations (Teruzzi et al., 2014, 2018).

We tested a total of 17 classes of simulations (summarized in Tables 1 and 2).

• In the first set of simulations, the biogeochemical model was forced with PAR from BGC-Argo floats. Experimental values of temperature and density (computed from float profiles) were also taken into consideration. A simulation for each of the BGC-Argo float trajectories was performed with this set-up, hereafter abbreviated as REF.

• Four additional sets of simulations were performed on the REF configuration by applying different values of vertical eddy diffusivity coefficients (MLD1, MLD2, MLD3, and MLD4) in order to assess uncertainties due to different vertical diffusion parametrizations.

• Six additional sets of simulations were performed by forcing the biogeochemical model with PAR obtained by alternative bio-optical parametrizations (OPT1, OPT2a, b, c, d), one of which (OPT3) also considers the current modelling approach used in the Med-MFC. In this way, the possibility of using biogeochemical models in the absence of PAR measurements was assessed.

• A set of simulations was devoted to understanding the impact of using a constant light approximation (CL1 and CL2 configurations) rather than following the diurnal light variation on Chl distribution.

• Furthermore, we evaluated the impact on light propagation due to coloured phytoplankton degradation products, i.e. coloured dissolved organic matter (CDOM) (OPT4a, b, c, d and OPT5).

Initial conditions for all biogeochemical variables of BFM are provided by the CMEMS reanalysis of Mediterranean Sea biogeochemistry (period 1999–2015; Teruzzi et al., 2014) produced by the MedBFM model system. The initialization profiles for our 1-D configuration are extracted from the MedBFM model output array, taking the nearest model point to the BGC-Argo position in time and space.

The simulations' timescale corresponds to a typical BGC-Argo time-series length during the period 2012–2016, i.e. 11 months on average, with a vertical resolution of 1 m. After being initialized, the model evolves without further assimilation of biogeochemical data from the 3-D configuration.

Vertical eddy diffusivity coefficient profiles D_v(z) are represented as Gaussian-shaped functions, thus allowing a gradual increase in vertical mixing through the pycnocline. Approaches and impacts of using different parametrizations to reconstruct mixing along the water column are shown and discussed in Sect. 2.2.1.

2.2.1 Vertical mixing models

Vertical mixing is estimated from potential density (obtained from temperature and salinity data from floats) along the water column. Vertical eddy diffusivity coefficients (D_v) are defined as Gaussian-shaped functions in the form of

$$\begin{array}{}\text{(2)}& D_{\mathrm{v}}=D_{\mathrm{v}}^{\mathrm{MLD}}{e}^{-\mathrm{0.5}{\left(\frac{z}{\left({\mathit{\sigma }}^{*}\mathrm{MLD}\right)}\right)}^{\mathrm{2}}}+D_{\mathrm{v}}^{\mathrm{background}}.\end{array}$$

σ was identified after an initial tuning procedure and equals 0.3 in all simulations. Values in the REF model are equal to $D_{\mathrm{v}}^{\mathrm{MLD}}=\mathrm{1.0}$ m² s⁻¹ and $D_{\mathrm{v}}^{\mathrm{background}}$ = 10⁻⁴ m² s⁻¹.

The mixed-layer depth (MLD) was defined with the density criterion at the threshold value (de Boyer Montégut et al., 2004; D'Ortenzio and Prieur, 2012):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathit{\rho }}_{\mathit{\theta }}=\left|{\mathit{\rho }}_{\mathit{\theta }}\right(\mathrm{10}\mathrm{m})-{\mathit{\rho }}_{\mathit{\theta }}(z\left)\right|=\mathrm{0.03}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}.\end{array}$$

In simulations MLD1, MLD2, MLD3, and MLD4, $D_{\mathrm{v}}^{\mathrm{background}}$ values were perturbed by 2 orders of magnitude (from 10⁻⁶ to 10⁻⁴ m² s⁻¹) in order to estimate the impact such variations have on modelled Chl profile shapes compared to measured ones (see Table 1).

2.2.2 Bio-optical models

Alternative parametrizations to measured PAR profiles were used in models OPT1, OPT2abcd, OPT3, OPT4abc, and OPT5. They differ in methods used to evaluate the diffuse attenuation coefficient K_d(PAR), which is parametrized as a function of Chl concentration rather than being directly calculated from BGC-Argo irradiance data (see Tables 1 and 2).

OPT1 uses the relationship obtained by a statistical analysis performed by Riley and Conover (1956); Riley (1975):

$$\begin{array}{}\text{(4)}& K_{\mathrm{d}}\left(\mathrm{PAR}\right)=\mathrm{0.04}+\mathrm{0.0088}\left[\mathrm{Chl}\right]+\mathrm{0.054}[\mathrm{Chl}{]}^{\frac{\mathrm{2}}{\mathrm{3}}}.\end{array}$$

In OPT2 models, statistical regressions were carried out between K_d(PAR) and Chl measured by BGC-Argo floats at four different depth ranges: 150, 75, 45 and 30 m (OPT2a to OPT2d; see Table 2 for details):

$$\begin{array}{}\text{(5)}& K_{\mathrm{d}}\left(\mathrm{PAR}\right)=a[\mathrm{Chl}{]}^b+c.\end{array}$$

a and c represent regression coefficients and b the exponent (values reported in Table 2). Confidence intervals were calculated with a Student's two-sided t test, where the significance level α was set equal to 0.05. Diffuse attenuation coefficients K_d(PAR) were calculated for PAR measured by BGC-Argo floats as the local slope of the natural logarithm of downwelling irradiance for layers of 15 m thickness for the euphotic depth range, which corresponds to an attenuation of downward planar irradiance to 1 % of the subsurface value (Kirk, 1994).

Albeit the regression in the upper 30 m of the water column showed the highest correlation, all four bio-optical models were considered and adopted in simulations OPT2a, b, c, and d (Table 2).

In model OPT3, based on the BGC-Argo data set, K_d(PAR) is calculated for the first optical depth (Morel, 1988) and the layer of interest for satellite remote sensing (Gordon and McCluney, 1975), and then adopted as a constant parameter for the entire water column. Such a light extinction definition has also been used in the 3-D version of the OGSTM-BFM model, which integrates K_d(490) data from satellite sensors as the external optical forcing in the exponential formulation of downwelling irradiance (for more details see Lazzari et al., 2012, Sect. 2.2.3).

OPT4 and OPT5 models include CDOM dynamics, as in the Mediterranean Sea the latter can absorb more than 50 % of blue light (Organelli et al., 2014; Morel and Gentili, 2009), thus significantly impacting its attenuation along the water column. OPT4 assumes that CDOM is correlated to Chl production (Organelli et al., 2014) and that the light attenuation is therefore affected by a progressive accumulation of the latter (“dead” Chl, initialized at zero concentration). In OPT4, accumulation is compensated for by a linear decay set at different e-folding characteristic times: 1 d (OPT4a), 1 week (OPT4b), and 1 month (OPT4c).

OPT5 implemented a formulation of CDOM as described in Dutkiewicz et al. (2015): a 2 % fraction of all dissolved organic matter (DOM) fluxes is directed to CDOM, including both temperature-related decay and a photodegradation term based on PAR (Bissett et al., 1999). Given the mono-spectral nature of the current description of light, the attenuation of CDOM on PAR is computed by averaging the exponential law of CDOM absorption (Bricaud et al., 1981) in the visible range. Additional investigations are provided in Sect. 3.3 to discuss CDOM dynamics along the water column.

2.3 Statistical analysis

According to the work's objectives, four classes of simulations were considered, which correspond to the following subsections: the reference simulation, a subset with perturbed vertical mixing models, tests with different optical configurations, and a last group of additional analyses involving CDOM description and diurnal variability. Outputs are validated qualitatively and quantitatively in terms of profile shapes and the deep chlorophyll maximum (DCM) depth. The DCM definition is based on the absolute maximum of Chl, excluding results of DCM shallower than 40 m or deeper than 200 m, as well as the ones with concentrations lower than 0.1 mg m⁻³. All results, for both model and BGC-Argo floats, are averaged on a weekly basis. Model outputs are compared with a match-up shown as target and Taylor diagrams (Jolliff et al., 2009). The former evaluates results with root-mean-square distance (RMSD) as the main statistical parameter, which was calculated following Eq. (6):

$$\begin{array}{}\text{(6)}& \mathrm{RMSD}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\mathrm{\Sigma }}_{i=\mathrm{1}}^n(m_i-o_i{)}^{\mathrm{2}}},\end{array}$$

where n is the number of data, m are the model values, and o are the observables.

Due to the various sources of possible uncertainties in the fluorescence-to-Chl conversion of BGC-Argo profiles, we chose to focus our study on DCM depth rather than DCM magnitude. This is additionally justified with the BFM statistical sensitivity analyses (see Sect. S4 in the Supplement), which considered DCM width, DCM magnitude, and surface Chl. Results indicate that DCM depth is the most effective feature the model is able to reproduce.

3 Results and discussion

3.1 Reference simulation

The assimilation of PAR profiles helped to accurately estimate the DCM depth (Fig. 2). Measured and modelled DCM depth showed high correlation (r=0.83, p value < 0.005). Both model and measurements indicate that DCM depth varies typically between 50 and 70 m in western areas (ALB, SWM1, SWM2, NWM, TYR) and is generally deeper in eastern areas (ADR2, ION3, LEV2, LEV3), between 100 and 140 m. The model tends to slightly underestimate the DCM depth variability (Fig. 2, regression slope = 0.81<1): the deepest simulated DCMs are around 125 m depth, whereas data from floats reach 140 m (e.g. lovbio018c).

Figure 2 Match-up diagram comparing DCM depth obtained from BGC-Argo float data versus REF model results. Each dot corresponds to a weekly profile. The red line depicts the linear regression between data and model values, defined by its slope and intercept (Y-int) shown in the box. Units of RMSD, bias, and Y-int are in metres. The correlation coefficient r is significant, with p value < 0.005. The bottom sub-figure shows the residuals' histogram.

Chl patterns display high variability at both temporal and vertical scales (Figs. 3 to 6). The subsurface Chl pattern is formed by patchy structures and it is generally deepening eastward during stratification periods (Fig. 6). BGC-Argo observations indicate that DCM is further eroded by vertical mixing occurring predominantly in autumn and early winter. At the ocean surface, the increase in Chl is triggered by rather shallow mixing (≤75 m). Simulations provide an adequate reproduction of the Chl mixing timing and therefore of the DCM erosion. By comparing Hovmöller maps of all 31 floats (considering both depth and time variability) for measured and simulated Chl, a significant average correlation coefficient (r) of 0.75 is obtained: such a result quantitatively confirms that the alternation of mixing and stratification phases, as seen from BGC-Argo Chl measurements, is well reproduced.

Figure 3 Hovmöller diagrams of BGC-Argo float lovbio067c (WMO code 6901649) comparing measured results and simulated ones (REF). The six-image composite is organized as follows: panels (a), (b), and (c) show PAR, vertical eddy diffusivity, and the float trajectory; panels (d), (e), and (f) show Chl derived from fluorescence measurements, simulated Chl, and phosphate. The thick black–white line indicates the depth where PAR equals 5.8 µmol photons m⁻² s⁻¹ (Mignot et al., 2014). The number in parentheses in modelled Chl indicates point-by-point correlation with BGC-Argo float Chl.

Figure 4 As Fig. 3 but for the BGC-Argo float lovbio035b (WMO code 6901511).

Figure 5 As Fig. 3 but for BGC-Argo float lovbio016c (WMO code 6901510).

Figure 6 As Fig. 3 but for BGC-Argo float lovbio066d (WMO code 6901655).

Simulated Chl also reproduces episodic signals, such as Chl deepening due to specific mixing events. For example, a mixing event in the NWM subbasin, reaching approximately 200 m depth during winter 2015, triggers an intrusion of Chl (0.2 mg m⁻³) in deeper layers consistent with BGC-Argo float measurements (float lovbio067c, Fig. 3). Similar dynamics are reproduced in winter 2014 (Fig. 4) for the lovbio035b float drifting from NWM toward the ALB subbasin.

Considering float trajectories, two kinds of situations are possible: the BGC-Argo float trajectory is relatively stationary in the deployment area (Figs. 3, 5, and 6), or the float migrates extensively by following a given water mass (Fig. 4). Results confirm that in the second case, when lateral advection processes could play an important role in the float dynamics, the applied approach also allows an adequate representation of measured Chl patterns. The present multi-float simulation, however, does not include trajectories comprising both west and east Mediterranean basins, where strong gradients between deep water nutrient inventories could invalidate the approach. In such cases, nudging or more sophisticated techniques would be required (Kohlmeier and Ebenhöh, 2009).

REF results further demonstrate that irradiance along the water column is the driving mechanism controlling DCM depth in addition to mixing, which is proven by a significant correlation between DCM and euphotic depths for both measured and simulated Chl (Fig. 7a). Such findings were already established in Mignot et al. (2014), indicating that the DCM is located at a fixed PAR value, oscillating near the 5.8 µmol photons m⁻² s⁻¹ isolume (Fig. 7b, blue line). Data and model outputs in the present study show a higher variability of critical PAR values in the case of shallower DCM (Fig. 7b).

Figure 7 (a) DCM depth (z_DCM, x axis) compared to the euphotic depth (z_eu, y axis) for both modelled (red dot) and measured results (black dot). Red box (top left) reports statistics for model z_DCM versus z_eu, whereas the black box (bottom right) shows statistics for z_DCM derived from Chl data versus z_eu. (b) Irradiance values (y axis) at DCM depth (x axis) for both modelled (red dot) and measured results (black dot). Horizontal blue line marks the 5.8 irradiance threshold (units µmol photons m⁻² s⁻¹) as identified in Mignot et al. (2014).

The Mediterranean Sea is a nutrient-limited basin (e.g. Crispi et al., 2001; Lazzari et al., 2016; Powley et al., 2017); therefore an insight into the role played by nutrients requires further investigation. Phosphate dynamics show an increase in surface Chl driven by nutrient uptake in upper layers due to convective mixing (Figs. 3 to 6). During stratification periods, the phosphocline follows the euphotic layer threshold. Results from the sensitivity analyses (Sect. S4 in the Supplement) illustrate that the role of nutrients is significant in regulating Chl concentration in DCM.

The REF simulation is forced by PAR measurements; hence it is possible to evaluate the direct impact of nutrients' vertical fluxes (Cullen, 2015) compared to light on DCM properties. The effect of self-shading by Chl and CDOM can increase the role of nutrients in terms of DCM depth modulation, which can be evaluated only by using bio-optical models where attenuation is regulated by Chl or CDOM as presented in Sect. 3.3.

3.2 Vertical mixing models

As shown in the previous section, the vertical distribution of Chl displays a distinct variability, which can be at least partially ascribed to mixing. Typically, higher vertical eddy diffusivity values imply smoother structures. During the stratification phase, when DCM forms, the controlling mixing parameter is the background diffusivity $D_{\mathrm{v}}^{\mathrm{background}}$.

Simplified theoretical models, such as KiSS (Kierstead and Slobodkin, 1953; Skellam, 1951), can provide rough quantitative scales in order to determine the minimum vertical length scales (L₀) that allow the formation of stable biomass patches (Ryabov and Blasius, 2008), including the DCM, in a steady-state hypothesis:

$$\begin{array}{}\text{(7)}& L_{\mathrm{0}}\propto \sqrt{{\displaystyle \frac{D_{\mathrm{v}}}{\mathit{\mu }}}},\end{array}$$

where D_v is the vertical diffusivity coefficient and μ is the growth rate (in stratified conditions $D_{\mathrm{v}}=D_{\mathrm{v}}^{\mathrm{background}}$). An increase in background diffusion over a critical value will produce a dispersal of patchy structures (i.e. a relative maximum of Chl concentration), whereas an increase in growth rate μ can drive the formation of finer-scale structures.

The dynamics presented in this study are, however, more complex than KiSS, in terms of both BGC-Argo data and the 1-D BFM model. Vertical eddy diffusivity can simultaneously affect nutrients, phytoplankton, and mesozooplankton with intricate interactions, which in turn make it difficult to derive analytical solutions. Moreover, unlike KiSS, both the model and environment are hardly ever in a steady-state condition, as a result of daily and seasonal oscillations in physical forcings, which are essentially due to variability in diel irradiance and vertical mixing.

Several simulations, labelled MLD1, MLD2, MLD3, and MLD4, were carried out by changing the background vertical eddy diffusivity coefficient values ($D_{\mathrm{v}}^{\mathrm{background}}$) by 2 orders of magnitude, i.e. from 10⁻⁶ to 10⁻⁴ m² s⁻¹ (Table 1). This subset of simulations (with float-derived PAR) clusters at a correlation of approximately 0.8 with RMSD of DCM depth between 10 and 15 m, the same order of what was found by Salon et al. (2019). Perturbing $D_{\mathrm{v}}^{\mathrm{background}}$ over 2 orders of magnitude (from REF to MLD4) shows that the impact on DCM position is lower than 10 m (Fig. 8b), with an uplift of DCM depth with higher $D_{\mathrm{v}}^{\mathrm{background}}$ (Cullen, 2015). The direct impact of eddy diffusivity appears lower compared to direct light modulation on DCM depth.

Figure 8 (a) Taylor diagram showing model skill in reproducing DCM depth compared to data. Correlation is represented by the angle with a positive x axis, whereas distances from the origin depict standard deviations. Green circles illustrate iso-contours of RMSD levels; all units are in metres. (b) Target diagram showing model skill in reproducing DCM depth compared to data. Distance to the origin defines the RMSD; all units are in metres. The position on the x axis is positive if the model standard deviation is higher than the one from data results and negative in the opposite situation. For the sake of completeness, all models considered are reported in these summarizing skill diagrams.

3.3 Bio-optical models

The alternative bio-optical models (OPT1, OPT2, OPT3) were slightly less accurate compared to REF, with correlation decreasing from 0.8 to 0.6–0.5 (Fig. 8a). The OPT3 simulation showed a bias very close to zero, thus suggesting an intermediate skill compared to assimilated PAR simulations (e.g. REF) and the bio-optical models (OPT1 and OPT2). OPT1 and OPT2 cluster of simulations shows slightly lower correlations and an increase in bias (almost zero for OPT1 and from 6 to −14 m for OPT2a to OPT2d) with a RMSD of approximately 20 m in all cases.

Some of the bio-optical models considered, in particular OPT1, OPT2a, and OPT2b, reproduce the DCM depth gradient between western and eastern subbasins with a tolerance of ±10 m (Fig. 9). In previous studies (Crispi et al., 2002; Lazzari et al., 2012), the correct simulation of the DCM depth longitudinal gradient was obtained by forcing the system with a space-time-dependent light attenuation parameter based on Secchi disc climatology or on satellite K_d(490) data. Both empirical approaches prevent us from understanding whether the origin of such gradients is directly related to external forcings or if it can be interpreted as a self-emerging property (i.e. related to the appearance of features which are not directly and explicitly imposed from the choice of boundary conditions or model parameters used in the numerical experiment; de Mora et al., 2016). Results suggest that a gradient in DCM depth could be partially reproduced and explained in terms of internal biogeochemical processes and partially due to external forcings (i.e. downward irradiance and nutrient initial conditions), even without considering lateral dynamics (Fig. 9a, b).

Figure 9 (a) Monthly average of DCM depth for west (blue) and east (orange) profiles derived from REF simulation. Circles and crosses are the mean; triangles are medians. (b) Scatter plots of the residual difference between measured and modelled DCM. The x axis reports model configurations listed in Table 1. On the y axis, residuals' median values for west (blue) and east (orange) profiles are shown. Triangles indicate the 25th and 75th percentiles. (c) K_d for west and east subbasins during the stratified period; diamonds indicate the median over the vertical column; 25th and 75th percentiles are the horizontal lines. Crosses show the maximum over the vertical column. Panel (d) is the same as (c) but with double initial nutrient concentrations for the western basin simulations.

A direct analysis of the impact of alternative bio-optical models on light attenuation (Fig. 9c) indicates that the simulated eastern basin waters present generally lower K_d values (and lower dispersion around the median) for REF and OPT3. In other cases, where self shading is included, the variability is driven by Chl (from OPT1 to OPT4c) or by Chl and CDOM (OPT5), as bio-optical model parameters do not depend on space and time explicitly. West–east gradients are higher for maximum light attenuation along the water column (cross mark, Fig. 9c) where Chl concentration is higher. For OPT3, the average and maximum K_d overlap since K_d is parametrized as constant along the water column.

The average surface PAR of the data set considered is higher in eastern areas, especially during the months of January (40 %), September (15 %), October (22 %), November (36 %), and December (16 %), probably due to clearer atmospheric conditions (image not shown). During summer, when DCM stabilizes, the west–east differences in measured surface PAR are lower and oscillate around 10 %; however they still contribute to increasing irradiance penetration at deeper layers.

The western and eastern subbasins are also different in terms of nutrient regimes that in turn impact biogeochemical dynamics and the DCM depth gradient in non-trivial ways. The role of nutrients can be evaluated by perturbing initial conditions for the trajectories starting in the western subbasin (see Sect. S4 in the Supplement). Results indicate that increased nutrients in the western subbasin cause an amplification of the west–east light attenuation gradients (Fig. 9d) due to the increase in Chl.

The emerging conceptual scheme is that the first-order controlling mechanism for DCM depth is related to light propagation along the water column, as shown in REF and OPT3 simulations. Other tests indicate that nutrients modulate K_d consistently with gradients simulated in REF. The decadal temporal scale of subsurface nutrient variability (Crispi et al., 2001) controlling self-shading mechanisms is longer than that of simulations, suggesting that the role of nutrients in DCM positioning is especially regulated through initial conditions chosen for the present simulations.

Another key factor pertains to shorter wavelengths (400–450 nm) in the visible part of the spectrum: when light penetrates deeper along the water column, compounds like CDOM are more effective in absorbing light and might in turn enhance spatial gradients in irradiance regimes, which could synergistically contribute to a deeper DCM in eastern subbasins. However, with a current mono-spectral formulation, such aspects still cannot be addressed. Multispectral configurations linked with specific PFT and CDOM absorption terms are thus needed for future in-depth studies of the questions raised in the present work (Dutkiewicz et al., 2015).

3.4 Daily variable versus constant PAR forcings

The use of daily averaged irradiance (i.e. with continuous light, CL1, and CL2) was compared against REF that includes the diurnal variability. A consistent reduction of surface Chl concentrations was observed in the former case (Fig. 10), with a correlation lower than REF, affecting (in relative terms) the values around DCM much less (Fig. 11).

Figure 10 Scatter plot comparing 0–25 m average surface Chl versus BGC-Argo float data for the stratified period condition (DCM > 40 m).

Figure 11 Example of a 6-week time series of vertical profiles from the lovbio035b BGC-Argo float (Fig. 4, from week 28 to week 33) based on diel variability and constant daily light descriptions, compared to BGC-Argo float Chl values (thicker blue line). The horizontal dashed blue line represents the euphotic depth z_eu, whereas the dashed black line indicates the depth where measured PAR equals 5.8 µmol photons m⁻² s⁻¹ as identified in Mignot et al. (2014). The legend reports model configurations listed in Table 1.

Near the surface, phytoplankton is limited by low nutrients (especially in eastern subbasins), whereas closer to DCM the trophic limitation is weaker, sometimes nonexistent (Behrenfeld and Boss, 2003; Behrenfeld et al., 2004). One possible explanation could be that light limitation at the DCM at low irradiance values is almost linear; thus the PAR daily averaging effects have a larger impact at the surface, where light limitation is highly non-linear due to saturation. Furthermore, the BFM formulation for Chl acclimation (Geider et al., 1998) in the case of diurnal variability generates an increase in Chl-to-carbon (Chl:C) ratio. This could in turn have important consequences in operational applications, where data assimilation is employed for model skill improvement: at the surface, the adoption of a diurnal cycle formulation could reduce corrections made by the assimilation scheme and therefore minimize possible spurious trends introduced by it (Gehlen et al., 2015).

Combining daily-averaged irradiances with the lowest diffusivity rates ($D_{\mathrm{v}}^{\mathrm{background}}={\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹, simulation CL2) results in additional relative Chl maxima at surface layers (Fig. 11, panel “T=33 weeks”), as well as in increased patchiness of the whole vertical profile. Similar Chl profiles with multiple subsurface maxima were identified in a comprehensive fluorescence data analysis in the Mediterranean Sea (Lavigne et al., 2015). Theoretical considerations predict different maxima along the water column based on the Tilman resource competition theory applied to a heterogeneous system (Ryabov and Blasius, 2011). At this stage, however, it is difficult to assess whether the patchy structures observed in data and model are, for various reasons, realistic or artefactual. Nonetheless, it can be ascertained that the background diffusion needed to maintain such structures in model simulations is very low. As a result, within the framework of currently used mathematical formulations in the 1-D BFM model, the inclusion of diurnal variability tends to reduce the formation of fine-scaled structures that could be interpreted in terms of a reduction in diel growth (μ) or seen as a possible perturbation that has an equivalent effect of an increased diffusion.

3.5 Bio-optical models with CDOM formulation

OPT4 and OPT5 simulations take into consideration CDOM dynamics by including an additional term in OPT2a, where light attenuation by PAR was described only in terms of Chl. In OPT4a, b, and c, CDOM is parametrized as dead Chl by changing only the rate of Chl decay from 1 d to 1 month. Such a simplified dynamics description derives from the high correlation observed between Chl and CDOM in Morel and Maritorena (2001). However, no analysis was carried out within the present data set to corroborate findings from Morel and Maritorena (2001) due to a lack of information on CDOM fluorescence. In all three model configurations, the dead Chl accumulation results in higher turbidity levels that in turn reduce light penetration depths. This is quantified by significantly negative DCM biases (over 40 m in OPT04c), which result in shallower DCM compared to BGC-Argo-derived profiles since the attenuation of Chl is overestimated even when considering the fastest degradation rates (Fig. 8). The experiment OPT5 mimics the CDOM dynamics described in Dutkiewicz et al. (2015) where a lower bias is observed compared to the (over)simplified OPT4 tests (where correlation coefficients range from 0.6 to less than 0.1 for OPT4a to OPT4c respectively). OPT5 still results in a negative bias of around 10 m compared to the values from −25 to −40 m for OPT4a to OPT4c. The model, regardless of initial conditions, correctly drives CDOM absorption coefficients in deeper layers to low values, while an enhanced surface production reinforces mineralization and bleaching (Fig. 12). Results of CDOM variability from the BOUSSOLE site (northwest Mediterranean; Antoine et al., 2008) show that CDOM absorption ranges to a maximum value of 0.07 m⁻¹ and indicate that there is a temporal delay between phytoplankton bloom and a maximum in CDOM absorption (Fig. 3 in Organelli et al., 2014), whereas deeper layers (below 100 m) have generally lower CDOM absorption. The data set shown in Organelli et al. (2014) evidences that cycles of CDOM accumulation are followed by depletion in the upper 10 m due to photodegradation in summer. In the model results presented here, bleaching has a deeper effect over the entire CDOM “productive” layer (see red and blue lines, Fig. 12), while the subsurface CDOM maximum is not reproduced. The lack of CDOM accumulation in deeper layers for the OPT5 configuration hinders a proper analysis of mechanisms related to the emergence of CDOM from subsurface dark layers. Improving model dynamics calibrations could possibly be achieved by utilizing information on CDOM light absorption from BGC-Argo float measurements (Xing et al., 2012; Organelli et al., 2017b).

Figure 12 Hovmöller diagrams for BGC-Argo float lovbio068d (WMO code 6901648) deployed in the northwest Mediterranean showing PAR (a), total Chl (b), and CDOM (c) simulated by model configuration OPT5. The white, red, and blue lines depict the euphotic, 100 %, and 10 % bleaching depths respectively.

4 Conclusions

The coupled modelling–experimental approach presented here provides a robust and accurate reproduction of the DCM depth variability across the Mediterranean Sea. Such a combined configuration of this kind can integrate multi-data measurements provided by BGC-Argo floats in a single framework. DCM is a ubiquitous feature of the Chl vertical structure in the Mediterranean, and different forcing conditions generate geographical gradients in DCM characteristics (i.e. shallower DCM in western regions, deepening eastwards). Second-order features, such as impulsive vertical spikes or specific patterns observed in BGC-Argo profiles, are also qualitatively reproduced. Results for the reference simulation, where measured PAR is adopted, are summarized as follows:

• mixing and irradiance propagation control Chl dynamics,

• DCM position is mostly controlled by PAR,

• nutrients control the amount of biomass at DCM.

It was demonstrated that vertical processes considered in the 1-D model, such as irradiance regimes and vertical mixing, allow us to properly reconstruct a large part of Chl dynamics, which was also quantified by skill diagrams. Moreover, the role of nutrients in modulating self-shading (as inferred with bio-optical alternative experiments) appears relevant to shape west–east heterogeneity of vertical light attenuation.

The emerging conceptual scheme is that DCM gradients are directly controlled by irradiance modulation, in turn controlled through bio-optical processes which change attenuation according to optically active substances (e.g. Chl, CDOM). Nutrients can impact attenuation by regulating Chl concentrations. The timescale of the subsurface nutrient inventory variability is longer than the ones considered in the present simulations; therefore initial conditions have an impact on west–east gradients.

Such data-rich experiments, combined with a 1-D numerical model, could also be considered a useful tool for a broader community, rather than only for biogeochemical modellers, in particular to address process studies.

The presented approach might also be strategical to quantify the amount of measured signal related to vertical dynamics and the one derived from other processes, like horizontal advection and subduction of water masses. The usage of PAR measured from BGC-Argo floats (used in REF, CL1, CL2, MLD1, MLD2, MLD3, and MLD4) provides higher correlations compared to configurations with alternative bio-optical models (used in OPT1, OPT2, OPT3, OPT4, and OPT5). The comparison of different bio-optical models indicates that, when lacking direct measurements of PAR in subsurface layers, the most fitting alternatives would be OPT3, OPT2a, and OPT1, resulting in lower bias and higher correlation coefficients (between 0.5 and 0.7), as well as lower RMSD values compared to REF. Our analysis can also help determine how the use of light fully integrated in the visible range of the spectrum (400 to 700 nm, REF) improves predictions when compared to simplified approaches (i.e. all the OPT simulations here considered).

Moreover, we show that on the timescales here considered (months), vertical processes are more relevant than horizontal ones, and the parametrizations used in the biogeochemical model are adequate to describe the main processes taking place on these scales.

Results also highlight the strategic relevance of BGC-Argo data: temperature, salinity, and radiometric parameters encapsulate fundamental information for the reconstruction of primary producer dynamics and are paramount to investigate hypotheses concerning DCM formation. CDOM fluorescence data measured by BGC-Argo floats could be integrated in simulations to further infer and reconstruct the observed biogeochemical processes.

Considering a general 3-D biogeochemical model, it is not possible to have a full data coverage of the in-water PAR field without a fully coupled radiative transfer model. The proposed approach could thus be exported and generalized at a global scale.