Retrieval of Particulate Backscattering Using Field and Satellite Radiometry: Assessment of the QAA Algorithm

: Particulate optical backscattering (b bp ) is a crucial parameter for the study of ocean biology and oceanic carbon estimations. In this work, b bp retrieval, by the quasi-analytical algorithm (QAA), is assessed using a large in situ database of matched b bp and remote-sensing reﬂectance (R rs ). The QAA is also applied to satellite R rs (ESA OC-CCI project) as well, after their validation against in situ R rs . Additionally, the e ﬀ ect of Raman Scattering on QAA retrievals is studied. Results show negligible biases above random noise when QAA-derived b bp is compared to in situ b bp . In addition, R rs from the CCI archive shows good agreement with in situ data. The QAA’s functional form of spectral backscattering slope, as derived from in situ radiometry, is validated. Finally, we show the importance of correcting for Raman Scattering over clear waters prior to semi-analytical retrieval. Overall, this work demonstrates the high e ﬃ ciency of QAA in the b bp detection in case of both in situ and ocean color data, but it also highlights the necessity to increase the number of observations that are severely under-sampled in respect to others environmental parameters.


Introduction
Retrieval of water inherent optical properties (IOPs) from both field and ocean color radiometry is at the base of several biogeochemical and physical oceanographic studies [1,2]. IOPs of algal and non-algal particles can be derived from remote sensing reflectance spectra (R rs ; units of sr −1 ) by using appropriate algorithms [3][4][5]. Among IOPs, the particulate optical backscattering coefficient (b bp ; in m −1 ) is related to the particle concentration in seawater, on their size distribution, refractive index, shape and structure [6][7][8]. Former research suggested that b bp is mostly influenced by submicron non-algal particles [9][10][11]. However, it has been recently shown that most of b bp is due to particles with equivalent diameters between 1 and 10 µm [8], thus including the contribution of phytoplankton cell and supporting the use of b bp for the retrieval of: (i) particulate organic concentration (POC) [12,13]; (ii) particle size distribution [14,15]; and (iii) phytoplankton carbon biomass concentration (C phyto ; mg m −3 ) [16][17][18], a reference wavelength of λ = 555 nm and then the calculation is extended to other wavelengths by assuming a power law b bp = b bp (λ 0 )(λ/λ 0 ) −η for the b bp with a spectral slope. η = p 1 1 − p 2 exp −p 3 r rs (443) r rs (555) (1) Equation (1) is widely used for QAA retrievals of b bp at multiple wavelengths. Nevertheless, we use the in situ dataset presented here to evaluate the accuracy of analytical η. The functional form of Equation (1) is used and the default numerical coefficients p 1 = 2.0, p 2 = 1.2 and p 3 = 0.9 [22] are replaced by unknown variables established by non-linear regression. To this aim, we used the iterative bi-square method, which minimizes a weighted sum of squared errors, where the weight given to each data point decreases with the distance from the fitted curve [37]. This procedure makes the curve sensitive to the bulk of the data and the effect of outliers is reduced. The error function is minimized through the trust region algorithm [38]. In addition, the 95% confidence prediction bounds are also computed.
It is known that for oligotrophic waters the Raman scattering plays a significant role that is not accounted for in the semi-analytical R rs modeling [39]. Therefore, a pertinent question is how much this phenomenon affects the semi-analytical b bp retrievals. With this aim, Lee et al. [40] developed empirical correction formulas to compensate the Raman scattering on the R rs . Here, we assess the effects of this compensation on the difference between the in situ b bp and R rs -derived b bp . The statistical assessment was also replicated in other different cases: (i) validation of satellite CCI R rs against in situ R rs ; and (ii) b bp retrievals after the application of QAA to satellite R rs (Raman corrected), were compared to in situ measurements at the different available wavelengths.
Estimated data y i are compared to reference data x i by using the following statistical indexes: relative bias (units of %), relative root-mean square error (RMS, units of %) and determination coefficient (r 2 )

In Situ Data
The in situ database is composed of three distinct datasets containing multi-spectral R rs and b bp : the recently updated global in situ database ( [33] hereafter V19 dataset), an in situ dataset collected by the Italian National Research Council (CNR) during two field campaigns in the Mediterranean Sea ( [29] hereafter CNR dataset) and the time-series of data acquired by the BOUSSOLE buoy in the northwestern Mediterranean Sea ( [34,41]; hereafter BOU dataset [42]). The three in situ databases were quality-checked as described below. All the R rs data were band-shifted to the CCI bands (those of the NASA SeaWiFS instrument, namely 412, 443, 490, 510, 555, and 670 nm). The band-shifting procedure [43] is a technique to compensate small band differences in multispectral R rs data. It takes into account the spectral shape of the absorption and scattering that contribute to R rs and constitutes a more accurate approach than a simple linear interpolation. Considering every wavelength an independent measurement, the final dataset accounts for a total of N = 2881 R rs and b bp co-located measurements around the global ocean ( Figure 1). As shown in Figure 2, the total R rs and b bp spectra cover from oligotrophic open-ocean to more eutrophic coastal waters as the range of R rs and b bp values vary between 0-0.02 sr −1 and 10 −4 -10 −1 m −1 respectively. Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 16 Figure 1. Geographical distribution of the in situ Rrs vs. bbp matchups. Some areas (A1, A2, and A3) concentrate a high point density and are highlighted in zoomed maps. Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively.

Figure 2.
Rrs and bbp spectra for the three-different datasets: V19, BOU, and CNR. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.

V19 Dataset
Rrs and IOPs, aggregated within ±6 nm, were downloaded. V19 is a global compilation of in situ data that was acquired from many sources (e.g., MOBY, AERONET-OC, SeaBASS, NOMAD, MERMAID, AMT, and many others), motivated by the validation of the ocean-color products from the ESA OC-CCI products. Methodologies were implemented for homogenization, quality control  Rrs and bbp spectra for the three-different datasets: V19, BOU, and CNR. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.

V19 Dataset
Rrs and IOPs, aggregated within ±6 nm, were downloaded. V19 is a global compilation of in situ data that was acquired from many sources (e.g., MOBY, AERONET-OC, SeaBASS, NOMAD, MERMAID, AMT, and many others), motivated by the validation of the ocean-color products from the ESA OC-CCI products. Methodologies were implemented for homogenization, quality control 2.2.1. V19 Dataset R rs and IOPs, aggregated within ±6 nm, were downloaded. V19 is a global compilation of in situ data that was acquired from many sources (e.g., MOBY, AERONET-OC, SeaBASS, NOMAD, MERMAID, AMT, and many others), motivated by the validation of the ocean-color products from the ESA OC-CCI products. Methodologies were implemented for homogenization, quality control and merging of all data. No changes were made to the original data, other than averaging of observations that were close in time and space, elimination of some points after quality control and conversion to a standard format [44].
In this study, data were selected only if valid and corresponding R rs and b bp measurements at all CCI bands were available. Such condition determines a total of N = 319 matchups. Remaining minor b bp wavelength mismatches were removed by linear interpolation to the CCI bands. Although V19 is a merged dataset from multiple datasets, the condition we set for the matchup left data that were originally from the NOMAD dataset only.

BOU Dataset
The BOUSSOLE (BOUee pour l'acquiSition d'une Série Optique a Long termE) project started in 1999, and its activities are developed on a site located in the northwestern Mediterranean Sea (7 • 54 E, 43 • 22 N, Figure 1, panel A2). Essential information about the site characteristics, the measurement platform, and the instrumentation was previously documented [34,41,42]. The b bp data were collected at 9 m nominal depth with a Hobilabs Hydroscat-4 (442, 488, 550, and 620 nm) and processed as in [45]. In addition, a quality control on b bp was applied that required a spectral b bp slope, calculated from every pair of two consecutive bands, within given bounds (more than −1 and less than 6). R rs data were derived with a set of Satlantic 200-series multispectral radiometers ( [46] and references therein). The R rs is available at a varying number of the following bands, depending on the time period: 412.5, 442.5, 490, 510, 555, 560, 665, 670, and 681.25 nm. Since the application of the QAA requires R rs at 443, 490, 555, and 670 nm, only R rs records whose native bands matched those needed by the QAA algorithms were selected (within a ±6 nm range). Data at 412.5 nm and 442.5 nm were band-shifted to 412 nm and 443 nm [43], respectively. In the green region, if the R rs at 555 nm was available, it was directly sampled and the R rs at 560 nm was ignored. If the R rs at 560 nm was available when the R rs at 555 nm was missing, the R rs at 560 nm was band-shifted to 555 nm. Similarly, in the red region, between the R rs at 665 nm and 670 nm, preference to R rs at 670 nm was given. R rs data at 681.25 nm was not considered for the analysis. Data was generally available within two hours from the local noon. The time series at sub-daily resolution were reduced by calculating the daily medians.

CNR Dataset
Data belong to two field campaigns conducted in 2013 and 2015 in Italian seas, encompassing a high optical range between open and coastal waters. Measurements were performed between 8:30 h and 16:00 h UTC. IOPs and R rs were collected sequentially at each station, with a maximum delay of 1 h and ship drift of maximum few hundreds of meters.
Backscattering was measured with an ECO-VSF3, manufactured by WET Labs, Inc., at the wavelengths 470, 530, and 660 nm. This instrument measures the volume scattering function at three backward angles and calculates b b by integration of a polynomial fit. Final data are the result of a binning across the first optical depth.
Radiometry was performed with OCR-507 radiometers, manufactured by Satlantic, Inc., measuring at the center bands 412, 443, 490, 510, 556, 665, and 865 nm. In-water upwelling radiance at nadir (L u ) sensor was mounted onto a free-falling T-shaped structure, with the multicast technique. Above-water downwelling irradiance (E s ) data were collected by a reference sensor, mounted at the top of the ship's deck. R rs was computed using the SERDA software developed at CNR [47]. All the R rs data were band-shifted to the CCI bands for consistency with the satellite R rs . Further details about this dataset are provided in Pitarch et al. [29].

Satellite ESA OC-CCI R rs Data
The ESA OC-CCI version 4.0 global daily R rs data at 4 km resolution for the period 1997-2017 were downloaded [48]. CCI products are the result of the merging of SeaWiFS, MERIS, MODIS, and VIIRS data in which the inter-sensor biases are removed [49]. This version 4.0 includes the latest NASA reprocessing R2018.0 that mostly accounts for the aging of the MODIS sensor. ESA OC-CCI provides the daily R rs data and associated uncertainty maps in terms of bias and RMS, which were generated with a procedure that included comparison to in situ data and optical water type analysis [48].
In this work, a conservative extraction procedure was followed, in which the center R rs data within a 3 × 3 pixels box was extracted only if all the 9 pixels were not flagged, therefore minimizing possible land border, cloud or other environmental contaminations, and obtaining the highest quality of matchups. Finally, for each single R rs , the bias was also extracted and then compensated pixel-by-pixel.

QAA Performance for b bp Retrievals from In Situ Data
QAA-retrieved b bp from in situ R rs is here compared to in situ measured b bp . Comparisons are made at the native bands of each in situ b bp instrument for the cases of BOU and CNR and at the CCI bands for V19. Statistics are also presented for each band and dataset.
A first assessment consists of applying the QAA without performing the compensation for Raman scattering. Here, results show a general overestimation of around +43.4% for V19 ( Table 1) that is not significant given the overall noise expressed by the RMS (152%). This high RMS is the likely consequence of the different protocols, instrumental and geophysical noises affecting all single contributors to the V19 dataset (Table 1). In the case of the BOU dataset, an overall overestimation of +49.2% is found for all the bands which is statistically significant given the related RMS (58.7%). On the other hand, the QAA applied to the CNR dataset showed the highest performances, with a bias of +3.3% and a RMS below 23%. Table 1. Statistical descriptors of the difference between the QAA-derived b bp and in situ b bp for each dataset, without Raman scattering compensation. Figure A1 provides a graphical representation of this To understand the importance of the Raman scattering correction in semi-analytical b bp retrievals, the analysis is repeated with corrected R rs [40]. The application of the Raman scattering correction reduces both bias and RMS nearly for all the b bp at all bands (Table 2 and Figure 3). Indeed, for the V19 dataset, the bias decreases to 12% with respect to the retrievals obtained without correction of the Raman scattering ( Table 2). The RMS reduction is around 34%. For the BOU data, the RMS and bias improve of about 11% and 12%, respectively. In the case of the CNR dataset, statistics show a modest increase in accuracy except for λ = 660 nm, which is likely influenced by chlorophyll-a fluorescence.
Although fluorescence peaks at around λ = 660 nm, the ECO-VSF 3 sensor, used to collect the CNR dataset, has a full width at half maximum (FWHM) of about 20 to 30 nm, so a fluorescence interference may not be excluded.
Overall, these results are somewhat expected as the Raman scattering correction produces a smaller effect in coastal waters [50], which represent a significant part of the CNR dataset with respect to the two other datasets (Figure 2). The overall statistics are in agreement with previous comparisons that showed negligible biases over noise at global scale [5] and at regional level [29]. Results in this section highlight the importance of applying the Raman scattering correction to the source R rs prior to semi-analytical b bp retrieval in order to increase the accuracy. Table 2. Statistical descriptors of the difference between the b bp -QAA derived and in situ b bp for each dataset with Raman scattering compensation. Figure A2 provides a graphical representation of this sensor, used to collect the CNR dataset, has a full width at half maximum (FWHM) of about 20 to 30 nm, so a fluorescence interference may not be excluded.
Overall, these results are somewhat expected as the Raman scattering correction produces a smaller effect in coastal waters [50], which represent a significant part of the CNR dataset with respect to the two other datasets (Figure 2). The overall statistics are in agreement with previous comparisons that showed negligible biases over noise at global scale [5] and at regional level [29]. Results in this section highlight the importance of applying the Raman scattering correction to the source Rrs prior to semi-analytical bbp retrieval in order to increase the accuracy. Table 2. Statistical descriptors of the difference between the bbp-QAA derived and in situ bbp for each dataset with Raman scattering compensation. Figure A2 provides a graphical representation of this

Estimation of the bbp Spectral Slope from Rrs Data
The in situ dataset described above is used (see Section 2.2) to assess the proposed relationship in the QAA and perform a model to data fit that is compared to the common QAA v6 equation [50]. Figure 4 shows a comparison of the independent variable (i.e., the blue-to-green band ratio rrs(443)/rrs(555)) with respect to derived from the in situ bbp. Fitting a functional form of Equation

Estimation of the b bp Spectral Slope from R rs Data
The in situ dataset described above is used (see Section 2.2) to assess the proposed relationship in the QAA and perform a model to data fit that is compared to the common QAA v6 equation [50]. Figure 4 shows a comparison of the independent variable (i.e., the blue-to-green band ratio r rs (443)/r rs (555)) with respect to η derived from the in situ b bp . Fitting a functional form of Equation (1) returns a curve (p 1 = 2.2, p 2 = 0.9 and p 3 = 0.5) and a 95% prediction interval, which is around ±1 wide, caused by the high scatter of the data cloud. The difference between η computed here and the one derived via QAA is much smaller than the width of the prediction interval, thus making them equivalent for prediction purposes. Therefore, by the principle of parsimony, the operational η functional form (dashed line in Figure 4) remains valid. However, one must keep in mind that the low predictive value of this relationship may result in b bp extrapolations to bands outside the reference one (usually 555 nm) that accumulate significant errors. In particular, within a worst-case scenario, an error in η estimation, ∆η = 1, will lead to a~26% error when extrapolating b bp from 555 nm to 412 nm.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 16 (1) returns a curve (p1 = 2.2, p2 = 0.9 and p3 = 0.5) and a 95% prediction interval, which is around ±1 wide, caused by the high scatter of the data cloud. The difference between computed here and the one derived via QAA is much smaller than the width of the prediction interval, thus making them equivalent for prediction purposes. Therefore, by the principle of parsimony, the operational functional form (dashed line in Figure 4) remains valid. However, one must keep in mind that the low predictive value of this relationship may result in bbp extrapolations to bands outside the reference one (usually 555 nm) that accumulate significant errors. In particular, within a worst-case scenario, an error in η estimation, Δη = 1, will lead to a ~26% error when extrapolating bbp from 555 nm to 412 nm.  (1) to all the data (p1 = 2.2, p2 = 0.9 and p3 = 0.5). The 95% confidence prediction bounds are represented by the grey shaded area. The dashed curve is the estimation from Rrs as defined in Equation (1). Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively.

Validation of CCI Rrs
Prior to applying to satellite data an algorithm that has been developed with in situ data, assessing the quality of the satellite Rrs with respect to in situ measured Rrs is desirable in order to identify possible biases. Therefore, this section uses the in situ Rrs contained in the three datasets to evaluate the CCI Rrs. There is a total of 882 matchups for V19, 581 for BOU, and 252 for CNR. Good agreement between in situ values and the CCI Rrs products is found ( Figure 5, Table 3) at all wavelengths, rather consistently with other previous results [5]. Overall, all datasets display similar performance, with negligible biases with respect to the overall noise expressed by the RMS. In the case of λ = 670, increased RMS is mostly due to the low values Rrs, except for CNR, that contains a higher data range. It is concluded that the CCI Rrs do not require adjustments at the studied wavelengths.
The magnitude of this RMS expresses a high bound for the overall uncertainty of the Rrs product as it is a measure of the errors in the comparison experiment, including those within the in situ data. The fraction of this error which is attributable to the satellite data only is likely to be much lower. To  (1) to all the data (p 1 = 2.2, p 2 = 0.9 and p 3 = 0.5). The 95% confidence prediction bounds are represented by the grey shaded area. The dashed curve is the η estimation from R rs as defined in Equation (1). Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively.

Validation of CCI R rs
Prior to applying to satellite data an algorithm that has been developed with in situ data, assessing the quality of the satellite R rs with respect to in situ measured R rs is desirable in order to identify possible biases. Therefore, this section uses the in situ R rs contained in the three datasets to evaluate the CCI R rs . There is a total of 882 matchups for V19, 581 for BOU, and 252 for CNR. Good agreement between in situ values and the CCI R rs products is found ( Figure 5, Table 3) at all wavelengths, rather consistently with other previous results [5]. Overall, all datasets display similar performance, with negligible biases with respect to the overall noise expressed by the RMS. In the case of λ = 670, increased RMS is mostly due to the low values R rs , except for CNR, that contains a higher data range. It is concluded that the CCI R rs do not require adjustments at the studied wavelengths.
The magnitude of this RMS expresses a high bound for the overall uncertainty of the R rs product as it is a measure of the errors in the comparison experiment, including those within the in situ data. The fraction of this error which is attributable to the satellite data only is likely to be much lower.
To have a measure of this fraction, a comparison to global in situ dataset with a traceable uncertainty budget would be desirable, though such option is presently not available.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 16 have a measure of this fraction, a comparison to global in situ dataset with a traceable uncertainty budget would be desirable, though such option is presently not available.

QAA Performance for bbp Retrievals from CCI Data
After assessing the quality of the QAA retrievals with in situ bbp and the quality of the CCI Rrs respect to in situ Rrs, the QAA is applied to the CCI Rrs to retrieve the bbp that are then compared to the in situ data. In agreement with our findings in Section 3.1, CCI Rrs are corrected for Raman Figure 5. Scatter plots of CCI R rs versus in situ R rs for the six different wavelengths. The dashed line represents the 1:1 ratio. Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively. Table 3. Statistical descriptors of the difference between satellite CCI R rs and in situ R rs for each dataset. Figure A3 provides a graphical representation of this

QAA Performance for b bp Retrievals from CCI Data
After assessing the quality of the QAA retrievals with in situ b bp and the quality of the CCI R rs respect to in situ R rs , the QAA is applied to the CCI R rs to retrieve the b bp that are then compared to the in situ data. In agreement with our findings in Section 3.1, CCI R rs are corrected for Raman scattering. Results of this comparison are shown in Figure 6 and Table 4. For V19, biases are not significant (less than 30%) in comparison of RMS values (less than 60%). On the other hand, similarly to the statistics derived from Section 3.1, QAA-derived b bp , as compared to the BOU data displays significant positive biases. Comparison with CNR data shows the highest performances, with bias of +2.7% and RMS of 48%. The conclusions from our analysis are consistent with previous comparisons to QAA, reporting negligible biases above noise level both at global and regional scales [25,30,51].
scattering. Results of this comparison are shown in Figure 6 and Table 4. For V19, biases are not significant (less than 30%) in comparison of RMS values (less than 60%). On the other hand, similarly to the statistics derived from Section 3.1, QAA-derived bbp, as compared to the BOU data displays significant positive biases. Comparison with CNR data shows the highest performances, with bias of +2.7% and RMS of 48%. The conclusions from our analysis are consistent with previous comparisons to QAA, reporting negligible biases above noise level both at global and regional scales [25,30,51]. Figure 6. Scatter plots of QAA derived bbp from CCI Rrs versus in situ bbp data for each wavelength and dataset considered and for the merged data set. The dashed line represents the 1:1 ratio. Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively. Table 4. Statistical descriptors of the difference between the QAA-derived bbp from satellite CCI Rrs and in situ bbp for each dataset with Raman scattering compensation. Figure A4 provides a graphical representation of this

Conclusions
The main findings of this work and their relevance for ocean color studies are summarized here: (1) Raman scattering compensation of Rrs prior to the application of the QAA significantly reduces errors in the retrieval of bbp with respect to in situ bbp. Inclusion of this processing step in operational schemes is recommended. (2) The QAA-derived bbp from in situ radiometry has negligible biases with respect to in situ bbp. Figure 6. Scatter plots of QAA derived b bp from CCI R rs versus in situ b bp data for each wavelength and dataset considered and for the merged data set. The dashed line represents the 1:1 ratio. Pink, yellow, and green dots refer to V19, BOU, and CNR data, respectively. Table 4. Statistical descriptors of the difference between the QAA-derived b bp from satellite CCI R rs and in situ b bp for each dataset with Raman scattering compensation. Figure A4 provides a graphical representation of this

Conclusions
The main findings of this work and their relevance for ocean color studies are summarized here: (1) Raman scattering compensation of R rs prior to the application of the QAA significantly reduces errors in the retrieval of b bp with respect to in situ b bp . Inclusion of this processing step in operational schemes is recommended. (2) The QAA-derived b bp from in situ radiometry has negligible biases with respect to in situ b bp .
(3) CCI R rs shows low biases but higher RMS differences with respect to in situ data, that could be excessive for the monitoring of natural change over short periods. Here, the standardization of in situ radiometry protocols is highly encouraged [52], in order to reduce the errors when in situ datasets formed by multiple contributors are merged and used for R rs matchup analysis. (4) In part as a consequence of the findings above, QAA-derived b bp from CCI R rs displays negligible biases respect to in situ b bp , with moderately low RMS errors. (5) The in situ radiometry-derived spectral backscattering slope (η) has low predictive value as compared to η derived from b bp matchups. In this context, the impact of using the best fitted curve instead of the widely used expression [22] is negligible, thus validating the application of the latter without its retuning.
Notwithstanding these results, one future challenge should be to evaluate the impact of two other sources of inelastic scattering before the application of QAA on R rs : (i) red fluorescence, caused by chlorophyll, that usually plays an important role around the peak close to 685 nm; and (ii) the blue fluorescence, caused by CDOM, that can be relevant close to the peak at 425 nm [53].
In addition, there is the need of increasing the amount of spatial and spectral coverage of high-quality in situ b bp observations. As of today, available multispectral b bp is limited to a small number of ship-borne data, or longer datasets but in fixed points (i.e., buoy). On the other hand, Biogeochemical-Argo floats cover large areas but their data are mainly given at a single band. Therefore, there is need to significantly increase the amount of b bp data at multiple bands, seasons and geographical regions. New technological developments on autonomous platforms will aid to enhance data density across many water types, to extend the CCI uncertainty derivation approach to b bp as well, thus allowing the mapping of uncertainties for every b bp product.
Lastly, in situ b bp measurements lag behind the standards on protocols and uncertainty characterization with respect to other quantities such as the radiometry [52]. Only when in situ uncertainty-characterized datasets, from instrument characterization to deployment [54], become available, more detailed algorithm validation could be performed and this will help to better evaluate the influence of optically active constituents (e.g., CDOM, chlorophyll). Acknowledgments: This is a contribution to the ESA Ocean Colour Climate Change Initiative of the European Space Agency. ESA and CNES are thanked for funding the BOUSSOLE project. We are grateful to all the contributors to the V19 dataset. J.P. thanks CNR-ISMAR for the stay at the Rome CNR location in which this manuscript was elaborated. Finally, we wish to thank four anonymous reviewers for their criticisms and suggestions that helped the manuscript to be improved.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
This appendix includes the graphical representation of the information in Tables 1-4, for a quicker interpretation of the results and the derived conclusions. The error bars are made by taking the mean bias in tables as central values and the standard deviation (σ) as bar width, calculated as σ = √ RMS 2 − Bias 2 . When error bars intersect the zero-difference line, the differences are assumed to be not significant. Figure A1. Relative differences between the QAA-derived bbp and in situ bbp for each dataset, without Raman scattering compensation. Data taken from Table 1. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.  Table 2. Pink, yellow and green lines refer to V19, BOU, and CNR data, respectively.  Table 1. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively. Figure A1. Relative differences between the QAA-derived bbp and in situ bbp for each dataset, without Raman scattering compensation. Data taken from Table 1. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.  Table 2. Pink, yellow and green lines refer to V19, BOU, and CNR data, respectively.  Table 2. Pink, yellow and green lines refer to V19, BOU, and CNR data, respectively. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 16 Figure A3. Relative differences between satellite CCI Rrs and in situ Rrs for each dataset. Data taken from Table 3. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively. Figure A4. Relative difference between the QAA-derived bbp from satellite CCI Rrs and in situ bbp for each dataset with Raman scattering compensation. Data taken from Table 1. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively. Figure A3. Relative differences between satellite CCI R rs and in situ R rs for each dataset. Data taken from Table 3. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 16 Figure A3. Relative differences between satellite CCI Rrs and in situ Rrs for each dataset. Data taken from Table 3. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively. Figure A4. Relative difference between the QAA-derived bbp from satellite CCI Rrs and in situ bbp for each dataset with Raman scattering compensation. Data taken from Table 1. Pink, yellow, and green lines refer to V19, BOU, and CNR data, respectively.