Radiative Transfer Modeling With Biogeochemical-Argo Float Data in the Mediterranean Sea.

A radiative transfer model was parameterized and validated using Biogeochemical Argo float data acquired between 2012 and 2017 across the Mediterranean Sea. Fluorescence-derived chlorophyll a concentration, particulate optical backscattering at 700 nm, and fluorescence of chromophoric dissolved organic matter (CDOM) were used to parametrize the light absorption and scattering coefficients of the optically significant water constituents (such as pure water, non-algal particles, CDOM, and phytoplankton). The model was validated with in situ downwelling irradiance profiles and apparent optical properties derived both from irradiance profiles and satellite data, such as the diffuse attenuation coefficients and remote sensing reflectance. Results showed that by using regional parameterizations that are not only related to chlorophyll concentration and vertical distribution, the model was able to capture a more accurate spectral response in the examined wavelength range compared to chlorophyll-related (or Case 1) optical models. When using alternative models that incorporated also measurements of CDOM fluorescence or particulate optical backscattering, the model skill increased at all examined wavelengths. Finally, using a multi-spectral optical configuration also enabled the estimation of the relative contribution of separate water constituents in the examined spectral range. Simulations including non-algal particles and CDOM performed up to 61% and 79% better than when considering the optical properties of pure seawater alone. Moreover, a simulation including phytoplankton light absorption resulted in an error reduction of up to 42%, especially at 412 nm and with a more uniform response at the wavelengths considered.

Introduction

Physical processes in the ocean are fundamentally related to the biogeochemical and light environments and thus impact the growth of phytoplankton, which requires both light and nutrients. In order to get a complete grasp of the complex interaction of physics and biogeochemistry, it is essential to include also numerical models with a physically accurate representation of optics. This in turn enables to quantify with greater precision the impacts of changes in biogeochemical processes (such as primary productivity), circulation patterns (through heat transfer), or the sole nature of absorbing and scattering matter in the examined water body.

Despite this knowledge on the importance of optics, most biogeochemical models, notwithstanding their increasing spatial resolutions, shorter computational times, and improved complexity, still employ an oversimplified methodology for optical calculations, usually to predict photosynthetically available radiation without considering its spectral dependency. One of the few attempts to overcome these limitations have been achieved in Fujii et al. (2007), Dutkiewicz et al. (2015), Mobley et al. (2015), and Gregg and Rousseaux (2016, 2017). In order to improve prediction capabilities on marine biogeochemical features, the implementation of multi‐ or hyperspectral optical modeling solutions remains therefore essential, also in order to follow up with the pace of such approaches already adopted in remote sensing and in situ platforms used to successfully observe aquatic biogeochemical phenomena.

With the development of satellites and in situ autonomous platforms, models tend to integrate measurements, either through data assimilation to correct for the numerical drift, or for the validation of the model itself. Combining numerical approaches with experimental data, both from in situ sensors, as well as from satellites, compensates for the shortcomings of each of the separate methods. While model simulations have no spatio‐temporal limitations, they are after all numerical representations of natural phenomena that need to be validated accordingly. Satellite remote sensing in the visible range of the spectrum discards more than 90% of the total signal as atmospheric noise, and it is limited to a certain spatial resolution and temporal frequency, reaching only surface layers under clear skies. The emergence of autonomous platforms, such as floats and gliders, can, on the other hand, provide more information also at greater depths, regardless of the sky conditions, and thus surpass the greatest limitations of satellites. Furthermore, they can provide measurements of additional variables that at present cannot be retrieved by satellites. Such sensors are, however, still spatio‐temporally heterogeneously distributed and cannot fully replace the synoptic coverage that satellites have.

In order to resolve the distribution of light in the water, information on absorption and scattering properties of the optically significant matter is needed. As such measurements are difficult to obtain, unless field measurements are carried out, semi‐empirical relationships between biogeochemical quantities and inherent optical properties (i.e., absorption and scattering spectra) are widely used to facilitate calculations of in‐water light propagation. Measurements of biogeochemical quantities, especially phytoplankton chlorophyll a, are most widely available, hence the advantage in using them for parameterization purposes.

The Mediterranean Sea has been defined as a bio‐optically anomalous region (Bricaud et al., 2002; Corsini et al., 2002; D’Ortenzio et al., 2002; Gitelson et al., 1996; Lee & Hu, 2006; Loisel et al., 2011; Morel, Claustre, et al., 2007; Morel & Gentili, 2009; Organelli, Claustre, et al., 2017; Volpe et al., 2007), so that such global empirical algorithms, both for satellite remote sensing as well as in situ applications, are less accurate. With the adoption of global relationships, Mediterranean waters appeared greener for a given chlorophyll a level than waters in other regions (Claustre et al., 2002; Morel & Gentili, 2009).

Among the possible causes for such phenomena could be:

Specific phytoplankton community structure (cell size, pigment packaging, pigment composition, and photophysiology), which can affect phytoplankton absorption () and particle backscattering signals ().

Excess of non‐algal (mineral) particles (NAP), such as Saharan dust (influencing and ).

Excess of chromophoric dissolved organic matter (CDOM), influencing ;

either separately, or as a combination of several factors (all absorption and scattering coefficients are in units of ).

In the past decade, the development of new technologies for the acquisition and analysis of bio‐optical variables has brought new insights on CDOM dynamics, size and composition of algal communities, absorption by phytoplankton () and non‐algal particles (), as well as on particulate backscattering (). Since 2012, a large array of autonomous Biogeochemical Argo (BGC‐Argo) floats has been deployed, measuring a whole set of bio‐optical and biogeochemical variables (IOCCG, 2011), which could fill the gap between sample acquisitions and remote sensing measurements. With their high vertical resolution profiles, BGC‐Argo floats can serve as an additional tool for tackling the bio‐optically anomalous nature of the Mediterranean Sea, also due to their high horizontal and vertical spatial coverage (Organelli, Claustre, et al., 2017).

For this reason, an analysis was hereby carried out to show the possibility of using a large array of BGC‐Argo float measurements both for a radiative transfer model set‐up, as well as for validation purposes. More specifically, profiles of bio‐optical and biogeochemical parameters (i.e., fluorescence‐derived chlorophyll a concentration (Chl), particulate backscattering at 700 nm () and CDOM fluorescence (fDOM)) were used for inherent optical properties (IOP) parameterizations, testing several regionally adopted algorithms for , , and , and particle scattering . Radiometric measurements were, on the other hand, used for model validation.

The aim of this work is to:

Test different IOP configurations, especially Chl‐versus non‐Chl‐related bio‐optical models.

Verify the model performance comparing computation results with BGC‐Argo radiometric profiles and apparent optical properties (i.e., diffuse attenuation coefficients for downwelling irradiance and remote sensing reflectances), the latter derived from radiometric profiles and obtained from satellite remote sensing.

Check the model sensitivity to different absorption spectra within each group of optically significant constituents (pure water, NAP, CDOM, and phytoplankton)

Methods

BGC‐Argo Data Set

The BGC‐Argo data set used in this work was obtained from 39 floats operating between 2012 and 2017. The total number of profiles containing Chl measurements was 5,092, however for the sake of the analysis completeness, a few requirements needed to be met. First, only profiles containing the whole suite of following variables were considered: temperature (T, ), salinity (S, PSU), chlorophyll (Chl, ), particle backscattering at 700 mm (, ), fluorescent component of chromophoric dissolved organic matter (fDOM, ppb of quinine sulfate), and downwelling plane irradiance at 380, 412, and 490 nm (, ). Profiles lacking at least one of the required variables were excluded from further calculations (2,112). Then, only profiles acquired between 10:00 am and 2:00 pm local time were considered in order to obtain radiometric measurements at low solar zenith angles, thus removing additional 396 profiles. The total number of profiles left for the analysis was 2,584. The quality control (QC) procedure of radiometric data follows the steps described in Organelli, Claustre, et al. (2016), whereas the Chl and QC protocols are found in Schmechtig et al. (2014, 2018). Profiles of all variables were uniformly interpolated on a 1 m grid, starting at 0.5 m. As light is one of the key mechanisms controlling the deep chlorophyll maximum (DCM) depth (Cullen, 2015; Mignot et al., 2014), the choice of the maximum depth range corresponds to the maximum DCM values in the Eastern Mediterranean, which is up to 120 m (Christaki et al., 2001). Therefore, additional 370 profiles were discarded that had depths shallower than 150 m. For a successful calculation of the depth derivative of radiometric profiles the diffuse attenuation coefficients of downwelling irradiance (, ) with a non‐linear fit of an exponential function with the least squares method, further conditions needed to be met: the first depth measurement of had to be shallower than 1 m (thus discarding 130 profiles) and the number of measurements within the first 10 m had to be at least 5, which discarded another 757 profiles. Moreover, a condition of less than 30% difference between modeled and synthetic (obtained from the at the first optical depth) values was thus added as the existent quality‐control procedure for radiometric quantities still retains noisy behavior, which resulted in 147 profiles less. After having applied all the QC procedures, the final number of useful profiles for this work resulted in 1,126, spatially distributed as in Figure 1. In order to remove spikes and negative values, all variables except T and S were further corrected by applying a 5‐point median filter, followed by a 7‐point running mean. Negative values were assigned to zero. The fact that the applied smoothing procedure might remove some spikes which could be actually indicators of larger aggregates (Briggs et al., 2011), goes beyond the scope of the present study.

Figure 1

Spatial distribution of the complete data set, after the additional criteria applied to remove profiles which do not meet computation requirements.

In‐Water Radiative Transfer Model

The irradiance distribution along the water column was parameterized into three streams as described in Dutkiewicz et al. (2015) and Gregg and Rousseaux (2016): the direct () and diffuse () downwelling irradiance components and the upwelling diffuse irradiance (). The downwelling plane irradiance is equivalent to the sum of the two downward streams (). The light spectrum was discretized into 25 nm bands covering the range between 350 and 700 nm. For each band, , and were solved as a system of three differential equations: where is the depth, are the attenuation factors, and and the backward and forward scattering factors, respectively. The attenuation factors were calculated as the sum of absorption (), scattering (), and backscattering () coefficients normalized over cosines:

In the three‐stream approach, the shape factors were considered constant , , as well as the average cosines and (Aas, 1987), while where denotes the solar zenith angle corrected with water refraction index. Absorption, scattering and backscattering coefficients were defined as a linear combination of separate water constituents:

Different IOP models to determine , , and are further presented in Section 2.3.

Scattering factors were similarly defined: where is the ratio of backscattering () to total scattering ().

Solving the in‐water radiative transfer model requires three boundary conditions, one for each stream. and were derived from the multi‐spectral atmospheric radiative transfer model OASIM (Gregg & Casey, 2009), specifically validated with the surface irradiance values from the same BGC‐Argo data set in the Mediterranean Sea (Lazzari et al., 2021). The boundary conditions for the upward component were set as . The equations were discretized along depth using the same resolution of the BGC‐Argo data and integrated numerically following the methodology described in Dutkiewicz et al. (2015).

IOP Models

The in‐water radiative transfer analysis comprised of six bigger clusters of IOP simulations, as reported below. The aim of these tests was to show that the model does accurately take into consideration the spectral response based on the selection of appropriate IOPs (both absorption and scattering) and thus correctly resolves the radiative transfer equations. By considering one IOP at the time, it was possible to quantify how much does the IOP for a given optically significant constituent contribute to the relative improvement compared to the base. Furthermore, within each constituent, it was possible to assess the impact that each of the model range has at the output at separate wavelengths. In the following subsections, separate groups of IOP models are described in more detail, along with the upgrades that were tested.

Pure water absorption and scattering (, ).

1. + NAP absorption ().

1. + CDOM absorption ().

1. + phytoplankton absorption ().

1. + particle scattering ().

1. +  +  +  + .

Most models that link biogeochemical quantities with IOPs are assessed for Case 1 water optical types that can be defined as water bodies for which the inherent optical properties (of CDOM and NAP) co‐vary with phytoplankton and hence with Chl concentration (Morel & Prieur, 1977). Even though such empirical relationships can be quite useful for parameter estimations, there exists the tendency to oversimplify the optical response of a generally complex biogeochemical environment, as thoroughly discussed in Mobley et al. (2004). Hence, one of the goals of this paper was to try to compare Case 1 water types with alternative parameterizations that considered additional biogeochemical variables and are described in the following subsections.

Simulation results were verified in two different ways. First, modeled irradiance profiles were matched‐up with measured profiles at all 3 available wavelengths within the upper 150 m of depth. Second, diffuse attenuation coefficients of downwelling plane irradiances () for the first optical depth (i.e., the depth range for which the light at a specific wavelength attenuates e‐fold) were calculated for both modeled and measured profiles. as an apparent optical property (AOP) does have the advantage of conveying more information on IOPs and to a certain extent remove the impact of the external environment's variability (change in sun location, cloud cover, or surface waves, Mobley et al., 2010). The influence of the external factors is however still present, despite the quality‐control procedure introduced by Organelli, Claustre, et al. (2016), resulting in noisy or oddly‐shaped profiles which were discarded by including additional conditions described in Section 2.1. Moreover, at 490 nm it is possible to make a three‐platform comparison including model, float and remote sensing data, which is further elucidated in Section 2.4.

Pure Water Absorption and Scattering

The original in‐water modeling configuration, described in Gregg and Rousseaux (2016), resolves a pure water absorption spectrum based on data from various sources, as reported therein. However, the UV and blue part of the visible spectrum (from 250 to 550 nm) was improved also with more recent spectrophotometric measurements by Mason et al. (2016), which introduced lower values compared to the findings of Smith and Baker (1981), Morel, Gentili, et al. (2007), Pope and Fry (1997), and Lee et al. (2015). Moreover, pure water absorption was accounted also for the influence on seawater optical properties of T and S according to Sullivan et al. (2006). The original values for pure water scattering from Smith and Baker (1981) and Morel, Claustre, et al. (2007) were further upgraded by calculating values based on the method described by Zhang et al. (2009), thereby accounting for the contribution of T and S. The backscattering‐to‐total scattering ratio for water is kept as 0.5 as in Gregg and Rousseaux (2016), assuming an isotropic scattering regime.

Non‐Algal Particles (NAP) Absorption

The non‐algal particles are defined as a composite of living organic particles, such as bacteria, zooplankton, detrital organic matter, and suspended inorganic particles (Mobley et al., 2010). The absorption spectrum, despite its heterogeneous biogeochemical composition, is described with an exponentially decreasing shape from UV to the red part of the spectrum:

The absorption at the reference wavelength, , can be estimated in two ways: either as a function of Chl (a Case 1 optical water type—see Equation 4 in Bricaud et al., 2010), or by considering the range of values measured in the Mediterranean Sea, that is, between 0.0087 and 0.8 (Babin et al., 2003), with the higher values corresponding to highly turbid waters. The slope varies from 0.0178 and 0.0104 , with a mean value of 0.0129 (Babin et al., 2003). It should be noted, however, that the data collected in the work were from coastal regions, therefore the minimum values could also overestimate the contribution of compared to the open ocean. To better reproduce the vertical distribution of NAP, different profile shapes are considered when estimating model IOPs: Case 1 optical types assume a co‐variability with Chl, and additional tests were performed by considering as a better proxy for non‐algal particle vertical distribution.

The signal is comprised both of organic and inorganic particles, however, the separation of the two fractions is at present still not possible to achieve. As it has already been demonstrated that the contribution of detrital non‐algal particles to the total signal can be very high in Mediterranean waters (Bellacicco et al., 2019), a hypothesis was placed to consider as a better parameter than Chl from the BGC‐Argo set of measurements in terms of NAP depth variability. The summary of models is shown in Table 1.

Table 1

List of All the Models Used in This Study

a_NAP
Name	Model	Equation	Profile shape	S_NAP range	a_NAP(443) range
a_NAP_Case1_Chl a_NAP_Case1_b_bp	Bricaud et al. (2010)	0.013Chl^0.615e^−S_NAP(λ−440)	Chl b_bp(700)	0.0104–0.0178	‐
a_NAP_Babin_Chl a_NAP_Babin_b_bp	Babin et al. (2003)	a_NAP(443)e^−S_NAP(λ−443)	Chl b_bp(700)	0.0104–0.0178	0.0087–0.8
Final	Babin et al. (2003)	a_NAP(443)e^−S_NAP(λ−443)	b_bp(700)	0.0129	0.0087

Note. Either Case 1 from Bricaud et al. (2010) for , Morel & Gentili (2009) for and Morel et al. (2002) for or the non‐Case 1 from Babin et al. (2003) for , Organelli & Claustre (2019) for where and Antoine et al. (2011) for . The profile shape follows either Chl, or fDOM, with different ranges of the model parameters for each IOP: , , , and .

fDOM

CDOM Absorption

Similarly to , the spectral response of is also parameterized with a decreasing exponential function:

can be also estimated as a function of Chl from a regional Case 1 model presented in Morel and Gentili (2009) which is based on spectral coefficients of pure water as measured by Pope and Fry (1997). However, given the substantial modification of the absorption spectra in the UV/blue range when following Mason et al. (2016) compared to originally adopted values from Pope and Fry (1997), a set of simulations was tested by subtracting the former with the updated one, , as shown in Equation 15. With previous values, amounted to a higher water absorption, which would have led to a significant underestimation of .

The remaining parameter to estimate was the slope , which can be taken from Babin et al. (2003) and Organelli et al. (2014), that is, ranging between 0.015 and 0.02 , with a mean value of 0.017 . As with NAP, the Case 1 model for was upgraded by considering the fDOM profile shape instead of a vertical parameterization depending only on Chl. Following Organelli and Claustre (2019), was approximated with , which was in turn calculated from the diffuse attenuation coefficient . The latter was derived from the BGC‐Argo irradiance profiles at 380 nm, both for the mixed layer as for the first optical depth. The former was obtained from a potential density threshold value criterion (de Boyer Montégut et al., 2004), whereas the latter corresponds to the e‐folding depth at the specific wavelength. was then estimated from a non‐linear fit with the least squares minimization of an exponential function for both depth ranges, and can be separated into pure water and biogenic components (Morel & Maritorena, 2001): where:

After having subtracted the pure water contribution as estimated in Morel and Maritorena (2001) (i.e., 0.0151 ), the remaining item, , serves as a proxy for . As discussed in Organelli and Claustre (2019), there are several previous studies in the clearest oligotrophic world oceans that have shown that CDOM dominates the light absorption budget at 380 nm (see references therein). at 380 nm contributes less than 20% to total non‐water absorption in clear oligotrophic waters (Bricaud et al., 2010). In the absence of coincident light absorption data to prove that these conclusions hold true also in the present case, other possible sources that affect light attenuation in the UVs, such as light absorption by mycosporine‐like amino acids and NAP, might be excluded or considered negligible. Instead, an additional set of simulations was performed by changing the relative contribution of by assigning a factor ranging between 0.5 and 1. In this way, different relative contributions of other constituents could be assessed while leaving some uncertainty in the method to use as a proxy for only.

Given the fact that the IOP models used here for pure water absorption are the new measurements of Mason et al. (2016), and a T‐S correction is applied (Sullivan et al., 2006; Zhang et al., 2009), different tests were tried in order to calculate as a function of and rather than adopting a constant value. The entire spectrum is then estimated with the slope range of values as described above, with the depth variability analogous to the fDOM shape. The summary of models is shown in Table 1.

Phytoplankton Absorption

For phytoplankton Chl‐specific absorption spectra, data for seven different algal species of varying size were used, with organisms adequate for surface applications, and several strains suitable for both surface and mixed layer (more details in Organelli, Nuccio, et al., 2017). Absorption spectra were obtained for species cultured at the light regime of 100 photons .

The total phytoplankton absorption is computed as the sum of separate phytoplankton functional types (PFT) spectra as shown in Equation 18:

The relative contribution of each PFT to the total Chl concentration, , followed the regional empirical algorithm introduced by Di Cicco et al. (2017, Table 4). For that purpose, seven algal species were merged into six PFTs: Diatoms, Dinoflagellates, Cryptophytes, Green Algae and Prochlorococcus, Prochlorococcus and Synechococcus, Coccolithophores. The relative contribution of Prochlorococcus was divided into 0.5 for the 2 PFTs containing the same species. Original spectra with a 1 nm frequency were converted to 25 nm bins, corresponding to the model spectral resolution.

The regional algorithm of Di Cicco et al. (2017) was validated with in situ data for first 50 m, with the majority of samples in the Western Mediterranean. Apart from the spatio‐temporal bias inherent to ship‐borne measurements with which the relationship was obtained, it is suitable for Chl values in the range between 0.02 and 5.52 . Therefore, Chl values higher than 5.52 or lower than 0.02 have been limited to Chl = 5.52 and 0.02 , respectively. The lower limit was placed in order to avoid numerical instabilities, whereas the higher limit was reached in only 5 profiles out of 1,126, all of them present in the North‐Western Mediterranean during spring blooms (i.e., 5.71, 5.77, 5.82, 5.96, 5.53 ). No special features were observed in any of the limiting cases.

Particle Scattering

Unlike the model set‐up in Gregg and Rousseaux (2016), the particle scattering is resolved as a total sum, and not partitioned into the relative scattering contributions of separate PFTs plus NAP. Following Equation 14 in Morel et al. (2002), is expressed as a function of Chl: where if and if . values are between −1 and 0. Commonly used in earlier models, the value of is derived from Mie theory and is known to be valid only for non‐absorbing particles with a Junge particle size distribution slope of −4 (Van de Hulst, 1981) with a particle size range between  = 0 and  =  (Boss et al., 2001). Similarly to the PFT regional algorithm modification, is calculated as if Chl were equal to 0.02 for values lower than the minimum concentration. Both Chl and vertical profiles were taken into consideration to account for the depth variability. Alternatively, from BGC‐Argo floats can be also used to estimate . A spectrum of can be obtained from Equation 20: where describes the backscattering spectral slope and can be related to particle size distributions, assuming that the particles are non‐absorbing. Lower slope values (around 0–1) indicate the presence of larger particles and vice‐versa. The range of tested values was between 0 and 4, where the highest slope value agrees with measurements at the BOUSSOLE buoy (Antoine et al., 2011), and a mean value of 2 was found according to Organelli, Bricaud, et al. (2016). The relative contribution of back‐to total particle scattering can be quantified with a known backscattering ratio : in the present set of simulations ranging between 0.002 and 0.015 and spectrally constant (Antoine et al., 2011). The depth variability follows profiles.

Remote Sensing Data

Both AOPs, that is, remote sensing reflectance and , can be described as functions of absorption and (back)scattering coefficients (Gordon, 1989; Morel & Gentili, 1993).

In order to compare model data with satellite measurements, the calculation of in‐water remote sensing reflectance from the model was carried out by following: where the calculation of Q, a function of wavelength , Chl, and solar zenith angle , follows the procedure introduced by Morel et al. (2002):

Values of and are interpolated from the look‐up Table 2 in the Case 1 model from Morel et al. (2002). Surface Chl values were taken from float measurements at the shallowest depth. In case of Chl concentrations below 0.03 , and were taken from the minimum value.

Table 2

Summary of the Skill of Simulations (Only Pure Water IOPs, Pure Water IOPs and , Pure Water IOPs and , and Pure Water IOPs Plus Separate IOPs) Quantified in Terms of Relative Changes in Root Mean Square Error (), Bias (), and Pearson Correlation Coefficient () With Respect to the Previous Model Configuration (i.e., One Line Above in the Table)

Name	Model	Profile shape	Metrics
			380 nm			412 nm			490 nm
			ΔRMSE[%]	Δbias[%]	Δr[%]	ΔRMSE[%]	Δbias[%]	Δr[%]	ΔRMSE[%]	Δbias[%]	Δr[%]
Pure water
a_w_Gregg	Gregg and Rousseaux (2016)	–	–	–	–	–	–	–	–	–	–
a_w_Mason	Mason et al. (2016)	–	40.6	40.0	−41.0	22.0	22.2	−15.0	–	6.7	–
a_w_Mason_TS	Mason et al. (2016) T‐S corrected	–	−5.0	−5.0	11.4	−14.0	−14.0	15.8	−5.0	–	–
Final	relative to the initial configuration	–	34.4	33.0	−32.4	5.0	6.0	−2.5	−5.0	6.7	–
Pure water IOPs and a_NAP
a_NAP_Case1_Chl	Bricaud et al. (2010)	Chl	‐	‐	‐	‐	‐	‐	‐	‐	‐
a_NAP_Case1_b_bp	Bricaud et al. (2010)	b_bp(700)	−12.5	−11.2	11.2	−10.0	−7.5	4.0	‐	−7.2	‐
a_NAP_Babin_Chl	Babin et al. (2003)	Chl	‐	−4.2	−2.9	−3.4	−8.0	−1.3	‐	‐	−1.1
a_NAP_Babin_b_bp	Babin et al. (2003)	b_bp(700)	−40.0	−39.2	22.4	−31.0	−34.8	9.1	−10.6	−15.4	1.10
Final	relative to the initial configuration	b_bp(700)	−47.0	−49.2	32.2	−39.4	−44.5	12.0	−10.6	−21.5	‐
Pure water IOPs plus a_CDOM
a_CDOM_Case1_Chl	Morel and Gentili (2009)	Chl	–	–	–	–	–	–	–	–	–
a_CDOM_Case1_fDOM	Morel and Gentili (2009)	fDOM	−26.9	−38.9	13.2	−26.9	−44.4	6.3	−11.0	−25.0	1.1
a_CDOM_K_bio_Morel	Organelli and Claustre (2019)	fDOM	−26.3	−18.2	9.1	−15.8	–	3.6	6.3	–	–
a_CDOM_K_bio_Mason	Organelli and Claustre (2019)	fDOM	−21.4	−33.3	3.5	−30.0	−18.8	2.3	–	10.0	–
a_CDOM_K_bio_Mason_a_w_corr	Organelli and Claustre (2019)	fDOM	−18.2	−16.7	2.3	−7.7	−14.3	1.1	−6.7	−11.1	1.1
Final	relative to the initial configuration	fDOM	−57.7	−66.7	27.9	−50.0	−61.1	12.7	−16.7	−25.0	1.1
Pure water IOPs plus separate IOPs: a_NAP, a_CDOM, a_ϕ and b_p
a_w_Mason_TS	Mason et al. (2016)	–	–	–	–	–	–	–	–	–	–
a_NAP_Babin_b_bp	Babin et al. (2003)	b_bp(700)	−60.5	−65.0	86.4	−53.5	−60.5	27.3	−19.0	−31.2	2.2
a_CDOM_K_bio_Mason_a_w_corr	Organelli and Claustre (2019)	fDOM	−79.1	−87.5	102.0	−72.1	−84.2	36.4	−33.3	−50.0	4.5
a_ϕ	Di Cicco et al. (2017)	Chl	−32.6	−37.5	52.3	−41.9	−50.0	22.7	−28.6	−43.8	4.5
b_p	Antoine et al. (2011)	b_bp(700)	−4.7	−5.0	11.4	−11.6	−13.2	9.1	−19.0	−31.3	2.2

The conversion from in‐water to above‐water remote sensing reflectance (hereafter ) follows the relationship from Lee et al. (2002):

Satellite data were obtained from Copernicus Marine Environment Monitoring Service the Ocean Color Level 3 products (OCEANCOLOUR MED OPTICS L3 REP OBSERVATIONS 009 095), comprising of data at 6 wavelengths: 412, 443, 490, 510, 555, 670 nm, as well as of the diffuse attenuation of downwelling irradiance at 490 nm, . Given the fact that no upwelling component of irradiance measurements is available from BGC‐Argo floats, a more in‐depth study of most appropriate scattering regimes is left for similar tests with data from multi‐spectral platforms as ProVal (Leymarie et al., 2018)

Locations of floats were matched‐up with daily satellite data of a 1 km grid space resolution. A total of 445 points were left for the period corresponding to the simulations considered. Due to a reduced number of matched‐up quantities, the values of and were spatially averaged to Western and Eastern Mediterranean basins, and temporally in the form of monthly climatological values.

Results and Discussion

IOP Model Validation

In order to verify the improvement of various modeling configurations, simulations were clustered into groups of separate IOPs, each with its own selection of tests and modifications. The model skill was quantified with three statistical parameters (root mean square error or RMSE, bias, and Pearson correlation coefficient ), resulting from a point‐by‐point match‐up of modeled and measured downwelling irradiance values for the first 150 m at three wavelengths.

Starting with a pure water IOP model, the updated absorption spectrum in the UV/blue range (, Mason et al., 2016) reveals a skill deterioration due to a much lower water absorption in the tested range of the spectrum. Most noticeably at shorter wavelengths, RMSE increases by 0.13, 0.09, and 0.01 (40.6%, 22%, and ) at 380, 412, and 490 nm, respectively, compared to the “reference” simulation (, with the pure water absorption spectrum used in Gregg & Rousseaux, 2016). Similarly, bias increases by 0.12, 0.08, and 0.01 (40%, 22.2%, and 6.7%) compared to the reference configuration. The correlation coefficient decreases by 0.26, 0.10, and 0.001 (41%, 15%, and ), respectively, Figure 2a and Table 2. values from Mason et al. (2016) are chosen as reference in subsequent simulations, assuming that the latest technology development enabled more accurate spectral measurements.

Figure 2

Bar plots resulting from a point‐by‐point match‐up with modeled and measured values. The three figures (from left to right) show bias, root mean square error, and Pearson correlation coefficient, respectively. Each simulation type has three bar plots, representing different wavelengths (purple: 380 nm, cyan: 412 nm, and orange: 490 nm): (a) Pure water IOPs. (b) Pure water IOPs and . (c) Pure water IOPs and . (d) Pure water plus separate IOPs: , , , and .

The inclusion of T and S data both for absorption (Sullivan et al., 2006) and scattering spectra (Zhang et al., 2009) displays a smaller, however notable improvement in the model’s skill. Compared to the simulation with absorption values from Mason et al. (2016), RMSE decreases by 0.02, 0.07, and 0.01 (5%, 14%, and 5%) and the bias by 0.02, 0.06, and (5%, 14%, and ). The last configuration (), with the modification from Mason et al. (2016) and T‐S corrected models of and , was therefore chosen for subsequent tests.

A series of tests were performed for parameterizations, Figure 2b and Table 2. Starting with a Case 1 model that follows the Chl profile shape ), the consecutive improvements incorporated the inclusion of depth variability (). Moreover, the range of and values from Babin et al. (2003) were tested with both Chl and shapes ( and ). Among the tests with a varying range of values ( between 0.0104 and 0.0178, and between 0.0087 and 0.08 , the latter corresponding to highly turbid waters), the minimum value for was chosen, assuming that floats are located in open waters with a low or negligible contribution of sediments. The slope is selected from a mean value of 0.0129 from Babin et al., 2003.

Considering vertical profile () instead of Chl () in the Case 1 configuration from Bricaud et al. (2010) significantly increases the skill, especially toward the UV, as the RMSE decreases by 0.04, 0.03, and 0.001 (12.5%, 10%, and less than 1%). The bias, on the other hand, decreases by 0.03, 0.02, and 0.01 (11.2%, 7.5%, and 7.2%). Shifting toward non‐Case‐1 representations, with the inclusion of the range of values observed in in situ measurements, gives an overall better match‐up statistics, which especially improves when considering the vertical profile (), Figure 2b. Comparing the ‐shaped model with values from Babin et al. (2003) () and the analogous Case 1 model (), RMSE decreases by 0.11, 0.11, 0.02 (40%, 33%, and 10.6%—values not directly shown in Table 2), and the bias by 0.1, 0.1, and 0.02 (42%, 40%, and 15%) for the three measured wavelengths, respectively. increases by 0.13 (19%), 0.07 (7.6%), and at 490 remains the same. Therefore, according to the present data, the best agreement is achieved using the vertical profile with the model suggested by Babin et al. (2003), that is, the model.

Similarly to NAP simulations, CDOM absorption models were also compared considering three aspects: the Case 1 versus alternative parameterizations, Chl () versus fDOM IOP depth variability (), and additional spectral corrections due to modifications in the pure water spectrum shown in Equation 15. As in , considering shapes alternative to Chl, such as profiles of fDOM, reveals a drastic improvement in the match‐up statistics. fDOM‐shaped Case 1 model from Morel and Gentili (2009) () introduces a RMSE decrease 0.07, 0.07, and 0.02 (26.9%, 26.9%, and 11.0%) and a reduction of bias amounting to 0.07, 0.08, and 0.03 (38.9%, 44.4%, and 25%). increases by 0.09, 0.05, and 0.01 (13.2%, 6.3%, and 1.1%), Figure 2c and Table 2. As in Figure 2b, the significant impact on the lowering of bias and RMSE values was also due to a deviation from Case 1 models. This was achieved by adopting the approach presented in Organelli and Claustre (2019), described in Section 2.3.3, with the difference that the first optical depth range was rather considered as it resulted in a better performance compared to the MLD (not shown). Relative to the fDOM‐shaped Case 1 model (), in the RMSE decreases by 0.05 and 0.03 (26.3% and 15.8%) at 380 and 412 nm and increases by 0.01 (6.3%) at 490 nm, respectively. Subsequent simulations result in an upgrade in the calculation of : from the original value of 0.0151 (Morel & Maritorena, 2001), was calculated by taking into consideration the T‐S corrections for both absorption and scattering values (). Moreover, the was modified for the spectral correction of (). Compared to the constant value simulation (), the final configuration resulted in a decrease in RMSE by 0.05, 0.04, and 0.03 (35.7%, 25.0%, and 6.7%) and in a bias decrease by 0.04, 0.04, and 0.02 (44%, 40%, and 20.0%) at 380, 412, and 490 nm, respectively (not directly shown in Table 2).

The contribution of remaining IOPs, phytoplankton absorption and scattering by particles , are shown alongside the skill of the chosen models for separate IOP groups described above, Figure 2d and Table 2. The PFT modeling configuration described in Section 2.3.4, compared to the pure water simulation results, resulted in a RMSE decrease by 0.14, 0.18, and 0.06 (32.6%, 41.9%, and 28.6%) and in a bias decrease by 0.15, 0.19, and 0.07 (37.5%, 50%, and 43.8%). The correlation increased by 0.23, 0.15, and 0.04 (52.3%, 22.7%, and 4.5%). Based on the phytoplankton absorption curves adopted in the model, the highest decrease in RMSE and bias at 412 nm can be explained by the proximity to the chlorophyll absorption peak in the blue, which can also explain a more uniform spectral change of skill. Moreover, the absorption values of most PFTs (except Cryptophytes and Synechococcus) are similar at 380 and 490 nm, with slightly higher values at 380 nm.

Even though several configurations were tested, their impact on the match‐up was negligible, leading to small differences between simulation results. The chosen scattering model was a non‐Case 1 type derived from measurements following Equations 20 and 21, with a maximum backscattering ratio of 0.015 and a spectral slope of 3. The selection of these values is motivated by examining and climatologies as discussed in Section 3.2.

The final modeling configuration, including all optically significant constituents considered in this study, that is results in a RMSE ranging from 0.05 to 0.09 , a negative bias of −0.02 at 412 nm and positive values of 0.01 at 380 and 490 nm, while is 0.93 at 380 nm, 0.94 at 412 nm, and 0.95 at 490 nm. Figure 3. The slope is closest to 1 at 490 nm, with the highest value observed at 380 nm (1.03), signifying a model overestimation, and lowest slope at 412 nm (0.78), with model values lower than float measurements.

Figure 3

Match‐up with the final modeling configuration including all IOPs in the first 150 m, displaying the root mean square error (RMSE), bias, Pearson correlation coefficient (r), slope, intercept (Y‐int), and number of points (N) for each of the three wavelengths considered.

Examples of the model results for west and east are displayed in Figures 4a and 4b, respectively. The top row shows model forcings that were used for IOP parameterization and their depth variability (Chl, and fDOM), whereas the bottom row displays both the model output and radiometric measurements (). Modeled and measured irradiance values are quite in agreement, both in terms of vertical shapes and the first optical (e‐folding) depth ranges at all three wavelengths considered. As the model is configured in such way that it takes the IOP depth variability based on the local float biogeochemical measurements, no regional bias was observed in the results. This is quite encouraging given the fact that two different bio‐optical regimes are clearly at play in the Western and Eastern Mediterranean. More specifically, the Western Mediterranean subsurface Chl maximum overlaps with the shape of , which might suggest that the modifications in Chl are related to actual changes in phytoplankton biomass and community structure, Figure 4a (Barbieux et al., 2018). On the contrary, no such co‐variance is observed for the Eastern Mediterranean, and values of are much higher, which might be explained more due to physiology and photoadaptation than to changes in the actual biomass, Figure 4b (Barbieux et al., 2018). Moreover, in the Eastern Mediterranean the particles might be also of more mineral origin due to episodic dust deposition events (Claustre et al., 2002). Such findings clearly demonstrate the importance of having synoptic biogeochemical, bio‐optical and radiometric measurements in order to validate the IOP metrics with available irradiance profiles. Furthermore, simulations show that the inclusion of additional parameters (i.e., fDOM and ) in the IOP parameterizations results in a significant improvement of the match‐up statistics compared to Chl‐shaped (Case 1) IOP models.

Figure 4

Example of a BGC‐Argo profile with modeled and measured radiometric values. Pale dots are values prior to the QC additionally applied to this study. (a) Western Mediterranean. (b) Eastern Mediterranean.

Comparison With In Situ and Remote Sensing Apparent Optical Properties

Results were further assessed also in terms of diffuse attenuation coefficients for downwelling irradiance and remote sensing reflectance . Both AOPs were calculated for the first optical (i.e., e‐folding) depth range, following methods described in Sections 2.3 and 2.4. Additionally, for the wavelength 490 nm, a three‐platform comparison was possible with model‐ and float‐derived versus satellite data. The match‐up of satellite and float observations amounted to 445 co‐located measurements which were spatio‐temporally aggregated into climatological months and grouped according to western and eastern basins. More information on the number of profiles per month per region is shown in Table 3.

Table 3

Number of Matched‐Up Profiles per Month per Region (i.e., Western and Eastern Mediterranean)

month	West	East
1	58	42
2	49	48
3	86	69
4	58	70
5	75	110
6	43	48
7	32	38
8	19	26
9	22	16
10	17	22
11	31	25
12	70	42

Seasonal Taylor diagrams of , divided in west (darker points) and east (lighter points), for RMSE and values are shown in Figure 5. The aim was to assess the impact of changing the relative contribution of within the term as described in Section 2.3.3. The factor was plotted for the extreme ranges of tested values: (opaque) and (transparent). In terms of the scattering modeling configurations, there was no difference in the skill between the two ranges of backscattering ratio, hence only is shown in Figure 5. Different values of slopes are displayed with different markers. Results convey two distinct clusters for the western and eastern basins, which could imply regionally different bio‐optical regimes, with RMSE values always lower for the Eastern Mediterranean. Such findings are in line with Terzić et al. (2019), which also shows zonal gradients in modeled and float‐derived values. Moreover, changing values does not seem to have a strong impact on the skill in terms of statistics, as points are in most cases concentrated. factor proves to have the largest impact on the model skill (f = 0.5), especially at shorter wavelengths during winter, when RMSE is reduced for more than a factor of 5 at 380 nm. Even though worsens the model match‐up in terms of (), as shown in Figure 2c and Table 2, this might suggest that optically active constituents absorbing in the UV/blue (such as CDOM) are more important at greater depths, as is assessed only at the first optical depth. This might in turn explain why works better at greater depths, where fDOM most likely has a larger impact.

Figure 5

Taylor diagrams of seasonally divided (rows) values of RMSE of at 380, 412, and 490 nm (columns). Dark colors represent the Western Mediterranean, while lighter colors depict the Eastern Mediterranean. Transparent points indicate values of multiplied by a factor of 1 and the opaque ones stand for a factor of 0.5. The symbols show different values of the slope .

Modeled coefficients replicated the monthly dynamics computed from float measurements, Figure 6. At 380 nm, maximum discrepancy is seen in winter and spring months for western basins, with a largest difference in the month of April, with mean values of 0.125 and 0.11 for model and data, respectively (top figure in Figure 6). At longer wavelengths the difference diminishes, with good consistency achieved also between the three comparing platforms at 490 nm. However, satellites do not seem to capture highest values in spring for the Western Mediterranean, which is shown from BGC‐Argo floats and model results. Overall, values are larger for western than eastern basins, as shown in Terzić et al. (2019). Modeling IOPs as functions of available biogeochemical and bio‐optical measurements therefore provides a significant reproduction of the zonal gradients. The similar magnitude of error bars from all platforms demonstrates also that the model and data variabilities are close.

Figure 6

Monthly climatology of values with , where . The top two figures display model (purple points) and float (orange points) values for western (darker color) and eastern (lighter color) basins. At 490 nm, additional green scatter points from satellite data are included (bottom figure).

Considering another AOP, however, a different result is obtained. is related to IOPs in such way that is directly proportional to backscattering and inversely proportional to the sum of absorption and backscattering (Morel & Gentili, 1993; Morel & Prieur, 1977). Raman scattering was accounted for by correcting values according to Lee et al. (2013), and its inclusion amounts to 3.5%, 4%, and 8.5% difference in terms of at 412, 443, and 490, respectively (not shown).

Seasonal Taylor diagrams of , divided in west (darker points) and east (lighter points), for RMSE and values are shown in Figure 7. Unlike in Figure 5, no zonal gradients are observed, instead, points seem to be quite dispersed. is generally less than 0.6 during most seasons, with exceptions seen in autumn (at 443 nm) and winter (at 490 nm), up to a maximum of around 0.8 and 0.7, respectively. In terms of the factor ( displayed as opaque and as transparent), the greatest impact is seen at higher slope values (η = 4), for which RMSE decreases at all seasons and wavelengths, however still having highest values among all . As in Figure 5, changing the factor has a greater impact than changing the backscattering ratio, as the range of tested values is quite small (between 0.2% and 1.5%), and thus only values of are shown in Figure 7.

Figure 7

Taylor diagrams of seasonally divided (rows) values of RMSE of at 412, 443, and 490 nm (columns). Dark colors represent the Western Mediterranean, while lighter colors depict the Eastern Mediterranean. Transparent points indicate values of multiplied by a factor of 1 and the opaque ones stand for a factor of 0.5. The symbols show different values of the slope .

Given the lack of in situ upwelling radiometric measurements, as well as the uncertainty of remote sensing in the blue part of the spectrum, no definite conclusions can be placed on the most adequate scattering model parameters. However, using Case 1 from Equation 19 leads to an underestimation of modeled for both west and east, resulting in up to a 60% discrepancy with satellite values, even when using the shape, and especially during summer months (not shown). By rather focusing on the ‐derived scattering models while looking at monthly climatological data shown in Figures 8a and 8b, certain range of slope values can be preferred for certain seasons and regions. Lower values () seemed to work best during winter/early spring in the Western Mediterranean (seen from Figure 7), with for the spring in the Eastern Mediterranean, as well as for late spring and summer at west, Figure 8a. resulted in a better agreement with remote‐sensing data for summer and autumn in the Eastern Mediterranean, Figure 8b. Slopes of 2 and 3 are also consistent with the range of values from Antoine et al. (2011). This finding might suggest that there are different scattering regimes at play in the two basins, most likely stemming from a different particle size distribution (Antoine et al., 2011), which can provide information also on the dominant phytoplankton (Kostadinov et al., 2009; Organelli et al., 2020). Lower slope values thus imply larger particles, which is consistent with the results in the west during usual spring bloom events with larger, microphytoplankton assemblages (20–200 ). On the other hand, higher slope values could suggest smaller particles, consistent with the pico‐ or nanophytoplankton (0.2–2 and 2–20 , respectively) assemblages usually predominant at the basin level, with the former prevailing especially during spring/summer and the latter during winter. Such conclusions are in line with the previously detected patterns of phytoplankton distribution in Siokou‐Frangou et al. (2010) and Uitz et al. (2012), which were confirmed also by Sammartino et al. (2015), Di Cicco et al. (2017), and Navarro et al. (2017).

Figure 8

Monthly climatology of values with . All three figures display model (purple points) and satellite (orange points) values for western (darker color) and eastern (lighter color) basins at 412, 443, and 490 nm, respectively. (a) Where and . (b) Where and .

Conclusions

BGC‐Argo floats prove to be an essential observing system to further explore the possibility of integrating data in numerical modeling of physical, as well as biogeochemical and optical properties. Due to the high number of profiles with synoptic measurements of physical and bio‐optical parameters, it is possible to use the almost complete suite of measured variables (T, S, Chl, fDOM, , and ) to test various state‐of‐the‐art parameterizations of absorption and scattering properties of sea water constituents. The current wavelength selection of measurements constitutes an ideal tool to explore the part of the spectrum that is least understood, mostly for the contributions from CDOM and NAP. This is particularly true in the Mediterranean Sea, where the blue‐to‐green reflectance ratio‐based algorithms are known to have low performances (Morel & Gentili, 2009) because of the higher‐than‐expected contribution of CDOM for a given Chl concentration. The major findings of this work can be summarized as follows: the inclusion of T and S data is recommended to account for the small, but significant spectral modulation of seawater compared to pure water, which also improves the model skill. Furthermore, the tests performed on Case 1 IOP models reveal that the inclusion of additional biogeochemical measurements in the IOP parameterizations results in improved match‐up statistics, both when comparing with irradiance profiles, as well as with in situ and remote‐sensing derived AOPs. The shape of for variability increases the skill compared to Chl‐shaped models by up to 40% in the case of RMSE. Moreover, it was demonstrated that the use of fDOM shape and the estimation of for CDOM absorption from , as well as the spectral correction of the updated spectrum, all contribute to an upgrade in CDOM modeling of up to 57.7% in terms of RMSE. Different relative contributions of as an indicator of were shown for different metrics, matching up values at a 150 m depth range versus at the first optical depth. Results implied a lower relative contribution of to () at shallower depths, and a higher one () at greater depths, suggesting a major importance of CDOM dynamics also at depths which cannot be captured by satellites. Therefore, partitioning the contributions of NAP and CDOM to the total absorption with additional experiments would be also advantageous, as well as the assessment of relative contributions of different constituents to the total signal, thus separating the organic and inorganic parts.

The inclusion of PFTs has demonstrated the importance of accounting also for phytoplankton, resulting in a more uniform spectral response in the blue, decreasing the RMSE up to 41.9% compared to the pure water simulations. Observations of other biogeochemical parameters, such as oxygen, nitrate, and pH, can be possibly integrated with a coupled biogeochemical model. All of these variables are already available on the BGC‐Argo float standard configuration (Claustre et al., 2020). This could offer the opportunity, with an existing validation data set, to consider also the phytoplankton ecology and dynamics of separate functional groups. Such work demonstrates the advantages of combining data with numerical models, which can pave way to a better understanding of biogeochemical processes in the examined regions.

The focus of this study is also more on the absorption models rather than scattering due to the lack of measurements and the uncertainty of remote sensing in the blue part of the spectrum. However, despite these limitations, the model is still able to capture the spatio‐temporal variability of slope values, indicating different phytoplankton and particle size distributions. With the integration of multi‐spectral data from platforms like ProVal (Leymarie et al., 2018), it will be possible to further examine this issue in more detail. This will enable also the calculation of in situ and remote‐sensing reflectance estimations , thus surpassing the current limitation of quantifying the skill between satellite data with model values due to the scarcity of satellite observations spatio‐temporally co‐located with BGC‐Argo float profiles, supporting further the three‐platform comparisons.

To conclude, the key point raised in this study is that the inclusion of multi‐spectral measurements is essential to tackle the proper biogeochemical response, surpassing the most‐commonly PAR‐related parameterizations of phytoplankton growth. With the advancement of satellite sensors and their algorithms it would be necessary to make a comparison of radiative transfer models of different degrees of complexity, and perform similar tests with hyperspectral models which are able to solve a full radiative transfer equation resulting in a complete radiance distribution (Hedley et al., 2020).

Erratum

In the originally published version of this article, the third group of simulations in Table 2 (Pure water IOPs plus a) inadvertently repeated the same numbers and simulation names as the fourth group of simulations (Pure water IOPs plus a, a and a). This has been corrected, and this may be considered the official version of record.