Improving the Retrieval of Carbon-Based Phytoplankton Biomass from Satellite Ocean Colour Observations

: Phytoplankton is at the base of the marine food web and plays a fundamental role in the global carbon cycle. Ongoing climate change signiﬁcantly impacts phytoplankton distribution in the ocean. Monitoring phytoplankton is crucial for a full understanding of changes in the marine ecosystem. To observe phytoplankton from space, chlorophyll- a concentration (Chl) has been widely used as a proxy of algal biomass, although it can be impacted by physiology. Therefore, there has been an increasing focus towards estimating phytoplankton biomass in units of carbon (C phyto ). Here, we developed an algorithm to quantify C phyto from space-based observations that accounts for the spatio-temporal variations of the backscattering coe ﬃ cient associated with the fraction of detrital particles that do not covary with Chl. The main ﬁndings are: (i) a spatial and temporal variation of the detritus component must be accounted for in the C phyto algorithm; (ii) the reﬁned C phyto algorithm performs better (relative bias of 23.7%) than any previously existing model; and (iii) our algorithm shows the lowest error in C phyto across areas where picophytoplankton dominates (relative bias of 14%). In other areas, it is currently not possible to accurately assess the performance of the reﬁned algorithm due to the paucity of in situ carbon data associated with nano- and micro-phytoplankton size classes.


Introduction
Phytoplankton is responsible for approximately half of global primary production and is at the base of the marine food web [1]. Phytoplankton is consequently a fundamental actor in the global carbon cycle [2]. Moreover, these organisms are regarded as sentinels of changes in the ocean because of their capacity to respond rapidly to environmental perturbations. Several factors, such as ocean circulation, anthropogenic activities, and climate affect phytoplankton abundance and distribution. In particular, phytoplankton spatio-temporal patterns are expected to vary with climate change [3]. Monitoring the distribution of phytoplankton is crucial for a full understanding of the oceanic biogeochemical cycles.
Ocean colour observations have significantly improved our capability to map phytoplankton global distribution since the 1970s. Chlorophyll-a concentration (Chl; in mg m −3 ) is a consolidated proxy of algal biomass [4]. However, Chl variability may not always correspond to actual changes in algal biomass but rather to cellular physiological adjustments in response to environmental stressors such as light and nutrient limitation [5,6]. Therefore, there has been an increasing focus on the estimation of phytoplankton biomass in units of carbon (C phyto ; mg C m −3 ). C phyto has found application in studies of: (i) primary production [7,8]; (ii) phytoplankton physiology [9]; (iii) phytoplankton growth/loss rates [10,11]; (iv) comparison with marine ecosystem models [12]; (v) pools of carbon in the ocean (i.e., their stocks and turnover rates); and (vi) ocean-atmosphere interactions in marine phytoplankton feedback in atmospheric aerosols [13].
Martinez-Vicente et al. (2017) recently reviewed algorithms for deriving C phyto from satellite ocean colour observations [14]. Three groups of algorithms are currently available: (i) backscatter-based models developed using satellite [7] or in situ data [14][15][16]; (ii) Chl-based models [17,18]; and (iii) allometric conversion-based models [19][20][21]. Among these algorithms, those relying on the particulate optical backscattering coefficient, b bp (in m −1 ), raised high interest because b bp is related to the concentration and composition of suspended particles in seawater, and to their size and shape [22]. Although former research suggested that b bp was mainly influenced by submicron detrital particles [22], it has recently been demonstrated that most of the b bp signal in the surface oligotrophic ocean is due to particles with equivalent spherical diameters between 1 and 10 µm [23], thus reinforcing the role of phytoplankton as a source of the open-ocean b bp [24], and thus the usefulness of b bp to retrieve C phyto . However, satellite b bp -based algorithms show large errors in predicting C phyto mainly because the current formulation does not take into account the spatio-temporal changes of non-algal particles (NAP) [25].
Here we refine the first and widely applied satellite b bp -based C phyto model [7] and make crucial arrangements to improve its performance. Beh05 uses the linear relationship between Chl and b bp to estimate the background contribution of NAP to b bp (443), b k bp , which corresponds to b bp when Chl is zero (i.e., the intercept of the linear fit). Beh05 is based on the following equation: In Beh05, SF is the scaling factor used to convert b bp into C phyto and is set to 13,000 mg C m −2 [7]. b k bp is assumed to represent a stable surface population of heterotrophic and detrital particles that is not strictly dependent on phytoplankton dynamics and physiology. Thus, in Beh05, b k bp was assumed as a constant value (3.5·10 −4 m −1 ). However, NAP can include various types of particles ranging from heterotrophic organisms such as bacteria, micrograzers, and viruses, to fecal pellets and cell debris, and mineral particles and bubbles [22] that vary both in space and time. The nature, definition, and evaluation of b k bp is still under debate and has been investigated using in situ [24][25][26] and satellite data [27,28]. Understanding the spatial and temporal dynamics of b k bp in the open ocean will improve estimations of C phyto and, more generally, our knowledge on the fate of marine carbon.
Beh05 estimated b k bp as the offset (intercept) of the least-square regression between monthly satellite-derived b bp and Chl computed via the Garver-Siegel-Maritorena (GSM) semi-analytical algorithm [29]. The regression analysis used only Chl values higher than 0.14 mg m −3 to separate the changes in Chl due to physiology (i.e., photoacclimation) from actual changes in biomass, and thus in C phyto . Bellacicco et al. (2018) used satellite data from the Mediterranean Sea to estimate b k bp for various bioregions and seasons, and showed that b k bp varied regionally and seasonally, thus rejecting the hypothesis of invariance of the heterotrophic and detrital components at the sea surface [27]. Global variability in b k bp was observed [28], with a reported median global value of 9.5 10 −4 m −1 , nearly threefold that found by Beh05. The newly evaluated b k bp resulted in C phyto that were half those estimated by Beh05.
The Chl-b bp relationship is influenced by phytoplankton specific composition and diversity (e.g., size, shape, internal structure), physiology (e.g., photoacclimation), and the nature of NAP itself [22,30,31]. For these reasons, Brewin et al. (2012) (hereafter Bre12; [26]) developed an analytical non-linear fit between b bp and Chl that accounted for changes in phytoplankton cell size. In Bre12, the offset of the non-linear fit (i.e., b k bp ) between b bp and Chl in clear waters was interpreted as a constant background of NAP or also partly influenced by very small phytoplankton (e.g., prochlorophytes) in addition to NAP. In Bre12, the reported b k bp values were 7.0·10 −4 m −1 at the wavelength of 470 nm and 5.6·10 −4 m −1 at 526 nm. More recently, the Bre12 model was used with a large database of 0-1000 m depth profiles of Chl and b bp acquired by the global Biogeochemical-Argo (aka BGC-Argo) float array. b k bp was demonstrated also to vary over depth and to account only for a small fraction of total b bp in productive areas, while being the main source of b bp in oligotrophic waters [25].
The former evidence suggests that if a realistic spatio-temporal variation of b k bp is introduced in the C phyto algorithm its precision and accuracy can be largely improved. Specifically, the hypothesis that b k bp varies in space and time is here applied monthly at the scale of satellite pixel rather than using a unique b k bp value or pre-defined bioregions. Algorithm outputs are validated against the only available, to the best of our knowledge, in situ dataset and compared to the performances of Beh05, Bel18, and Bre12 with single b k bp values. The performance of the new algorithm is also compared to empirical approaches [14,16]. The analysis was developed globally by partitioning the data among Optical Water Class (OWC; [32]) to explore the efficiency and applicability of the different tested approaches.

Ocean Colour Data
The full European Space Agency (ESA) Ocean Colour-Climate Change Initiative (OC-CCI) time series (1997-2019) of global daily R rs and Chl version 4.2 data at 4 km resolution was downloaded from the ESA OC-CCI FTP server (https://esa-oceancolour-cci.org; ftp://oceancolour.org/occci-v4.2/ geographic/netcdf/daily/rrs/). ESA OC-CCI data products are the result of the merging of SeaWiFS, MERIS, MODIS, and VIIRS observations in which the inter-sensor biases are removed. Version 4.2 includes the latest NASA reprocessing R2018.0, which mostly accounts for the ageing of the MODIS sensor. The ESA OC-CCI provides the daily R rs data and associated uncertainty maps in terms of bias and root mean square error. In this work, for each daily R rs , the bias was also extracted and then corrected pixel-by-pixel, as recommended in the Product User Guide [33].
Chl was estimated with a blending of the OCI (as implemented by NASA, itself a combination of CI and OC4), OC5, and OC3 algorithms. For each daily Chl value, the bias was also extracted and compensated at the pixel scale. Daily b bp maps at 443 nm were produced from daily R rs for the same period (1997-2019) by applying the Quasi-Analytical Algorithm (QAA; [34,35]) with prior correction for Raman scattering [27,36]. The accuracy of the QAA in retrieving b bp was fully demonstrated [27,[36][37][38][39]. In this work, the QAA algorithm was thus selected for its high efficiency in the b bp retrievals, and to show consistency and coherency with the OC-CCI program in which the QAA is the designated algorithm [37,40].
The daily datasets were then under-sampled to 25 km resolution to resolve only the broader oceanographic scales of variability.
2.2. Computation of b k bp b k bp was estimated for every pixel and month by using the following basic relationship [7,25,27,28]: where k is the slope of the least-square regression fit between daily Chl and b bp , and b k bp is the intercept of the fit [7]. The term k·Chl thus identifies the fraction of b bp that covaries with Chl [7,25,27,28]. A scheme of the algorithm implementation is shown in Figure S1.
We computed monthly climatological b k bp maps at 25 km by applying Equation (2) to daily Chl and b bp data within each month of the year from 1998 to 2019. The estimation of b k bp , and consequently of C phyto , relies on good relationships between Chl and b bp which then constitutes the fundamental condition to exploit this method. To this aim, for each b k bp map, the significance S (using Student's t-test), the Pearson correlation coefficient (r; not shown), and the 1σ uncertainty maps (due to the linear fit between Chl and b bp ; see Section 15.2 of [41]) were computed to give an estimate of the robustness of the fit.
We then applied for each pixel of every resulting monthly b k bp map an additional moving average of 1000 km to remove any source of noise and smaller scale variability, while conserving the large-scale oceanic patterns [42].
Finally, the monthly b k bp maps computed with the new approach were interpolated to obtain daily climatological b k bp maps at 25 km resolution. Daily pixel-based estimates were then used to assess C phyto by applying a revised Equation (1) [7,27,28], such as: In some pixels, C phyto may show values less of, or close to 0. This occurred for pixels where the Chl-b bp relationship had a significance S < 0.95 and r ≤ 0. In these pixels, b k bp may be higher than b bp thus giving non-reliable, negative estimations of C phyto (less than 0; e.g., in the subtropical gyres). For these areas, we applied a threshold of C phyto equal to 0.13 mg C m −3 to highlight that C phyto is expected to be low [14]. All match-up points falling in the areas where the Chl-b bp relationship showed a significance S < 0.95 and r ≤ 0 were not included in the validation and algorithm performance analysis ( Figure S1). These areas thus need to be interpreted and managed with caution.
For comparative purposes, C phyto was also computed following the approaches of Beh05, Bel18, and Bre12 with single

In Situ Reference Data
The in situ C phyto database [14] is a compilation of data acquired from many sources (e.g., MAREDAT and Atlantic Meridional Transect, AMT, cruises) for a total of N = 557 data points and consists of carbon biomass of picophytoplankton organisms (i.e., cell size < 2 µm). It was downloaded from http://www.zenodo.org (doi:10.5281/zenodo.1067229). For complete details about the in situ dataset see [14]. Only pixels with a good (S ≥ 0.95 and r > 0) satellite relationship between Chl and b bp were retained so that the original 557 data points decreased to a total of 396 matchups ( Figure 1). The final matchup database encompassed from oligotrophic to mesotrophic waters (Chl from 0.035 to 3.13 mg m −3 and C phyto from 1.80 to 60.25 mg C m −3 ) and OWC from 1 to 13. OWC from 1 to 6, corresponding to less productive waters, represented 56% of the in situ data.

Statistical Assessment
Estimated satellite C phyto (y i ) computed with Equation (3) was compared to reference in situ data (x i ) by computing the bias (δ; mg m −3 ), the relative bias in percentage (∇; in %), and the standard deviation of the difference (σ ∆ ; mg m −3 ). The assessment was made for each OWC individually and for all matchups at once:

Algorithm Performance for Cphyto Retrievals
Results and statistics of the comparison between satellite Cphyto estimates obtained with our algorithm and the previously published values, and in situ Cphyto, are reported in Table 1 and Figure  2. Overall, algorithms perform with  values spanning from 23.7% to 203%. The algorithm here developed shows the smallest bias ( of 23.7%). By comparison, the algorithms based on a constant bp k value have lower performance, and Gra15, Beh05, and MV17 showed a systematic overestimation for low Cphyto values (Table 1 and Figure 2).
The bbp-based algorithms that use a single bp k constant value (Beh05, Bre12, and Bel18) show similar  ∆ . Specifically, Bre12 has a bp k value that is twice as high as that used in Beh05, and Bel18 has a value nearly three times larger than that of Beh05.  ranges from 23.7% (Bel18) and 80.1% (Bre12), to 130.3% (Beh05), resulting in different performance in relation to the selection of the single value.
In the case of Beh05, a source of discrepancy leading to such a bias may be linked to the bbp input from the algorithm as also highlighted in [27]. Indeed, the original equation was derived using the GSM algorithm [29] for obtaining a relationship between Chl and bbp and thus for the bp k estimation, instead of the QAA, applied here. In addition, it has recently been found that Raman scattering (here corrected prior to the QAA application) plays a fundamental role for the retrievals of bbp in very clear waters because, if not corrected, it results in overestimation of bbp [36,43]. However, these results highlight how the use of a single bp k determines a strong overestimation of satellite Cphyto, in line with the results of MV17 (Table 1 and Figure 2). The comparison between our

Algorithm Performance for C phyto Retrievals
Results and statistics of the comparison between satellite C phyto estimates obtained with our algorithm and the previously published values, and in situ C phyto , are reported in Table 1 and Figure 2. Overall, algorithms perform with ∇ values spanning from 23.7% to 203%. The algorithm here developed shows the smallest bias (∇ of 23.7%). By comparison, the algorithms based on a constant b k bp value have lower performance, and Gra15, Beh05, and MV17 showed a systematic overestimation for low C phyto values (Table 1 and Figure 2).
The b bp -based algorithms that use a single b k bp constant value (Beh05, Bre12, and Bel18) show similar σ ∆ . Specifically, Bre12 has a b k bp value that is twice as high as that used in Beh05, and Bel18 has a value nearly three times larger than that of Beh05. ∇ ranges from 23.7% (Bel18) and 80.1% (Bre12), to 130.3% (Beh05), resulting in different performance in relation to the selection of the single value. Table 1. Basic statistics for each single approach. For each approach is reported the bias (δ; in mg m −3 ), the standard deviation (σ ∆ ; in mg m −3 ), and the relative percentage bias (∇; in %). Note that OWCs 1 and 2 are grouped in a single class to increase the number of observations available for the statistics. This is also replicated for OWCs 11, 12, and 13. The statistics using log 10 transformed data for all OWCs are reported in Table S1.   Table 1.

Spatial and Temporal Distribution of Cphyto
The global mean climatology of Cphyto ( Figure 3a) agrees with the expected geographical distribution [14]. Cphyto highs are found in the most productive regions such as the high latitude regions and coastal areas, and in upwelling systems such as those off Mauritania and Perù, the Benguela and Kuroshio currents, and along the equatorial belt of the Pacific Ocean (Figure 3a). Cphyto lows are found in mid-latitude areas such as the subtropical gyres or the Eastern Mediterranean Sea and Black Sea. Cphyto varies within the same order of magnitude of recent estimates made using satellite data [14] and biogeochemical models [44]. The two main oceanic regimes as defined by [45] are also highlighted from our analysis: "photoacclimation-dominance" and "biomass-dominance". The former is typical of oligotrophic areas (e.g., subtropical gyres) where the variability of Chl is uncoupled from biomass and the process of acclimation to light and nutrients drives Chl variations [45]. In these areas, Cphyto shows the lowest mean values (<10 mg C m −3 ). Conversely, the latter is typical of most productive areas where Chl and bbp strongly covary [31,46]. The high Chl-bbp covariability indicates that particles (and biomass) co-vary with phytoplankton abundance, whereas the physiological photoacclimation process plays a secondary role in determining the Chl variations; e.g., in the North Atlantic and Southern Oceans, Cphyto shows the highest mean values (>30 mg C m −3 ).

Figure 2.
Matchups between satellite C phyto versus in situ C phyto data for each algorithm (N = 396). The dashed line represents the 1:1 ratio. Gold dots represent in situ C phyto data corresponding to OWCs 1-6; blue dots represent in situ C phyto data for OWCs 7-13. Statistics are reported in Table 1.
In the case of Beh05, a source of discrepancy leading to such a bias may be linked to the b bp input from the algorithm as also highlighted in [27]. Indeed, the original equation was derived using the GSM algorithm [29] for obtaining a relationship between Chl and b bp and thus for the b k bp estimation, instead of the QAA, applied here. In addition, it has recently been found that Raman scattering (here corrected prior to the QAA application) plays a fundamental role for the retrievals of b bp in very clear waters because, if not corrected, it results in overestimation of b bp [36,43].
However, these results highlight how the use of a single b k bp determines a strong overestimation of satellite C phyto , in line with the results of MV17 (Table 1 and Figure 2). The comparison between our approach and previously published models underlines the necessity to account for the spatio-temporal variability of the b k bp for satellite-based C phyto calculation. However, because Bre12 and Bel18 show relatively low bias, such approaches might be used alternatively to the new approach proposed here because it may occur in the most oligotrophic waters (e.g., subtropical gyre areas). Note that in Bre12, the b k bp was at 470 nm so that the overestimation of C phyto can decrease if the b k bp is reported to the band of 443 nm.
The empirical algorithms of Gra15 and MV17 yield the worst results in terms of ∇. In the case of Gra15, this is the consequence of the few in situ measurements used for the algorithm definition, located in the Equatorial Pacific Ocean and in the Atlantic Ocean. Extrapolation of their algorithm, developed at 470 nm, to 443 nm, may also have had an effect. Regarding MV17, the results found are consistent with those reported in Tables 3 and 5 of their work.
MV17 reported that the geographical distribution of the in situ data, the median C phyto concentrations, and the carbon-to-chlorophyll ratios are in agreement with previous observations in oligotrophic oceanic conditions [16,18]. The fact that the in situ C phyto dataset used here contains only data from the picophytoplankton population must be taken into account for the interpretation of the algorithm evaluation. In most oligotrophic waters (OWCs from 1 to 6), δ shows the lowest values for all of the algorithms, with a maximum value of 14.4 mg m −3 (Table 1). However, when analyzing the ∇ index, the C phyto computed by considering the spatio-temporal variability of b k bp is the most efficient (Table 1), with the exception of OWCs 5 and 6, for which Bel18 shows the best performance. Indeed, within OWCs in which picophytoplankton is expected to dominate (OWCs from 1 to 6), the performance of the majority of the tested algorithms generally improves (Table 1), with the exceptions of Gra15 and MV17. As expected, there is a decrease in the accuracy for all of the algorithms in the OWCs from 7 to 13. This could be influenced by the dominance of nano-and micro-phytoplankton (cell size > 2 µm) in these waters, whereas the picophytoplankton contribute in low relative proportions [14]. Because most of the b bp signal is due to particles with equivalent diameters larger than that those associated with picophytoplankton [23], a decreasing performance of the b bp -based algorithms in water dominated by larger cells (e.g., nano-and micro-phytoplankton) was expected (Table 1).

Spatial and Temporal Distribution of C phyto
The global mean climatology of C phyto (Figure 3a) agrees with the expected geographical distribution [14]. C phyto highs are found in the most productive regions such as the high latitude regions and coastal areas, and in upwelling systems such as those off Mauritania and Perù, the Benguela and Kuroshio currents, and along the equatorial belt of the Pacific Ocean (Figure 3a). C phyto lows are found in mid-latitude areas such as the subtropical gyres or the Eastern Mediterranean Sea and Black Sea. C phyto varies within the same order of magnitude of recent estimates made using satellite data [14] and biogeochemical models [44]. The two main oceanic regimes as defined by [45] are also highlighted from our analysis: "photoacclimation-dominance" and "biomass-dominance". The former is typical of oligotrophic areas (e.g., subtropical gyres) where the variability of Chl is uncoupled from biomass and the process of acclimation to light and nutrients drives Chl variations [45]. In these areas, C phyto shows the lowest mean values (<10 mg C m −3 ). Conversely, the latter is typical of most productive areas where Chl and b bp strongly covary [31,46]. The high Chl-b bp co-variability indicates that particles (and biomass) co-vary with phytoplankton abundance, whereas the physiological photoacclimation process plays a secondary role in determining the Chl variations; e.g., in the North Atlantic and Southern Oceans, C phyto shows the highest mean values (>30 mg C m −3 ).
The amplitude of seasonal variations (i.e., interannual variability) of C phyto is shown in Figure 3b. Most of the variability is observed at high latitudes (e.g., Baltic Sea, Irish Sea, Norwegian Sea, North Sea), in the North Pacific and North Atlantic Oceans, and in the Southern Ocean, because σ std values are higher than 12 mg m −3 . The North Atlantic and Southern Oceans are characterized by high seasonal variations due to intense productive months followed by unproductive periods [47][48][49] (Figures S2-S13). In the North Atlantic Ocean, C phyto starts to increase in May, reaching the maximum values in June and July (>50 mg C m −3 ) following the spring bloom. Then, C phyto decreases to the lowest values for the remainder of the year [47]. In the Southern Ocean, C phyto follows the typical seasonal cycle with high values during the most productive months which occur from November to March (>50 mg C m −3 ). Highs are mostly observed in the Antarctic Circumpolar Current (ACC) and the Patagonia Shelf. Later, there is a decrease in C phyto from May to September (Figures S2-S13). Additionally, the northwestern Indian Ocean, known to be an Oxygen Minimum Zone, as well as the Benguela upwelling system, also show significant variability in C phyto , with σ std greater than 12 mg m −3 . Conversely, mid-latitude areas (e.g., in the subtropical gyres), display much lower variability, with σ std less than 2 mg C m −3 . Figure S14 clearly shows how the single b k bp of Beh05 is located at the lowest limit of possible values, confirming that reported in [28], whereas the Bel18 and Bre12 single values are representative of a wider area. In the subtropical gyres, b k bp does not significantly deviate from the coefficient found by Beh05, but in the most productive regions (e.g., North Atlantic and Southern Oceans, high latitude seas, and coastal upwelling regions), they largely differ with great impact on C phyto estimations. This discrepancy in C phyto is still valid in the case of the use of Bel18 and Bre12 single values, but for other areas (e.g., in the subtropical gyres or in the Mediterranean Sea), the use of the single value of Bel18 or Bre12 yields an underestimation of C phyto , whereas in the productive Arctic and Southern Oceans C phyto is generally overestimated. This is also indicated in the seasonal b k bp variation as reported in the Figures S2-S13. In the North Atlantic Ocean, b k bp shows high values during the spring bloom and low values from December to April. However, in this area, the b k bp is only a small fraction of the entire b bp signal as the consequence of the high variability in the b bp and phytoplankton abundance [25,27] thus following the phytoplankton dynamics. In the Southern Ocean, and especially the region influenced by the ACC and Patagonia Shelf, b k bp shows highs from December to April (i.e., due to coccoliths) and lows in the period May-September in accordance also with the C phyto seasonal cycle [25]. In less productive seas, the seasonal cycle of b k bp is smooth, suggesting low NAP variations. In these areas, the seasonal b k bp change is weak throughout the year but also dominates the b bp signal [25]. These areas are characterized by limited nutrient availability determining low C phyto , mostly of picophytoplankton class, high bacterial, and detrital concentrations [50][51][52]. It follows that the use of a varying b k bp enables accounting for its seasonal variations in the C phyto computation, giving the model more flexibility in relation to the change of the biogeochemical and trophic conditions.  The amplitude of seasonal variations (i.e., interannual variability) of Cphyto is shown in Figure 3b. Most of the variability is observed at high latitudes (e.g., Baltic Sea, Irish Sea, Norwegian Sea, North Sea), in the North Pacific and North Atlantic Oceans, and in the Southern Ocean, because std values are higher than 12 mg m −3 . The North Atlantic and Southern Oceans are characterized by high seasonal variations due to intense productive months followed by unproductive periods [47][48][49] ( Figures S2-S13). In the North Atlantic Ocean, Cphyto starts to increase in May, reaching the maximum values in June and July (> 50 mg C m −3 ) following the spring bloom. Then, Cphyto decreases to the lowest values for the remainder of the year [47]. In the Southern Ocean, Cphyto follows the typical seasonal cycle with high values during the most productive months which occur from November to Annual global mean (a) and standard deviation (b) of C phyto . Note that at high latitudes the number of satellite observations used for the mean and standard deviation computations are lower than those used at low or mid-latitudes due to the winter darkness (i.e., high sun zenith angle).

Caveats of the C phyto Algorithm
Three main caveats of the algorithm should be noted. The first is that C phyto assumes a constant scaling factor, SF, equal to 13,000 mg C m −2 , following [7]. This can be an additional possible source of error in the satellite retrievals of C phyto . Therefore, one important future effort should be to investigate a refined scaling factor relating b bp to C phyto coupled with the b k bp space-time variability as presented in this work. With this perspective, additional laboratory work should be done to evaluate if change in SF values can affect the C phyto estimations.
The second caveat of the algorithm is that it relies on a tight relationship between b bp and Chl, which is influenced by the algorithms (including the atmospheric correction) used for Chl and b bp retrievals, in addition to the environmental conditions. In some areas of the subtropical gyres, the Chl-b bp relationship showed S < 0.95 (Figures S2-S13) and r smaller than or close to 0, as also reported by [28]. The main reason for this relationship may be the photoacclimation process, known to be dominant in such areas over the year, and which introduces variability in Chl that is unrelated to biomass. In these areas, our refined method must be applied with caution. With this perspective, one future challenge is to solve this limitation in b k bp computation and, consequently, in the satellite C phyto . Indeed, as the subtropical gyres cover about 60% of the global ocean, it can be of a great importance to compute C phyto and Chl:C phyto to understand oceanic biogeochemical cycles. Conversely, the remaining areas of global ocean (e.g., the North Atlantic Ocean) present significant positive r and highly significant S (Figures S2-S13) throughout the year, indicating that either b bp or Chl can be used for determining the phytoplankton dynamics and distribution, as expected in open-ocean waters [27]. In these areas, the algorithm performance is robust.
The third caveat is that the algorithm validation is restricted only to in situ C phyto data associated with picophytoplankton carbon. It determines that the algorithm performance has to be interpreted with caution in those areas where other phytoplankton size classes dominate. This means that, currently, it is not possible to define the accuracy for OWCs from 7 to 13. Indeed, Table 1 shows that there is a decrease in the accuracy for all of the algorithms in the OWCs where nano-and micro-phytoplankton dominates. Thus, one future need is to improve the in situ C phyto dataset with new measurements representative of all of the phytoplankton size classes.

Conclusions and Future Perspectives
In this work, a revisited version of the original b bp -based algorithm for C phyto retrieval from space [7] was proposed. The spatio-temporal variations of the background backscattering coefficient of non-algal particles (b k bp ) is the crucial part of this refinement, and builds upon a series of recent studies [25,27,28]. The main findings are: 1.
The new C phyto algorithm proposed here performs better than any previously published model, with a relative error of 24% (Table 1 and Figure 2) with respect to a reference in situ dataset.

2.
The new algorithm shows the lowest error in C phyto (14.0%) across most of the OWCs in which the picophytoplankton population dominates (Table 1). On the contrary, the highest errors (36.6%) occur in OWCs 7-13 in which larger phytoplankton cells are supposed to dominate (Table 1).
However, we acknowledge that the refined algorithm presented here could be improved by addressing the caveats mentioned above. Enlarging the databases of quality-controlled and freely accessible in situ C phyto at the different size-classes (pico-, nano-, and micro-), measured simultaneously with an array of optical variables (particularly b bp ) could help to reduce uncertainties in C phyto retrieval from space. Improving the accuracy of satellite C phyto could be of great importance to adding new knowledge regarding the contribution of phytoplankton to particulate organic carbon, and for the validation of ocean primary productivity and biogeochemical models at a wider scale [8,12].
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/21/3640/s1, Figure S1: Scheme of the algorithm developed, Table S1: Basic statistics of the validation analysis for all the matchup with log 10 transformed data: Figures S2-S13: global monthly climatologies of b k bp (a), 1-σ uncertainty (b), significance S (c) and C phyto (d); Figure S14: global mean and standard deviation b k bp maps.