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

**Jaime Pitarch 1,2, Marco Bellacicco 3,\*, Emanuele Organelli 2,4, Gianluca Volpe 2, Simone Colella 2, Vincenzo Vellucci <sup>5</sup> and Salvatore Marullo 2,3**


Received: 29 November 2019; Accepted: 20 December 2019; Published: 24 December 2019

**Abstract:** Particulate optical backscattering (bbp) is a crucial parameter for the study of ocean biology and oceanic carbon estimations. In this work, bbp retrieval, by the quasi-analytical algorithm (QAA), is assessed using a large in situ database of matched bbp and remote-sensing reflectance (Rrs). The QAA is also applied to satellite Rrs (ESA OC-CCI project) as well, after their validation against in situ Rrs. Additionally, the effect of Raman Scattering on QAA retrievals is studied. Results show negligible biases above random noise when QAA-derived bbp is compared to in situ bbp. In addition, Rrs 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 efficiency of QAA in the bbp 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.

**Keywords:** particulate optical backscattering; Raman scattering; QAA algorithm; ESA OC-CCI

#### **1. 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 (Rrs; units of sr<sup>−</sup>1) by using appropriate algorithms [3–5]. Among IOPs, the particulate optical backscattering coefficient (bbp; in m−1) is related to the particle concentration in seawater, on their size distribution, refractive index, shape and structure [6–8]. Former research suggested that bbp is mostly influenced by submicron non-algal particles [9–11]. However, it has been recently shown that most of bbp 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 bbp for the retrieval of: (i) particulate organic concentration (POC) [12,13]; (ii) particle size distribution [14,15]; and (iii) phytoplankton carbon biomass concentration (Cphyto; mg m−3) [16–18], a key parameter also for phytoplankton physiology studies [2,19,20]. Efficiency in the bbp retrieval is crucial for ocean biology and global ocean carbon estimations.

On one hand, radiative transfer theory provides the link between bbp and optical radiometry [21]. Therefore, inversion algorithms for bbp detection from optical radiometry can be developed. In particular, the quasi-analytical algorithm [3,22] is a multi-level algorithm that concatenates a sequence of empirical, analytical, and semi-analytical steps to retrieve spectral total non-water light absorption and backscattering (anw and bbp) first and to decompose anw into its CDOM, algal and non-algal contributions. Specifically, about bbp, some studies suggested some degree of bbp overestimation by the QAA [23,24], but their reference bbp data were sub-products of chlorophyll-a (Chl) measurements. QAA estimations from satellite Rrs showed a bias of +16.4% with respect to in situ bbp for the Adriatic Sea [25]. Using the in situ NOMAD dataset [26], a bbp overestimation of +38% by the QAA with respect to the observed value was reported [27]. Other results, based on in situ matchups, showed a bias of +2.5–8.8% for the QAA-derived bbp in Arctic waters, and of +9.5% to +16.4% in low-latitude waters [28]. Pitarch et al. [29] reported a slight underestimation within 10% in the Mediterranean Sea. Most recently, QAA-derived bbp from different satellite sensors (i.e., MODIS, VIIRS, OLCI) showed good performance with respect to a large in situ bbp dataset collected on biogeochemical (BGC)-Argo floats [30].

Currently, in the European Space Agency (ESA) Ocean Colour (OC) Climate Change Initiative (CCI), QAA is the selected algorithm to retrieve bbp. Specifically, the ESA OC-CCI project aims at creating a long-term, consistent, uncertainty-characterized time series of ocean color products, for use in climate-change studies [5,31]. In such a context, while in the case of Chl the uncertainties are fully provided, the bbp satellite products lack such information that is also crucial for POC and Cphyto estimations [1,32]. This absence of statistical assessment is influenced by the paucity of a sufficient number of in situ observations for the determination of uncertainties.

Nowadays, the uncertainties associated to QAA-based bbp retrievals globally are not known. In order to provide a best-effort bbp uncertainty assessment, this work aims at evaluating the efficiency of QAA for global bbp retrievals by using a large database of corresponding in situ Rrs and bbp data (*N* = 2881). In details, we use the updated version of the recent in situ global bio-optical dataset [33] together with field measurements from the BOUSSOLE buoy [34] and two different oceanographic cruises in the Tyrrhenian and Adriatic Seas. Unlike previous studies [29], here, the QAA performance is considered at multiple bands that further allow the evaluation of the bbp spectral slope retrievals against in situ measurements. The goals of this paper are thus: (i) to define the accuracy of the QAA for bbp retrievals using only in situ Rrs data; (ii) to validate the CCI Rrs with in situ corresponding data; and (iii) to evaluate the performance of the QAA using satellite CCI Rrs as input data.

#### **2. Data and Methods**

#### *2.1. Assessment of the Quasi-Analytical Algorithm (QAA)*

The original algorithm [3] has undergone many updates and developments by several researchers. The QAA version here used is based on the algorithm for the CCI bands which is currently integrated in the SeaWiFS data analysis system (SeaDAS) [22].

The sub-surface Rrs (named *rrs*) is calculated as *rrs* = *Rrs*/(0.52 + 1.7*Rrs*) and modeled as a function of the IOPs as: *rrs* = *g*0*u* + *g*1*u*2, with *u* = *bb*/(*a* + *bb*), g0 = 0.089 and g1 = 0.1245. This approach provided good results in the Mediterranean Sea in case of oligotrophic and coastal waters [29,35].

The QAA uses an empirical inversion of Rrs to retrieve absorption and then it solves total backscattering (bb) analytically. bbp is calculated by subtraction of pure seawater backscattering (bbw) for an average temperature of 14 ◦C and an average salinity of 38 PSU [36]. bbp is first estimated at

a reference wavelength of λ = 555 nm and then the calculation is extended to other wavelengths by assuming a power law *bbp* = *bbp*(λ0)(λ/λ0) <sup>−</sup><sup>η</sup> for the bbp with a spectral slope.

$$\eta\_{\parallel} = \,\_1p\_1 \left[ 1 - p\_2 \exp\left( -p\_3 \frac{r\_{rs}(443)}{r\_{rs}(555)} \right) \right] \tag{1}$$

Equation (1) is widely used for QAA retrievals of bbp 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*<sup>1</sup> = 2.0, *p*<sup>2</sup> = 1.2 and *p*<sup>3</sup> = 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 Rrs modeling [39]. Therefore, a pertinent question is how much this phenomenon affects the semi-analytical bbp retrievals. With this aim, Lee et al. [40] developed empirical correction formulas to compensate the Raman scattering on the Rrs. Here, we assess the effects of this compensation on the difference between the in situ bbp and Rrs-derived bbp. The statistical assessment was also replicated in other different cases: (i) validation of satellite CCI Rrs against in situ Rrs; and (ii) bbp retrievals after the application of QAA to satellite Rrs (Raman corrected), were compared to in situ measurements at the different available wavelengths.

Estimated data *yi* are compared to reference data *xi* by using the following statistical indexes: relative bias (units of %), relative root-mean square error (RMS, units of %) and determination coefficient (*r*2)

$$bias = 100 \frac{1}{N} \sum\_{i=1}^{N} \frac{y\_i - x\_i}{x\_i} \tag{2a}$$

$$RMS = 100\sqrt{\frac{1}{N} \sum\_{i=1}^{N} \left(\frac{y\_i - x\_i}{x\_i}\right)^2} \tag{2b}$$

$$r^2 = \frac{\sum\_{i=1}^{N} (\mathbf{x}\_i - \overline{\mathbf{x}})(y\_i - \overline{y})}{\sqrt{\sum\_{i=1}^{N} (\mathbf{x}\_i - \overline{\mathbf{x}})^2} \sqrt{\sum\_{i=1}^{N} (y\_i - \overline{y})^2}} \tag{2c}$$

#### *2.2. In Situ Data*

The in situ database is composed of three distinct datasets containing multi-spectral Rrs and bbp: 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 Rrs 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 Rrs data. It takes into account the spectral shape of the absorption and scattering that contribute to Rrs 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 Rrs and bbp co-located measurements around the global ocean (Figure 1). As shown in Figure 2, the total Rrs and bbp spectra cover from oligotrophic open-ocean to more eutrophic coastal waters as the range of Rrs and bbp values vary between 0–0.02 sr−<sup>1</sup> and 10<sup>−</sup>4–10−<sup>1</sup> m−<sup>1</sup> respectively.

**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.

#### 2.2.1. 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 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 Rrs and bbp measurements at all CCI bands were available. Such condition determines a total of *N* = 319 matchups. Remaining minor bbp 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.

#### 2.2.2. 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 bbp 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 bbp was applied that required a spectral bbp slope, calculated from every pair of two consecutive bands, within given bounds (more than −1 and less than 6). Rrs data were derived with a set of Satlantic 200-series multispectral radiometers ([46] and references therein). The Rrs 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 Rrs at 443, 490, 555, and 670 nm, only Rrs 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 Rrs at 555 nm was available, it was directly sampled and the Rrs at 560 nm was ignored. If the Rrs at 560 nm was available when the Rrs at 555 nm was missing, the Rrs at 560 nm was band-shifted to 555 nm. Similarly, in the red region, between the Rrs at 665 nm and 670 nm, preference to Rrs at 670 nm was given. Rrs 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.

#### 2.2.3. 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 Rrs 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 bb 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 (Lu) sensor was mounted onto a free-falling T-shaped structure, with the multicast technique. Above-water downwelling irradiance (Es) data were collected by a reference sensor, mounted at the top of the ship's deck. Rrs was computed using the SERDA software developed at CNR [47]. All the Rrs data were band-shifted to the CCI bands for consistency with the satellite Rrs. Further details about this dataset are provided in Pitarch et al. [29].

#### *2.3. Satellite ESA OC-CCI Rrs Data*

The ESA OC-CCI version 4.0 global daily Rrs 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 Rrs 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 Rrs 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 Rrs, the bias was also extracted and then compensated pixel-by-pixel.

#### **3. Results and Discussion**

#### *3.1. QAA Performance for bbp Retrievals from In Situ Data*

QAA-retrieved bbp from in situ Rrs is here compared to in situ measured bbp. Comparisons are made at the native bands of each in situ bbp 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 bbp and in situ bbp for each dataset, without Raman scattering compensation. Figure A1 provides a graphical representation of this table.


To understand the importance of the Raman scattering correction in semi-analytical bbp retrievals, the analysis is repeated with corrected Rrs [40]. The application of the Raman scattering correction reduces both bias and RMS nearly for all the bbp 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 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 table.


**Figure 3.** Scatter plots of QAA-derived bbp vs. in situ bbp data for each wavelength and dataset considered and for the merged dataset. Raman correction is applied to Rrs. The dashed line represents the 1:1 ratio. Pink, yellow and green dots refer to V19, BOU and CNR data, respectively.

#### *3.2. 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 (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.

**Figure 4.** η calculation considering all the in situ data available: V19 (pink dots), BOU (yellow dots), and CNR (green dots). The solid curve is the best fit of Equation (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.

#### *3.3. 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 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.

**Figure 5.** Scatter plots of CCI Rrs versus in situ Rrs 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 Rrs and in situ Rrs for each dataset. Figure A3 provides a graphical representation of this table.

#### *3.4. 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 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 table.

#### **4. Conclusions**

The main findings of this work and their relevance for ocean color studies are summarized here:


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 Rrs matchup analysis.


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 Rrs: (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 bbp observations. As of today, available multispectral bbp 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 bbp 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 bbp as well, thus allowing the mapping of uncertainties for every bbp product.

Lastly, in situ bbp 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).

**Author Contributions:** Conceptualization, J.P., M.B., and S.M.; Methodology, J.P., M.B., and S.M.; Data curation, J.P., E.O., G.V., S.C., V.V.; Writing-review and editing, all authors contributed equally. All authors have read and agreed to the published version of the manuscript.

**Funding:** J.P. has received partial funding from the 'Coastal Ocean Darkening' project funded by the Ministry for Science and Culture of Lower Saxony, Germany (VWZN3175). M.B. has a postdoctoral fellowship by the European Space Agency (ESA). This work was supported by the ESA Living Planet Fellowship Project PHYSIOGLOB: Assessing the inter-annual physiological response of phytoplankton to global warming using long-term satellite observations, 2018–2020.

**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 <sup>σ</sup> <sup>=</sup> <sup>√</sup> *RMS*<sup>2</sup> − *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.

**Figure A2.** Relative differences between the bbp-QAA derived and in situ bbp for each dataset with Raman scattering compensation. Data taken from Table 2. Pink, yellow and green lines refer to V19, BOU, and CNR data, respectively.

**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.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
