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

The present work is based on a dataset comprised of 31 Biogeochemical (BGC) Argo floats that collected 0-1000 m vertical profiles of biogeochemical and optical data from 2012 to 2016 in the Mediterranean Sea. The dataset was integrated in 1-dimensional model simulations following the trajectories of each float and considering measured photosynthetically avail5 able radiation (PAR) profiles as the reference light parameterization. The simulations were aimed to be consistent with data measured by float sensors, especially in terms of the deep chlorophyll maximum (DCM) depth. Moreover, we tested several light models in order to estimate their impact on modeled biogeochemical properties, including self-shading dynamics based on chlorophyll and colored dissolved organic matter (CDOM) concentrations. The results, evaluated with the corresponding in-situ BGC-Argo chlorophyll data, indicate that the proposed approach allows to properly simulate the chlorophyll dynamics 10 and illustrate how PAR and vertical mixing are essential environmental regulation factors driving primary producers dynamics. The higher skills are reached using in-situ PAR, but some of the alternative bio-optical models here presented show comparable skill in reproducing DCM depth spatial variability. Simulation results show that during the stratification phase the diel cycle has significant impact on the surface chlorophyll regimes. The approach here presented serves as a computationally smooth solution to analyse BGC-Argo floats data and to corroborate hypotheses on their spatio-temporal variability. 15

and 2016 into a one-dimensional model to simulate the vertical and temporal variability of algal chlorophyll concentrations. In addition to PAR as input, alternative light and vertical mixing models were considered. We evaluate the models' skill to reproduce the spatial and temporal variability of deep 15 chlorophyll maxima as observed by BGC-Argo floats. The assumptions used to set up the one-dimensional model are validated by the high number of co-located in situ measurements. Our results illustrate the key role of PAR and vertical mixing in shaping the vertical dynamics of primary produc-20 ers in the Mediterranean Sea. Moreover, we demonstrate the importance of modelling the diel cycle to simulate chlorophyll concentrations in stratified waters at the surface.CE2

Introduction
In most biogeochemical models the description of optics is 25 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 30 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 35 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 dras-40 tic 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(Organelli et al., , 2017a) have made 45 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 50 for understanding marine biogeochemistry processes. At present, no studies have 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 5 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. 10 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 15 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 20 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 25 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 30 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. Gen-35 eral remarks are illustrated in the Conclusions section.
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 55 (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 60 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 of the Supplement).
Vertical profiles of Chl concentration were quality-65 controlled according to the procedure of the international BGC-Argo program 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 se-70 ries Chl fluorometers, Chl concentrations were corrected by a factor of 0.5 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). 75 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 80 (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 85 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 of the Supplement.

1-D biogeochemical model
Biogeochemical processes have been simulated according to the voxel approach ("volume element with biological con-5 tent 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 10 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 pro-15 cesses 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 20 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 25 model results and experimental data. The model is formulated through a system of partial differential equations: where C i is the ith biogeochemical tracer simulated (i = 1, 30 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 35 programmed at around local noon, the variability related to diurnal variation in solar irradiance is taken into consideration according to Kirk (1994). The biogeochemical model BFM (Vichi et al., 2013) is a biomass-based numerical model that simulates the biogeo-40 chemical 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;Lazzari et al., 45 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, phos-50 phorus, and nitrogen concentrations. Particulate and dis-solved 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 fu-55 ture analysis (according to data availability and optical model complexity) a study of plankton functional type (PFT) resource competition dynamics Blasius, 2011, 2014).

70
-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 75 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, 80 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, 85 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 90 impact of using a constant light approximation (CL1 and CL2 configurations) rather than following the diurnal light variation on Chl distribution.
Initial conditions for all biogeochemical variables of BFM are provided by the CMEMS reanalysis of Mediterranean Sea biogeochemistry (period 1999Teruzzi et al., 100 2014) produced by the MedBFM model system. The initialization profiles for our 1-D configuration are extracted from as OPT2a + Chl degradation to CDOM − timescale 1 d OPT4b as OPT2a + Chl degradation to CDOM − timescale 1 week OPT4c as OPT2a + Chl degradation to CDOM − timescale 1 month OPT5 as OPT2a + CDOM following Dutkiewicz et al. (2015) 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. 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. Ap-10 proaches and impacts of using different parametrizations to reconstruct mixing along the water column are shown and discussed in Sect. 2.2.1.

Vertical mixing models
Vertical mixing is estimated from potential density (obtained 15 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 (2) 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): In simulations MLD1, MLD2, MLD3, and MLD4, D background v values were perturbed by 2 orders of magnitude (from 10 −6 to 10 −4 m 2 s −1 ) in order to estimate the impact such variations have on modelled Chl profile shapes com-30 pared to measured ones (see Table 1).

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 35 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): In OPT2 models, statistical regressions were carried out between K d (PAR) and Chl measured by BGC-Argo floats at 5 four different depth ranges: 150, 75, 45 and 30 m (OPT2a to OPT2d; see Table 2 for details): 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 15 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 mod-20 els 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 25 (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 expo-30 nential 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), 35 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, 40 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 45 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 50 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.

Statistical analysis
According to the work's objectives, four classes of simu-55 lations 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 60 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 65 than 0.1 mg m −3 . 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 70 parameter, which was calculated following Eq. (6): 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 75 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 Supplement), which considered DCM width, DCM magnitude, and surface Chl. 80 Results indicate that DCM depth is the most effective feature the model is able to reproduce.

Reference simulation
The assimilation of PAR profiles helped to accurately esti-85 mate 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 90 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).

95
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 100 mixing occurring predominantly in autumn and early winter. At the ocean surface, the increase in Chl is triggered by  (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 repre-25 sentation 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 sophis-30 ticated 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 signifi-35 cant 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 −2 s −1 isolume (Fig. 7b, blue line). Data 40 and model outputs in the present study show a higher variability of critical PAR values in the case of shallower DCM (Fig. 7b).
The Mediterranean Sea is a nutrient-limited basin (e.g. Crispi et al., 2001;Lazzari et al., 2016;Powley et al., 2017); 45 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. 50 Results from the sensitivity analyses (Sect. S4 in 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' 55 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 60 presented in Sect. 3.3.

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 dif-65 fusivity values imply smoother structures. During the stratification phase, when DCM forms, the controlling mixing parameter is the background diffusivity D background v . Simplified theoretical models, such as KiSS (Kierstead and Slobodkin, 1953;Skellam, 1951), can provide rough 70 quantitative scales in order to determine the minimum vertical length scales (L 0 ) that allow the formation of stable biomass patches (Ryabov and Blasius, 2008), including the DCM, in a steady-state hypothesis: where D v is the vertical diffusivity coefficient and µ is the growth rate (in stratified conditions D v = D background v ). 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 5 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 15 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 background v ) by 2 orders of magnitude, i.e. from 10 −6 to 10 −4 m 2 s −1 (Table 1). This 20 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 background v over 2 orders of magnitude (from REF to MLD4) shows that the impact on 25 DCM position is lower than 10 m (Fig. 8b), with an uplift of DCM depth with higher D background v (Cullen, 2015). The direct impact of eddy diffusivity appears lower compared to direct light modulation on DCM depth. 30 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 35 (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.

5
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 10 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 15 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 gra-20 dient 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).

25
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 30 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. 35 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 40 (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.

5
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 modu- 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 25 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 ab-30 sorption terms are thus needed for future in-depth studies of the questions raised in the present work (Dutkiewicz et al., 2015).

Daily variable versus constant PAR forcings
The use of daily averaged irradiance (i.e. with continuous 35 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). 40 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 5 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 applica-10 tions, 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).

15
Combining daily-averaged irradiances with the lowest diffusivity rates (D background v = 10 −6 m 2 s −1 , 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 20 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 heteroge-25 neous 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 30 model simulations is very low. As a result, within the frame-work 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 35 as a possible perturbation that has an equivalent effect of an increased diffusion.

Bio-optical models with CDOM formulation
OPT4 and OPT5 simulations take into consideration CDOM dynamics by including an additional term in OPT2a, where 40 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 45 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 turbid-50 ity 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 55 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 60 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 BOUS-5 SOLE site (northwest Mediterranean; Antoine et al., 2008) show that CDOM absorption ranges to a maximum value of 0.07 m −1 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 10 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 15 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 analy-  Table 1. sis 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).

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 can integrate in a single framework multi-data measurements provided by BGC-Argo floats. 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).

15
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:

20
mixing and irradiance propagation control Chl dynamics, It was demonstrated that vertical processes considered in the 25 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 30 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, 35 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.

40
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 modelers, in particular to address process studies.
The presented approach might also be strategical to quan-45 tify the amount of measured signal related to vertical dynamics and the one derived from other processes, such as hori-zontal 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 biooptical models (used in OPT1, OPT2, OPT3, OPT4, and 5 OPT5). CL1 (without diurnal cycle) shows the overall highest correlation, comparable with REF. 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 10 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 sim-15 plified approaches (i.e. all the OPT simulations here considered).
These results also highlight the strategic relevance of BGC-Argo data: temperature, salinity, and radiometric parameters encapsulate fundamental information for the recon-20 struction 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.CE3 25 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. Such an approach could thus be exported and generalized at a global scale. 30 Code and data availability. The BFM biogeochemical model and its documentation can be downloaded at the following address: http://bfm-community.eu/ (last access: 22 April 2019). The qualitycontrolled databases used in the present paper are publicly available from the SEANOE (SEA scieNtific Open data Edition) pub- 35 lisher at https://doi.org/10.17882/49388  and https://doi.org/10.17882/47142 (Organelli et al., 2016b) for vertical profiles and products within the first optical depth respectively.

40
Author contributions. ET and PL have designed the paper. ET performed the BGC-Argo data analysis, and PL performed the simulations. ET, PL, EO, SS, CS, FD, and PC contributed to the paper writing.
Competing interests. The authors declare that they have no conflict 45 of interest.
Acknowledgements. This study has been conducted using EU Copernicus Marine Service information. The simulations were performed in the framework of the ISCRA C project NOVBIOGE (HP10C8C9O6), created by CINECA, Italy. We acknowledge spon-50 sorship from the MISTRALS-MERMEX project. TS1 Financial support. This work is part of the PhD project of Elena Terzić (funded under the CMEMS contract for the Biogeochemistry Production Unit for the Mediterranean Sea) and of the BIOP-TIMOD CMEMS Service Evolution project. CMEMS is imple-55 mented by Mercator Ocean International in the framework of a delegation agreement with the European Union. This work was supported by the French "Equipement d'avenir" NAOS project (Novel Argo Ocean Observing System) funded by Agence Nationale de la Recherche (grant no. ANR J11R107-F); the "Remotely-sensed 60 biogeochemical cycles of the oceans -remOcean" project funded by the European Research Council (grant no. 246777); the Argo-Italy project funded by the Italian Ministry of Education, University and Research; and the French Bio-Argo program -Bio-Argo France funded by CNES-TOSCA, LEFE Cyber, and GMMC.

65
Review statement. This paper was edited by Christine Klaas and reviewed by Zarko Kovac and Maurizio Ribera d'Alcala.
Bissett, W., Walsh, J., Dieterle, D., and Carder, K.: Carbon cycling in the upper waters of the Sargasso Sea: I. Numerical simulation of differential carbon and nitrogen fluxes, Deep-Sea Res. Pt. I, 46, 205-269, 1999. Bricaud, A., Morel, A., and Prieur, L.: Absorption by dissolved organic matter of the sea (yellow substance) in the UV and visible domains 1, Limnol. Oceanogr., 26, 43-53, 1981. Crispi, G., Mosetti, R., Solidoro, C., and Crise, A.: Nutrients cycling in Mediterranean basins: the role of the biological pump in the trophic regime, Ecol. Model., 138, 101-114, 2001. Remarks from the language copy-editor CE1 In future rounds of proofreading, if possible please provide your comments and corrections using the highlight-andcomment function in Adobe Reader as this reduces the risk of missing or incorrectly inserting a comment and significantly reduces the time required for proofreading.

CE2
Please note that such extensive changes are no longer appropriate at this stage of the publication process and must be approved by the editor. Please also note that unnecessary stylistic changes have also not been inserted, as these are also no longer appropriate. If the corrections need to be inserted, please provide a short statement regarding these corrections that can be forwarded by us to the editor.

CE3
Please see previous remark regarding these corrections.
Remarks from the typesetter TS1 Please confirm placement.