Ocean color algorithm for the retrieval of the particle size distribution and carbon-based phytoplankton size classes using a two-component coated-sphere backscattering model

The particle size distribution (PSD) of suspended particles in near-surface seawater is a key property linking biogeochemical and ecosystem characteristics with optical properties that affect ocean color remote sensing. Phytoplankton size affects their physiological characteristics and ecosystem and biogeochemical roles, e.g., in the biological carbon pump, which has an important role in the global carbon cycle and thus climate. It is thus important to develop capabilities for measurement and predictive understanding of the structure and function of oceanic ecosystems, including the PSD, phytoplankton size classes (PSCs), and phytoplankton functional types (PFTs). Here, we present an ocean color satellite algorithm for the retrieval of the parameters of an assumed power-law PSD. The forward optical model considers two distinct particle populations: phytoplankton and non-algal particles (NAPs). Phytoplankton are modeled as coated spheres following the Equivalent Algal Populations (EAP) framework, and NAPs are modeled as homogeneous spheres. The forward model uses Mie and Aden–Kerker scattering computations, for homogeneous and coated spheres, respectively, to model the total particulate spectral backscattering coefficient as the sum of phytoplankton and NAP backscattering. The PSD retrieval is achieved via spectral angle mapping (SAM), which uses backscattering end-members created by the forward model. The PSD is used to retrieve size-partitioned absolute and fractional phytoplankton carbon concentrations (i.e., carbon-based PSCs), as well as particulate organic carbon (POC), using allometric coefficients. This model formulation also allows the estimation of chlorophyll a concentration via the retrieved PSD, as well as percent of backscattering due to NAPs vs. phytoplankton. The PSD algorithm is operationally applied to the merged Ocean Colour Climate Change Initiative (OC-CCI) v5.0 ocean color data set. Results of an initial validation effort are also presented using PSD, POC, and picophytoplankton carbon in situ measurements. Validation results indicate the need for an empirical tuning for the absolute phytoplankton carbon concentrations; however these results and comparison with other phytoplankton carbon algorithms are ambiguous as to the need for the tuning. The latter finding illustrates the continued need for high-quality, consistent, large global data sets of PSD, phytoplankton carbon, and related variables to facilitate future algorithm improvements.

1 Introduction

Oxygenic photosynthesis by marine phytoplankton is a critical planetary-scale process supplying solar energy to the biosphere by fixing inorganic carbon; it is responsible for roughly half of global annual net primary productivity (e.g., Field et al., 1998). Ocean ecosystems play a key role in Earth’s carbon cycle and climate by affecting atmospheric CO₂ via the biological carbon pump, which sequesters some of the fixed carbon to the deeper ocean for longer timescales (e.g., Eppley and Peterson, 1979; Chisholm, 2000; Henson et al., 2011; Boyd et al., 2019; Brewin et al., 2021). The biological pump is influenced by the structure and function of oceanic ecosystems (e.g., Falkowski et al., 1998; Siegel et al., 2014); therefore, mechanistic, predictive understanding of ocean ecosystems is of high priority to Earth systems and climate research (e.g., Buesseler and Boyd, 2009; Siegel et al., 2016). Satellite remote sensing of ocean color is a key tool for the global characterization of ocean ecology (e.g., Siegel et al., 2013). This has led to large efforts to elucidate biological pump mechanisms using multiple platforms, including satellites, e.g., the EXPORTS program (Siegel et al., 2016).

Phytoplankton cell size (diameters varying from ≈ 0.5 to >50 µm (e.g., Clavano et al., 2007) is a key trait that affects multiple phytoplankton characteristics (Marañón, 2015), as well as sinking rates (e.g.,  Falkowski et al., 1998; Burd and Jackson, 2009; Stemmann and Boss, 2012; Siegel et al., 2014). Phytoplankton size classes (PSCs) thus tend to closely correspond to phytoplankton functional types (PFTs; e.g., Le Quéré et al., 2005). Importantly, phytoplankton cells also affect the inherent optical properties (IOPs) (e.g., absorption and backscattering coefficients) of the water column in a size-dependent manner (e.g., Mobley et al., 2002; Morel and Bricaud, 1986; Stramski and Kiefer, 1991; Kostadinov et al., 2009). This is because particle size (relative to the incident light wavelength) is one of the governing variables affecting the magnitude and spectral shape of light scattering and absorption caused by a particle (e.g., Bohren and Huffman, 1998). Therefore the particle size distribution (PSD) of phytoplankton (and other suspended particles in seawater) is a key property affecting both optical properties and cellular physiological and biogeochemical properties; i.e., it is a fundamental property linking ocean color remote sensing and ecosystem/biogeochemical characteristics. The size distribution of particles suspended in near-surface ocean waters is often described as a power law, given in differential form as follows (e.g.,  Bader, 1970; Sheldon et al., 1972; Jonasz, 1983; Boss et al., 2001; Twardowski et al., 2001; Kostadinov et al., 2009; Roy et al., 2017):

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{T}}}{\mathrm{d}D}}=N\left(D\right)=N_{\mathrm{0}}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\xi }},\end{array}$$

where D is particle diameter; N (m⁻⁴) is the differential number concentration of particles per unit volume seawater and per bin width of particle diameter; N₀=N(2  µm) is the particle number concentration at a reference diameter, here D₀=2 µm; and ξ is the power-law slope of the PSD. Equation (1) has to be integrated over a given diameter range to get the total particle number concentration in that range, N_T (m⁻³).

Ocean color is quantified by the spectral shape and magnitude of the remote-sensing reflectance, R_rs(λ) (sr⁻¹; also denoted simply R_rs for brevity below), where λ is the wavelength of light in vacuo. The Kostadinov–Siegel–Maritorena 2009 (KSM09, Kostadinov et al., 2009) algorithm retrieves the parameters of an assumed power-law PSD (ξ and N₀ in Eq. 1) from ocean color remote-sensing observations using the spectral shape (Loisel et al., 2006) and magnitude of the particulate backscattering coefficient, b_bp(λ) (m⁻¹). b_bp(λ) can be retrieved using existing inherent optical property (IOP) inversion algorithms; KSM09 uses the Loisel and Stramski (2000) IOP inversion. Subsequently, the retrieved PSD parameters allow the quantification of absolute and fractional PSCs: picophytoplankton, nanophytoplankton, and microphytoplankton based on bio-volume (Kostadinov et al., 2010) or phytoplankton carbon (Kostadinov et al., 2016a) (henceforth TK16) via allometric relationships (Menden-Deuer and Lessard, 2000). Phytoplankton carbon (phyto C) is the key variable of interest for carbon cycle and climate studies and modeling, and TK16 (data set available: Kostadinov et al., 2016b) represents a relatively unique carbon-based approach among PSC/PFT algorithms (Mouw et al., 2017) as it is based on knowledge of the PSD and allometric relationships to get at size-partitioned phyto C. Roy et al. (2013) and Roy et al. (2017) retrieve phytoplankton-specific PSD and size-partitioned phyto C based on the phytoplankton absorption coefficient.

The KSM09 PSD algorithm (and the TK16 phyto C/PSC derived from it) is built on the assumption of a single population of particles (approximated by homogeneous spheres), representing backscattering due to the entire oceanic particle assemblage: phytoplankton cells and non-algal particles (NAPs). However, particle internal composition and shape influence its optical properties (e.g Quirantes and Bernard, 2004; Quirantes and Bernard, 2006). Recent results suggest that the structural complexity of oceanic particles enhances backscattering significantly and can explain the so-called “missing backscattering” in the ocean (Organelli et al., 2018), i.e., the lack of optical closure between theoretically modeled and measured b_bp. Coated spheres (i.e., spheres consisting of concentric layers/shells of different material properties) can be used to better represent phytoplankton cells and their internal heterogeneity and composition (e.g., Bernard et al., 2009; Robertson Lain and Bernard, 2018), and they have significantly enhanced backscattering compared to their homogeneous equivalents (Duforêt-Gaurier et al., 2018; Organelli et al., 2018).

Here, we introduce a major improvement of the PSD algorithm formulation. Two separate particle populations are modeled: living phytoplankton cells and NAPs. Phytoplankton cells are modeled as coated spheres, following the Equivalent Algal Populations (EAP) framework (Bernard et al., 2009; Robertson Lain et al., 2014; Robertson Lain and Bernard, 2018). EAP explicitly models intracellular chlorophyll concentration, Chl_i, as governing the imaginary index of refraction and thus allows for bulk chlorophyll concentration (Chl) to be computed from a specific PSD. The coated-sphere EAP-based model is useful to better represent phytoplankton cells specifically; however, not all backscattering particles are phytoplankton (Stramski et al., 2004), and in fact, sub-micron NAPs even smaller than the smallest autotroph (≈ 0.5 µm in diameter) are critical for determining the spectral shape of b_bp, which is key for PSD retrieval with KSM09 and the algorithm presented here. Particles other than and smaller than phytoplankton are likely to significantly contribute to backscattering (Stramski et al., 2004; Zhang et al., 2020) in spite of evidence that phytoplankton/larger particles contribute more than Mie theory predicts, based on homogeneous spheres (e.g., Dall'Olmo et al., 2009). Thus, a two-component particle model is used here, separately modeling NAPs as homogeneous spheres of wider size range than phytoplankton, so that bulk b_bp of oceanic waters can be modeled (e.g.,  Stramski et al., 2001; Moutier et al., 2016; Duforêt-Gaurier et al., 2018). NAPs are modeled as having generally organic detrital composition, but with some allowance for higher indices of refraction to account for minerogenic particle contributions. The PSD forward model can thus also produce a first-order estimate of particulate organic carbon (POC) and the percent contribution of phytoplankton and NAPs to b_bp.

Subsequent sections present details of the two-component, EAP-based forward IOP model, the inversion methodology developed for operational application of the PSD algorithm, and the use of the satellite-derived PSD to retrieve derived products (following the methods of TK16 with some modifications), namely absolute and fractional size-partitioned phytoplankton carbon (henceforth phyto C) (i.e., carbon-based PSCs), as well as Chl and POC estimates. The novel algorithm is applied operationally to monthly data from the multi-sensor merged Ocean Colour Climate Change Initiative (OC-CCI) v5.0 data set (Sathyendranath et al., 2019, 2021); examples are shown in the paper, and the entire data set is publicly available and linked below (see “Data availability”). We then present and discuss an initial effort of validation of the new PSD algorithm and derived products using global compilations of PSD, picophytoplankton carbon, and POC in situ data. A comparison with other existing methods to retrieve phyto C is presented. We also discuss algorithm uncertainties, assumptions, and limitations as well as future work directions.

2 Data and methods

2.1 Particle optical model input specification for phytoplankton and NAPs

The contributions of two separate particle populations to bulk backscattering are modeled using Mie theory (Mie, 1908) for homogeneous spherical particles and the Aden–Kerker (Aden and Kerker, 1951) method for coated spheres. Living phytoplankton cells are represented by the first particle population, and all other suspended particles of any origin (i.e., NAPs) are represented by the second population. Living phytoplankton cells are modeled as coated spheres using the Equivalent Algal Populations (EAP) framework (Bernard et al., 2009; Robertson Lain et al., 2014; Robertson Lain and Bernard, 2018) for determining optical model inputs, in particular the complex indices of refraction of the particle core and coat. NAPs are modeled as homogeneous spheres meant to represent organic detritus, but also allowing for their real index of refraction to vary over a wider range to take into account the contribution of mineral particles.

A characteristic of the PSD algorithm presented here is that it is mechanistic to the extent feasible, i.e., based on first principles and causality, even at the expense of increasing complexity. For example, as in EAP, the imaginary refractive index (RI) of the cell is a function of intracellular chlorophyll concentration, Chl_i. We vary some optical model inputs in a Monte Carlo simulation in order to assess uncertainty and base the PSD inversion on an ensemble of forward runs rather than a single set of inputs. Details of uncertainty estimation and propagation are given in Supplement Sect. S1. Details of how each input parameter for phytoplankton cells and for NAPs is specified, as well as the statistical distributions from which the Monte Carlo simulation instances were picked, are specified in Tables 1 and 2.

As in the EAP model, the chloroplast is represented by the particle coat. Its relative volume, V_s, is picked from a distribution as shown in Table 1. The chloroplast's imaginary refractive index (RI) (relative to seawater) at 675 nm, n^′(675), is then computed as follows (Morel and Bricaud, 1986; Bernard et al., 2009; Robertson Lain and Bernard, 2018):

$$\begin{array}{}\text{(2)}& n^{\prime }\left(\mathrm{675}\right)={\displaystyle \frac{{\mathrm{Chl}}^{*}\times {\mathrm{Chl}}_{\mathrm{i}}\times {\mathrm{10}}^{\mathrm{6}}\times \mathrm{675}\times {\mathrm{10}}^{-\mathrm{9}}}{\mathrm{4}\mathit{\pi }\times V_{\mathrm{s}}\times n_{\mathrm{sw}}\left(\mathrm{675}\right)}},\end{array}$$

where Chl${}^{*}=\mathrm{0.027}$ m² mg⁻¹ is the theoretical maximum specific absorption coefficient of chlorophyll at 675 nm when dissolved in water (Bernard et al., 2009; Robertson Lain and Bernard, 2018), Chl_i is the intracellular chlorophyll concentration in kilograms of chlorophyll per cubic meter of cellular material, and n_sw(675) is seawater's absolute real RI at 675 nm. A hyperspectral basis vector from the EAP model (based on measurements; for details see Bernard et al., 2009; Robertson Lain and Bernard, 2018) is then scaled using the value at 675 from Eq. (2), obtaining a hyperspectral relative imaginary RI for the coat as chloroplast. In Eq. (2), Chl_i applies to the whole cell and is therefore scaled using V_s to obtain n^′(675) for the coat alone. The nominal chloroplast's relative real RI is then picked from a distribution as shown in Table 1 and modified as a function of its imaginary RI according to the Kramers–Kronig relations (implemented as a Hilbert transform) (Bernard et al., 2009; Robertson Lain and Bernard, 2018).

The cell cytoplasm is represented by the particle core. Its relative real RI is picked from a distribution given in Table 1, and it is modified by the Kramers–Kronig relations using a constant hyperspectral detritus-like imaginary RI, i.e., having a colored dissolved organic matter (CDOM)-like exponential spectral shape, resulting in a spectrally varying hyperspectral relative real RI. The phytoplankton particle population relative RIs and their Monte Carlo variability are summarized in Supplement Fig. S1.

Zhang et al. (2009)Boss et al., 2001(Morel and Bricaud, 1986; Bernard et al., 2009)Organelli et al., 2018(Bernard et al., 2009; Robertson Lain and Bernard, 2018)(Bernard et al., 2009; Robertson Lain and Bernard, 2018)Bernard et al., 2009Robertson Lain and Bernard, 2018(Bernard et al., 2009; Robertson Lain and Bernard, 2018)Morel et al., 1993

Table 1 Inputs for the coated-sphere Aden–Kerker optical scattering computations for the phytoplankton particle population. Modeling inputs common to both phytoplankton and NAPs (see Table 2) are given in the first three table rows. 𝒩(μ,σ) stands for a normal distribution with mean μ and standard deviation σ. For the indices of refraction, apart from Bernard et al. (2009) and Robertson Lain and Bernard (2018), see Morel and Bricaud (1986), Babin et al. (2003), Woźniak and Stramski (2004), and Duforêt-Gaurier et al. (2018).

(Bernard et al., 2009; Robertson Lain and Bernard, 2018)(Bernard et al., 2009; Robertson Lain and Bernard, 2018)Stramski and Kiefer, 1991Kostadinov et al., 2009Duforêt-Gaurier et al., 2018

Table 2 Inputs for the homogeneous-sphere Mie scattering code for the NAP population. Modeling inputs common to both phytoplankton and NAPs are given in the first three table rows of Table 1. 𝒩(μ,σ) stands for a normal distribution with mean μ and standard deviation σ.

The NAP population is represented by a homogeneous sphere, the relative RIs of which are picked so that its absorption spectrum is detritus-like (same as the core of phytoplankton), and its real RI is allowed to vary over a wider range of values, meant to represent mostly organic detritus, but with some minerogenic contributions, resulting in a mean nominal relative real RI of ≈ 1.06. The input RIs and other input parameters for NAPs are summarized in Table 2.

Specification of the input PSD parameters and the relationship of NAPs to phytoplankton PSDs is key to the construction of the forward and inverse models. Necessarily, some key simplifying assumptions are made here in order to construct an algorithm with operational application to modern multispectral ocean color sensors. The two key assumptions are that (1) phytoplankton and NAPs have a power-law PSD (Eq. 1) with the same slope ξ, and (2) the scaling parameter N₀ for NAPs is twice that of N₀ for phytoplankton (the forward model uses default values as in Tables 1 and 2). The latter assumption is chosen so that it results in a phyto C : POC ratio of 1:3 (see Kostadinov et al., 2016a; Behrenfeld et al., 2005; and Sect. 3.4 here) (as long as they are both estimated using the same size ranges). Together, these assumptions allow for the retrieval of one common PSD parameter set pertaining to the total particle population PSD (one ξ value and one total N₀ equal to the linear sum of the NAPs and phytoplankton N₀ values).

2.2 Backscattering calculations

The backscattering efficiencies, Q_bb(λ), for a single phytoplankton cell and a single non-algal particle were computed using the inputs described above in Sect. 2.1 and Tables 1 and 2. The coated-sphere code of Zhang et al. (2002) was used for both coated and homogeneous spheres. This code is included with the algorithm development scientific code of the PSD algorithm (see “Data availability”). Calculations were run for N=3000 instances of Monte Carlo simulations, each with a unique randomly picked combination of inputs for phytoplankton and NAPs. This resulted in 3000 sets of hyperspectral Q_bb values. High sampling resolution in diameter space was picked for the coated spheres (10 000 samples between minimum and maximum diameter) in order to minimize the influence of resonance spikes in Q_bb. For NAPs, 1000 samples of D were used.

Indices of refraction for both phytoplankton and NAPs are specified hyperspectrally (Supplement Fig. S1), and the computations are performed from 400 to 700 nm wavelength in vacuo with a step of 1 nm, allowing the resulting hyperspectral Q_bb(λ) values to be adapted for use with any combination of visible optical wavebands pertaining to recent and currently operating ocean color multispectral sensors or for planned (e.g., PACE Werdell et al., 2019) or existing hyperspectral sensors.

Before b_bp calculation, hyperspectral backscattering efficiencies, Q_bb, for each Monte Carlo run were first pre-processed by applying quality control and band-averaging using a moving-average 11 nm wide top-hat filter (using as central wavelengths the nominal bands of the following ocean color sensors: Sea-viewing Wide Field-of-view Sensor, SeaWiFS; Moderate Resolution Imaging Spectroradiometer, MODIS, Aqua; Medium Resolution Imaging Spectrometer, MERIS, and Ocean and Land Colour Instrument, OLCI; and the Visible and Infrared Imager/Radiometer Suite, VIIRS, on the Suomi National Polar-orbiting Partnership, S-NPP, plus 440 and 550 nm), resulting in 19 unique bands for band-averaged backscattering efficiencies, denoted here as $Q_{\mathrm{bb}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$. The band-averaged spectral particulate backscattering coefficient, $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$, was then calculated from the $Q_{\mathrm{bb}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ values and the input PSD as follows (e.g., van de Hulst, 1981; Kostadinov et al., 2009):

$$\begin{array}{}\text{(3)}& b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)=\underset{D_{\min }}{\overset{D_{\max }}{\int }}{\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}{D}^{\mathrm{2}}{Q}_{\mathrm{bb}}(D,\stackrel{\mathrm{‾}}{\mathit{\lambda }},m)N_{\mathrm{0}}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\xi }}\mathrm{d}D,\end{array}$$

where m is the complex index of refraction (specified separately for coat and core in the case of phytoplankton). Equation (3) is applied separately to the modeled phytoplankton and NAP Q_bb values and for each of the 3000 Monte Carlo runs. Band-averaged total $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ spectra are then calculated as the linear sum of phytoplankton and NAP backscattering.

2.3 PSD retrieval via spectral angle mapping

2.3.1 End-member construction

Band-averaged total $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ spectra were used to construct the backscattering end-members, $E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$, corresponding to specific input values of the PSD slope ξ. First, individual total b_bp spectra from each Monte Carlo run (N=3000) were normalized by the value at 555 nm. The median of all normalized spectra at each waveband was used as the end-member for each PSD slope, from ξ=2.5 to ξ=6 in steps of 0.05 (see Table 1). This approach allows the isolation of the b_bp spectral shape (dependent on ξ) and spectral magnitude (dependent on N₀) (Eq. 3). Using the hyperspectral underlying Q_bb values, end-members can be constructed for any desired set of wavelengths.

2.3.2 PSD parameter retrieval and operational application to OC-CCI ocean color data

The PSD parameters ξ and N₀ are retrieved using the backscattering end-members, $E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$, via the spectral angle mapping (SAM) technique (e.g., Dennison et al., 2004). Briefly, the end-members and satellite-observed b_bp spectra are treated as n-dimensional vectors, where n is the number of bands. The spectral angle between a given end-member and the observed spectrum is then calculated using the vector dot product as

$$\begin{array}{}\text{(4)}& \mathrm{\Theta }={\cos }^{-\mathrm{1}}\left({\displaystyle \frac{b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)\cdot E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)}{\parallel b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)\parallel \parallel E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)\parallel }}\right).\end{array}$$

Thus, spectral angle is an index of spectral shape similarity between two spectra, with more similar spectral shapes resulting in lower spectral angles. Equation (4) was used to calculate the spectral angle Θ between each of the 71 end-members, $E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$, and the input observed $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ spectrum. The value of ξ corresponding to the smallest spectral angle is then assigned as the retrieved PSD slope. Three wavebands were used, namely 490, 510, and 550 nm. For operational application to OC-CCI v5.0 (Sathyendranath et al., 2021) remote-sensing reflectance (R_rs(λ)) data (which do not have the 550 nm band), band-shifting was applied to the input R_rs(560) to estimate the corresponding R_rs(550), which is used in the Loisel and Stramski (2000) IOP inversion. The band-shifting was constructed using the band ratios between the respective original and target bands from a hyperspectral run of the Morel and Maritorena (2001) (MM01) model. No other bands were shifted.

The N₀ parameter is subsequently retrieved as the ratio of (1) the satellite-observed value of b_bp(443) and (2) the median value of the quantity $b_{\mathrm{bp}}\left(\mathrm{443}\right)/N_{\mathrm{0}}$ corresponding to the end-member class of the retrieved ξ and all statistically similar classes (see Supplement Sect. S1) across all Monte Carlo simulations.

2.4 Derived products: size-partitioned phytoplankton carbon, PSCs, POC, and chlorophyll

Once the PSD parameters are known, they can be used to compute derived products (Kostadinov et al., 2010, 2016a; Roy et al., 2017). Phytoplankton carbon in any size class spanning from cell diameter D_min to cell diameter D_max (m) can be estimated as

$$\begin{array}{}\text{(5)}& \mathrm{phyto}\mathrm{C}=\underset{D_{{min}_{\mathit{\varphi }}}}{\overset{D_{{max}_{\mathit{\varphi }}}}{\int }}{\mathrm{10}}^{-\mathrm{9}}a{\left({\mathrm{10}}^{\mathrm{18}}{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{D}^{\mathrm{3}}\right)}^b{N}_{\mathrm{0}\mathit{\varphi }}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\xi }}\mathrm{d}D,\end{array}$$

where $N_{\mathrm{0}\mathit{\varphi }}=\frac{\mathrm{1}}{\mathrm{3}}{N}_{\mathrm{0}}$, and N₀ (m⁻⁴) for the total PSD is the satellite-retrieved parameter from total particulate backscattering; the other PSD parameters are as in Eq. (1). Equation (5) was used to compute size-partitioned phyto C in three size classes – picophytoplankton (0.2 to 2 µm in diameter), nanophytoplankton (2 to 20 µm in diameter), and microphytoplankton (20 to 50 µm in diameter) – as well as total phyto C as the sum of the three classes. Carbon-based PSCs are defined as the fractional contribution of each of the three size classes to total phyto C (Kostadinov et al., 2016a). Given the first-order correspondence between PSCs and PFTs (e.g., Le Quéré et al., 2005), these PSCs can also be interpreted as PFTs. The allometric coefficients of Roy et al. (2017) are used here, namely a=0.54 and b=0.85; when cell volume V is expressed in cubic micrometers, cellular carbon is computed in picograms of carbon per cell using these coefficients (Eq. 5; see also Menden-Deuer and Lessard, 2000). Phyto C in Eq. 5 is given in milligrams per cubic meter; the conversion factors in Eq. (5) are used to convert from cubic meters to cubic micrometers and from picograms to milligrams of carbon (Kostadinov et al., 2016a; Roy et al., 2017). The factor of $\mathrm{1}/\mathrm{3}$ is an assumption of the model (Tables 1 and 2). Thus, an estimate of POC (computed using the same size limits as total phyto C) was calculated as 3× phyto C.

Chlorophyll concentration was estimated from the PSD retrievals and the input intracellular chlorophyll concentration, Chl_i (Table 1; Roy et al., 2017), as follows:

$$\begin{array}{}\text{(6)}& \mathrm{Chl}=\underset{D_{{min}_{\mathit{\varphi }}}}{\overset{D_{{max}_{\mathit{\varphi }}}}{\int }}{\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{D}^{\mathrm{3}}{\mathrm{Chl}}_{\mathrm{i}}{N}_{\mathrm{0}\mathit{\varphi }}{\left({\displaystyle \frac{D}{D_{\mathrm{0}}}}\right)}^{-\mathit{\xi }}\mathrm{d}D.\end{array}$$

In Eq. (6), Chl_i, D, D₀, and N_0ϕ all have to be expressed in consistent units so that Chl is obtained in milligrams per cubic meter. Here we use the median Chl_i across all Monte Carlo simulations to produce a single Chl estimate.

2.5 Validation and comparison

A data set of near-surface in situ PSD measurements was compiled for validation of the PSD parameter products, ξ and N₀ (Eq. 1). The data set consists of Coulter counter and Laser In-Situ Scattering and Transmissometry (LISST) measurements and a small set of PSDs derived from multiple instruments and modeling. Specifically, the compilation consists of the following data sets: (1) a compilation of several data sets of Coulter counter measurements, as used in the KSM09 algorithm validation in Kostadinov et al. (2009); (2) LISST-100X (Sequoia Scientific©) measurements from the Plumes and Blooms project (e.g., Toole and Siegel, 2001; Kostadinov et al., 2007) in the Santa Barbara Channel, as used in Kostadinov et al. (2012); (3) Coulter counter measurements from the Atlantic Meridional Transect cruise no. 26 (AMT26) (Organelli and Dall'Olmo, 2018), as compiled and used in Organelli et al. (2018, 2020); (4) in-line LISST 100-X measurements from the cruises North Atlantic Aerosols and Marine Ecosystems Study 3 and 4 (NAAMES; https://science.larc.nasa.gov/NAAMES/, last access: 25 January 2023; Boss and Haëntjens, 2017) and EXport Processes in the Ocean from Remote Sensing (EXPORTS; https://oceanexports.org/, last access: 7 October 2022) North Pacific (NP) and North Atlantic (NA) (Boss and Haëntjens, 2018) (data were downloaded from the NASA SeaBASS database; Werdell et al., 2003); and (5) PSDs obtained using a volume scattering function (VSF) inversion technique (Zhang et al., 2011, 2012) from VSFs measured during the NASA EXPORTS campaign (Siegel et al., 2016) in the North Pacific in 2018 (Siegel et al., 2021).

The compiled PSD data set was used to fit for the PSD parameters of Eq. (1) using the 2 to 20 µm diameter range. One data point was removed from the 2018 EXPORTS PSD data due to a poor fit to a power-law PSD. These in situ estimates were matched to satellite OC-CCI v5.0 (Sathyendranath et al., 2019, 2021) R_rs using the same matching methods described below for POC and picophytoplankton carbon data. Matched reflectances were used as input to the novel PSD algorithm presented here. The in situ and satellite PSD parameters were then compared using a type II linear regression and several additional algorithm performance metrics (e.g., Seegers et al., 2018), details of which are given in the caption of Fig. 8.

A large compilation of in situ POC data was collected from various public databases and private contributors and was used here to perform match-ups with satellite OC-CCI v5.0 data. In addition to the POC data (1997–2012) used in Evers-King et al. (2017) for algorithm validation (N=3891), this study also incorporated recent in situ POC data (2013–2020) from the SeaWiFS Bio-optical Archive and Storage System archive (Werdell et al., 2003) (https://seabass.gsfc.nasa.gov/, last access: 1 May 2021). The daily, 4 km, sinusoidal projection OC-CCI v5.0 data (1997–2020) (Sathyendranath et al., 2019, 2021) were used to extract the closest central satellite pixels to the in situ data points. If the central satellite match-up pixels were valid, the surrounding eight pixels (a 3×3 pixel box) were also extracted to estimate the mean, median, and standard deviation of all OC-CCI variables. The match-up data points were then averaged with respect to depth (0 to 10 m), location, and date. Moreover, a number of uncertain match-up data points were removed, as described below. A total number of 6041 match-up data points were obtained and used for analysis. Here, the median satellite R_rs(λ) matched-up spectra were used to compute the satellite-retrieved POC data using the PSD-based algorithm.

The in situ picophytoplankton carbon data set compiled and used for algorithm inter-comparison as part of the ESA POCO project (Martínez-Vicente et al., 2017) was used here to generate match-ups with satellite OC-CCI v5.0 R_rs data for further validation. Match-ups were generated in the same way as described above for POC.

All in situ data described above were excluded from the validation if any of the following conditions were met: (1) average bathymetric depth from a ≈ 9 km buffer around the in situ sample location was less than or equal to 200 m, or any grid cell elevations in that buffer were 0 m or higher, using a downsampled, 4 km version of the NOAA ETOPO1 data set (https://www.ngdc.noaa.gov/mgg/global/, last access: 23 October 2015); (2) the in situ sample depth was 15 m or greater; or (3) there were three or fewer satellite pixels available to use in the match-up, as detailed above.

All duplicate in situ match-ups (in the sense of multiple in situ data points that are close in space and time and receiving the same satellite match-up) were combined into a single match-up point as follows: for the PSD, the medians of the PSD measurements from the NAAMES and EXPORTS cruises in each LISST size bin for such duplicates were used for the calculation of in situ PSD parameters (a large number of duplicates since the data are in-line); for the rest of the PSD data, the fit PSD parameters themselves were averaged (a small number of duplicates). For the POC (large number of duplicates) and picophytoplankton carbon data (small number of duplicates), the averages of the duplicate in situ data were used.

In addition to validation against in situ measurements of the PSD, POC, and picophytoplankton carbon, satellite chlorophyll a (Chl) retrievals (using the standard algorithm of OC-CCI v5.0 at the match-up points) were compared with Chl estimated using the PSD retrieval (Eq. 6). Finally, global algorithm retrievals of total phyto C for May 2015 (using OC-CCI v5.0 data as input) were also compared with two alternative methods of retrieving phyto C: (1) the Roy et al. (2017) algorithm and (2) the Graff et al. (2012, 2015) algorithm, as implemented by NASA's Ocean Biology Processing Group (OBPG). Modeling and processing of results presented here are done using the sinusoidal projection images; maps presented here are given in equidistant cylindrical projection (i.e., unprojected latitude/longitude).

3 Results and discussion

3.1 Forward and inverse modeling

The first step in the algorithm development is the generation of 3000 Monte Carlo realizations of backscattering efficiencies as a function of particle diameter and wavelength, Q_bb(D,λ). The important differences between backscattering efficiencies of homogeneous and coated particles is discussed in Supplement Sect. S2. Here, we continue the discussion with the resulting integrated backscattering spectra (Eq. 3). Hyperspectral b_bp(λ) spectra modeled using a single forward optical model run are shown in Fig. 1. The computations use the median values of inputs varied in the Monte Carlo simulations (Tables 1 and 2). These normalized spectra illustrate the strong spectral shape dependence on the PSD slope ξ. Phytoplankton b_bp spectral shapes are complex, with various peaks and troughs near the absorption peaks of chlorophyll, but are more linear in the 490 to 550 nm range, which is the one used for the multispectral operational PSD algorithm. Regardless, the SAM methodology of retrieval used here allows for any spectral shape and does not impose a power-law fit to the shape of b_bp, as is done in KSM09 (Kostadinov et al., 2009) (see also Loisel et al., 2006). NAP backscattering exhibits smooth shapes due to the smooth shape of their absorption (Fig. 1b). Fundamentally, it is evident from Fig. 1a and b that the higher the PSD slope ξ, the steeper the b_bp spectral shape becomes, with higher values in the blue, since smaller particles dominate the signal. This dependence is at the root of the principle of operation of the PSD algorithm. For completeness, corresponding absorption spectra are illustrated in Supplement Fig. S3.

Figure 1 Modeled hyperspectral backscattering coefficient of (a) phytoplankton, using EAP-based coated-sphere scattering computations, and (b) NAPs, modeled as homogeneous spheres, as a function of the input power-law PSD slope (color-coded solid lines, as in legend). All spectra are shown normalized to the respective values at 555 nm. See Sect. 2 for more details.

The 71 end-members (EMs) created for operational application to existing major satellite ocean color missions and corresponding to PSD slope values between 2.5 and 6.0 with a step of 0.05 are displayed in Fig. 2a. They represent the modeled $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ spectra against which satellite-measured b_bp spectra are compared using the SAM method (Eq. 4). The spectral shape dependence on ξ demonstrates the theoretical ability to retrieve this parameter from space.

Figure 2 (a) Normalized spectral shapes of the b_bp end-members (EMs) developed for spectral angle mapping (SAM), shown at the 19 unique wavelengths used for band-averaging (see Sect. 2.2), and for PSD slopes as in the color legend. (b) The fraction of b_bp(λ) due to phytoplankton as a function of the PSD slope ξ. The wavelengths shown are indicated in the legend (nm). The means across all 3000 Monte Carlo simulations are shown. (c) Uncertainties in the PSD slope ξ retrieval using the SAM method, for each EM. Shown are the minimum and maximum value of the PSD slope for all end-member classes that are statistically similar to the given EM, according to the Kruskal–Wallis ANOVA (Supplement Sect. S1) (left y axis), and the resulting range of PSD slopes (right y axis) falling within these asymmetric uncertainty bounds. (d) Statistics of the parameter $\log \mathrm{10}\left(b_{\mathrm{bp}}\right(\mathrm{443})/N_{\mathrm{0}})$ for each EM, calculated for all 3000 Monte Carlo simulations and across all neighboring EM classes determined to be statistically similar to the given EM. μ in the legend stands for the mean, and σ is the standard deviation. The standard deviation of this parameter is used to estimate uncertainties in the N₀ retrieval.

An important question in bio-optical oceanography is determining the sources of backscattering in the ocean and their relative contributions. This is still an unresolved issue (Stramski et al., 2004), though progress has been made (e.g., Organelli et al., 2018; Koestner et al., 2020; Zhang et al., 2020). This issue is of central importance to the PSD model, as it assigns varying fractions of the b_bp signal to phytoplankton vs. NAPs under certain assumptions (Tables 1 and 2). Since in the two-component PSD model presented here phytoplankton and NAP b_bp are modeled separately, the fraction of b_bp due to phytoplankton vs. NAPs can be calculated. For a given PSD slope ξ and wavelength, the assumptions of the model dictate fixed fractional contributions by NAPs and phytoplankton to total b_bp, which are given in Fig. 2b. There is variability by wavelength, but the first-order variability is driven by the PSD slope, namely at low ξ values (ξ<4.0), phytoplankton contribution to b_bp is on the order of 30 % to 50 %, and it drops to near 0 % for higher slopes as ξ approaches 6.0. The curves are not monotonic, and peak phytoplankton contribution to b_bp occurs at ξ≈3.25.

The fractional contributions of Fig. 2b are derived from the forward theoretical modeling, and they are influenced by all model assumptions and are not validated independently. In particular, the decisions on integration diameter limits for NAPs and phytoplankton, as well as on the distributions of the indices of refraction for phytoplankton and NAPs, will have a strong influence on these values. Since NAPs are permitted here to have higher RIs than RIs typical of organic detritus only, if NAPs were strongly dominated by or composed only of organic particles, then NAP contribution to b_bp would be overestimated here. Of course these RIs are likely to be spatially and temporally variable, and the algorithm can be further improved by investigating and implementing such variability. Bellacicco et al. (2018) estimated global absolute b_bp due to NAPs and its fractional contribution to total b_bp using analysis of correlations with Chl. Qualitatively and to first order, their global pattern of percent b_bp due to NAPs agrees with the model results reported here, i.e., low relative NAP contributions at high latitudes and in eutrophic areas and higher relative contributions in more oligotrophic areas such as the fringes of the subtropical gyres (they exclude the gyres from their analysis) (cf. their Fig. 2c and 2b here). Note that the Bellacicco et al. (2018) estimate pertains only to NAPs non-covarying with Chl, making comparison harder. Further investigation is warranted to more rigorously compare their product to the values implicit in the PSD algorithm described here. Apart from analyzing the relative contribution of phytoplankton vs. NAPs to total b_bp, it is of interest to investigate the relative contributions of various size ranges to the modeled backscattering coefficient. This is illustrated in Supplement Fig. S4 and further discussed in Supplement Sect. S3.

The uncertainty in PSD slope retrieval is illustrated in Fig. 2c. These estimates are not symmetric about the ξ value and are derived via Kruskal–Wallis analysis of variance to determine class similarity (Supplement Sect. S1). As in KSM09, the general tendency is for the range of uncertainty in ξ to increase for lower PSD slopes, but it is always less than 0.5. The uncertainty in the $b_{\mathrm{bp}}\left(\mathrm{443}\right)/N_{\mathrm{0}}$ ratio used to retrieve the N₀ parameter is shown in Fig. 2d in log10 space. Mean and median values are similar, and the uncertainty about them does not vary much with PSD slope, also similarly to KSM09. The uncertainty in Fig. 2d at each ξ value includes all statistically similar classes of EMs.

3.2 Operational application of the PSD–phyto C algorithm to OC-CCI v5.0 merged satellite data

3.2.1 PSD parameters

The operational PSD algorithm presented here was applied to the monthly 4 km OC-CCI v5.0 R_rs(λ) data set (Sathyendranath et al., 2019, 2021). Both PSD parameters (ξ and N₀; Eq. 1) and derived products were generated (Sect. 2.3 and 2.4). These data and their monthly and overall climatologies (and associated uncertainties) are made publicly available (see “Data availability”). Here, we use May 2015 data to illustrate and discuss the new algorithm.

The PSD slope map (Fig. 3a) reveals a global spatial pattern consistent with expectations and with KSM09; namely the subtropical oligotrophic gyres are characterized by high PSD slopes (relatively high numerical dominance of small particles), whereas more eutrophic areas such as coastal areas, equatorial upwelling zones, and high latitudes exhibit lower slopes (increasing relative abundance of larger particles). This is consistent with oligotrophic ocean ecosystems being dominated by picophytoplankton, whereas microphytoplankton contribute significantly to the phytoplankton assemblage in eutrophic areas and during blooms (e.g., Kostadinov et al., 2009; Kostadinov et al., 2010; and references therein). PSD slope values retrieved by the SAM-based algorithm span the full modeled range of $\mathrm{2.5}\le \mathit{\xi }\le \mathrm{6.0}$. This is in contrast to KSM09, where values below 3.0 were not retrieved. The N₀ PSD parameter (Eq. 1) is, as expected, higher in coastal, high-latitude, and eutrophic areas (indicating higher particle loads) and lower in the oligotrophic subtropical gyres (Fig. 3b). N₀ varies over a few orders of magnitude. Note N₀'s units of m⁻⁴ (Eq. 1) and that care should be taken when comparing Eq. (1) and N₀ to other formulations of the PSD, e.g., the k parameter in Roy et al. (2017), as these are related but not equivalent (see also Vidondo et al., 1997).

Figure 3 Example operational retrievals of the PSD parameters (Eq. 1) and their uncertainties, using monthly OC-CCI v5.0 R_rs data for May 2015: (a) PSD slope ξ, (b) N₀ parameter (m⁻⁴ in log10 space), (c) uncertainty range for ξ, and (d) standard deviation (SD) of log10 of N₀.

Algorithm uncertainties are provided on a per-pixel basis. The uncertainty range estimates for ξ (Fig. 3c) (not necessarily symmetric about the ξ value) indicate that the gyres are characterized by lower uncertainties than the more eutrophic areas, as can be expected from Fig. 2c. These are partial uncertainty estimates, including those quantifiable and internal to the modeling, i.e., due to Mie parameter choices. Additional uncertainties inherent in the input OC-CCI R_rs values and those due to the IOP inversion algorithm used are not included in Fig. 3c and d and in subsequent propagated errors. Uncertainties in the N₀ parameter are more uniform spatially but higher in the gyres (Fig. 3d). Note that those are given in log10 space as a standard deviation, and a relatively small absolute value of the uncertainty translates to relatively large uncertainties in absolute particle concentrations.

3.2.2 Phytoplankton carbon and carbon-based PSCs, POC, and chlorophyll from the PSD

Global patterns of total phytoplankton carbon (phyto C) retrieved via the PSD and allometric relationships (Fig. 4a) exhibit the expected lower values in the oligotrophic gyres and higher values elsewhere. Similarly to the results of Kostadinov et al. (2016a), values range over approximately 3 orders of magnitude, which is a higher range than retrievals based on other methods, namely direct empirical algorithm POC retrieval (Stramski et al., 2008) or the Behrenfeld et al. (2005) method of scaling backscattering, and it is also higher than the range in CMIP5 model ensembles (cf. Fig. 1 in Kostadinov et al., 2016a). This putative underestimation in the gyres and overestimation in eutrophic areas suggests the need for algorithm tuning, which is discussed in Sect. 3.3 along with implications of validation results. Global validation of phyto C retrievals with analytical phyto C measurements is planned but is currently challenging as phytoplankton-specific carbon data are relatively novel (Graff et al., 2012, 2015) and still scarce. Here, an initial validation effort is undertaken using several other variables; see Sect. 3.3.

Figure 4 Example operational PSD-based retrievals of total and size-partitioned phytoplankton carbon, using monthly OC-CCI v5.0 R_rs data for May 2015: total phytoplankton carbon (a), picophytoplankton carbon (b), nanophytoplankton carbon (c), and microphytoplankton carbon (d). Units are milligrams per cubic meter, mapped in log10 space. The diameter limits for the three size classes are picophytoplankton (0.2 to 2 µm), nanophytoplankton (2 to 20 µm), and microphytoplankton (20 to 50 µm).

A key feature of the PSD-based algorithm is that phyto C can be partitioned into any number of size classes by choosing appropriate integration limits of Eq. (5). Absolute concentrations of pico-, nano-, and microphytoplankton are illustrated for May 2015 in Fig. 4b, c, and d, respectively. Picophytoplankton C is mapped on the same color scale as total phyto C (Fig. 4a), but pico- and nanophytoplankton C maps have differing scales, illustrating that while picophytoplankton C varies over ≈ 3 orders of magnitude spatially and globally, nanophytoplankton C varies over ≈ 4–5 orders of magnitude, and microphytoplankton varies over ≈ 7 orders of magnitude (see also Kostadinov et al., 2010; Kostadinov et al., 2016a). Note that empirical tuning will affect these ranges of variability; see Sect. 3.3. Fractional contributions of each of the three PSCs used here to total phyto C are illustrated in Fig. 5. Picophytoplankton dominate much of the open-ocean, lower-latitude, oligotrophic areas, contributing nearly 100 % of the carbon biomass there (Fig. 5a); nanophytoplankton (Fig. 5b) contribute up to ≈ 50 % of the biomass in the higher-latitude and more eutrophic areas; and microphytoplankton (Fig. 5c) contribute significantly only in the most eutrophic areas, e.g., during the North Atlantic bloom at ≈ 45–50^∘ N (May 2015 is shown). As previously noted (Kostadinov et al., 2010, 2016a), this general pattern is consistent with current understanding of ocean ecosystems.

Figure 5 Example operational retrieval of the percent contribution of each phytoplankton size class (PSC) to total phytoplankton carbon, determined via the PSD (Sect. 2.4). Retrievals use monthly OC-CCI v5.0 R_rs data for May 2015. The PSCs are picophytoplankton (a), nanophytoplankton (b), and microphytoplankton (c).

The fractional carbon-based PSCs (Fig. 5) are ratios of two integrals of Eq. (5); thus they are analytical functions of the PSD slope ξ and the b allometric coefficient (as well as the limits of integration used for each class and total phyto C). These functions are plotted in Fig. 6, together with the satellite-observed ξ histogram for May 2015, illustrating the most common values for the PSCs found in the ocean. Area-wise, the ocean is dominated by oligotrophic regions with high picophytoplankton contributions to C biomass.

Figure 6 Percent contribution of each PSC to total phytoplankton carbon (blue curves as in legend, left y axis) as a function of the PSD slope ξ. A histogram of the PSD slope from the (sinusoidal projection) OC-CCI v5.0-based image for May 2015 is shown in the background as brown bars (right y axis). The three PSC curves are analytically derived from the model, and no satellite data are used in producing them.

As an illustration of uncertainty propagation to derived products, the propagated uncertainty in total phyto C (Fig. 7a) and fractional picophytoplankton C biomass (Fig. 7b) is shown. Comparison of Fig. 4a with Fig. 7a indicates that absolute total phyto C uncertainties are of the same order of magnitude as the values themselves. This is a partial uncertainty estimation due to the assumed distributions of the Mie inputs (Tables 1 and 2) and due to the allometric coefficients. The Mie inputs are varied over wide ranges to accommodate various environments in the global ocean, with the goal of having a single first-principles-based operational algorithm applicable to first order globally. This increases the uncertainty estimates. The uncertainty for the fractional PSC products depends only on the uncertainties in ξ and b; thus they exhibit much lower internal algorithm uncertainty compared with absolute values. For picophytoplankton, they are $<\approx \mathrm{2}$ % for the oligotrophic gyres and do not exceed ≈ 7 % globally. This suggests that the fractional PSCs are more reliable products than the absolute values, and they can also be used with other products to partition them, e.g., total phytoplankton carbon estimated using the alternative methods shown in Fig. 9 (namely Graff et al., 2015, and Roy et al., 2017) or the Behrenfeld et al. (2005) or Sathyendranath et al. (2020) methods; POC products (e.g., Stramski et al., 2008) can be partitioned this way as well. Per-pixel uncertainties are estimated for all products and composite imagery (climatologies) as well and are provided with the OC-CCI-based data set associated with this paper (“Data availability”).

Figure 7 Propagated uncertainty in (a) total phytoplankton carbon given as 1 standard deviation (in milligrams per cubic meter), mapped in log10 space, and (b) fractional contribution of picophytoplankton to total phytoplankton carbon, given as 1 standard deviation in percent. (c) POC derived using the PSD retrievals (in milligrams per cubic meter), mapped in log10 space, and (d) chlorophyll concentration (Chl) derived using the PSD retrievals (in milligrams per cubic meter), mapped in log10 space. Monthly OC-CCI v5.0-based data for May 2015 are shown in all panels.

The formulation of the PSD algorithm allows for both POC and Chl (Eq. 6) to be estimated from the retrieved PSD. Due to the assumptions used, POC is phyto C multiplied by 3 (Fig. 7c). This is strictly true only if the POC estimate uses the same limits of integration as phyto C, which is an approximation of the usual POC operational definition (see, e.g., discussion of POC–PSD closure analysis in Kostadinov et al., 2016a). POC is thus estimated to first order, treating the retrieved NAPs as being composed of POC only and applying the same allometric relationship to NAPs as to phyto C, in spite of the fact that the assumed RI distribution of the NAPs is broader (Table 2). These are simplifying assumptions of the two-component model; a more accurate POC representation can be achieved if organic and inorganic NAPs are modeled as separate particle populations (e.g., Duforêt-Gaurier et al., 2018). This is a planned development of the model in the future; the goal here is to build an operational PSD–phyto C algorithm (based on first principles, as mechanistic as feasible) for use with multispectral satellite data of limited degrees of freedom. Hyperspectral sensors such as PACE (Werdell et al., 2019) should allow for some more degrees of freedom and thus for more independent particle components and their PSDs to be modeled separately. However, note that even hyperspectral data have limits on their degrees of freedom, which are expected to be much fewer than the number of sensor bands (Lee et al., 2007; Cael et al., 2020). An important benefit of POC is that it is a widely observed variable, available for global validation efforts (Sect. 3.3), as is Chl. Similarly to POC, there are benefits of the PSD-derived estimate of Chl (Fig. 7d); it can be used as additional verification/validation of model retrievals, and/or PSD-retrieved Chl can be used as a parameter to optimize (in algorithm tuning), as discussed shortly (Sect. 3.3). Next, we discuss validation/verification and tuning efforts in which both PSD-derived POC and PSD-derived Chl are used.

Figure 8 (a) Comparison of PSD slope derived from in situ measurements with the matched-up satellite retrieval (Sect. 2.5). Points are color-coded according to the corresponding matched-up satellite OC-CCI v5.0 Chl (colormap in milligrams per cubic meter in log10 space). In this figure (as well as in Figs. 10 and 11 and Supplement Figs. S6, S9, and S10), type II (reduced major axis, RMA) regression is used, and regression and validation statistics are given in the figure panels; “y-int” stands for the y intercept, RMS is root mean square (square root of the mean of squared differences between the in situ and satellite values), bias is the mean of the satellite minus in situ values, and MAE is the mean absolute error (the mean of the absolute values of the differences between the in situ and satellite values) (e.g., Seegers et al., 2018). The values in parentheses are the standard deviations of the slope and y intercept, respectively. (b) Same as in panel (a) for the N₀ parameter (Eq. 1) (axes in log10 space).

3.3 Algorithm validation, comparisons, and empirical tuning for the N₀ PSD parameter and absolute concentrations

In an initial validation effort, the novel PSD–phyto C algorithm is validated/verified using several variables. It is challenging to globally and thoroughly directly validate the major products of the algorithm – the PSD and size-partitioned phyto C – due to a paucity of globally spanning in situ observations, which are further reduced when performing satellite match-ups. Here, we validate or verify algorithm performance against compilations of the following variables: (1) in situ PSD observations (Sect. 2.5), (2) in situ POC observations, (3) in situ picophytoplankton C observations, and (4) concurrent satellite observations of Chl. Maps of the locations of in situ observations are shown in Supplement Fig. S5. In addition, we compare phyto C retrievals to several existing methods using the example May 2015 OC-CCI v5.0 image. Further, based on these results, we suggest an empirical tuning of the algorithm.

Validation results for the PSD slope ξ (Fig. 8a) indicate a statistically significant but noisy relationship between retrieved and observed slopes, with a positive bias for satellite retrievals and a regression slope substantially greater than unity. Most validation points are scattered in a cloud of data between 3.0 and 4.5 and do not exhibit much correlation, and there is a somewhat separate cluster of data centered about a slope of 5.25 in the satellite retrieval that have smaller corresponding in situ values of about 4.0. There is generally a clear tendency for points from more oligotrophic areas (as indicated by Chl color coding) to exhibit higher satellite values and more eutrophic areas to exhibit much lower satellite values. This tendency is weaker for the in situ observations, which tend to have a narrower range, mostly between 3.0 and 4.5. To first order, the satellite data are in the same range as in situ data, and the retrievals capture the in situ data trend; however, there is a pattern of having a bigger range of PSD slopes in satellite data than in the in situ match-ups, with the algorithm underestimating low values and overestimating high values.

Validation for the N₀ parameter (Fig. 8b) is statistically significant (somewhat higher R² than the ξ regression) but also quite noisy. Strong clustering of the in situ observations around 10^15.5–10^16.0 m⁻⁴ is observed, and the majority of these observations are somewhat overestimated in the satellite retrievals, which cluster around 10^16.25. Notably, a separate cluster associated with low Chl (and low values of both satellite and in situ N₀) shows that the satellite retrievals substantially underestimate these field data. Since N₀ is the PSD scaling parameter, which generally controls absolute number, volume, and carbon concentration variability to first order, this has implications for the global pattern of phytoplankton carbon retrievals (Fig. 4); namely it is consistent with underestimation in the oligotrophic gyres and overestimation in the eutrophic areas. Overall, both satellite and in situ data exhibit increasing values of N₀ with increasing Chl concentrations, as expected; i.e., more oligotrophic waters are associated with smaller overall particle number concentrations. Further discussion of the PSD validation by location of in situ data (Fig. S5a and b) is provided in Supplement Sect. S4 and illustrated in Supplement Fig. S6.

This pattern of under- and overestimation in the N₀ validation drives the slope of the validation regression to be much greater than unity and suggests an empirical tuning to absolute phytoplankton carbon estimates via a linear (in log10 space) tuning of N₀, as done in TK16 (Kostadinov et al., 2016a), who based the tuning on the validation regression. A similar approach is proposed here, but it is derived differently. Details of the tuning derivation procedure are given in Supplement Sect. S5. The following global tuning equation was obtained:

$$\begin{array}{}\text{(7)}& N_{\mathrm{0}}\mathit{\_}\mathrm{tuned}={\mathrm{10}}^{\mathrm{0.3859}\log \mathrm{10}\left(N_{\mathrm{0}}\right)+\mathrm{9.5531}},\end{array}$$

where N₀ is the original (untuned) PSD parameter. This tuning changes N₀ retrievals in a similar fashion to the TK16 tuning and is consistent with a tuning suggested by the N₀ in situ validation presented here (Fig. 8b); namely, low satellite N₀ values are increased, and high N₀ values are decreased, narrowing the overall range of variability of retrieved N₀ and thus the range of the retrieved derived variables. This addresses the low bias in oligotrophic gyres and the high bias in eutrophic areas. The goal of the tuning is to get more realistic absolute retrievals of POC and Chl (hypothesizing that this should also lead to more realistic phyto C retrievals; however, see discussion about picophytoplankton validation in Sect. 3.3.2).

The tuned N₀ parameter for May 2015 is mapped in Supplement Fig. S7a. The overall spatial pattern of higher values in more eutrophic areas is preserved, but the global range of values is reduced compared to the original formulation, increasing N₀ in the gyres and decreasing it in more productive areas. The resulting multiplicative factor to be applied in linear space to phyto C, POC, and Chl values is mapped in log10 space in Supplement Fig. S7b. Values in oligotrophic areas are generally greater than unity in linear space (mostly between 1 and 10), indicating that the tuning increases phyto C, POC, and Chl in these areas, up to about an order of magnitude (in limited areas mostly in the South Pacific Gyre) and more moderately elsewhere in the tropical and subtropical oligotrophic oceans. The equatorial upwelling areas and other transitional zones are not changed, and high-latitude oceans exhibit correction factors mostly less than unity in linear space (mostly between 0.1 and 1), which decreases phyto C and Chl up to an order of magnitude (rare, mostly less). This tuning is not applied to figures previously discussed here.

3.3.1 Comparison of the PSD-based phytoplankton carbon retrieval with existing satellite algorithms

In this section, we compare PSD-based phyto C retrievals presented here with two existing methods for its retrieval. The May 2015 original total phyto C retrieval is compared with the tuned total phyto C and the retrievals of the absorption- and PSD-based algorithm of Roy et al. (2017) and with the Graff et al. (2015) algorithm in Supplement Fig. S8. The histograms of these four images are compared in Fig. 9. The tuned retrievals are similar to those of Graff et al. (2015), whereas the original retrievals are similar to those of Roy et al. (2017), and the latter two have wider ranges globally compared to the former two. Of these algorithms, the simplest is that of Graff et al. (2015), as it is a direct scaling of b_bp, and it is based on in situ chemical analytical measurements of phyto C (Graff et al., 2012, 2015). These dichotomous inter-comparison results suggest that further algorithm inter-comparison and validation with direct in situ measurements of phyto C are needed to guide future algorithm developments; however these data are relatively novel and scarce globally. Validation results using in situ POC and picophytoplankton carbon (discussed next) exhibit a similar dichotomy.

Figure 9 Histograms of the images of Supplement Fig. S8, including the original PSD-based phyto C retrieval (Fig. 4a). Histogram counts are given on a log10 scale on the y axis, and the variable (x axis) is log-transformed as well. All four histograms are derived from the sinusoidal projection images for May 2015, using monthly OC-CCI v5.0 data.

3.3.2 Validation using POC and picophytoplankton carbon in situ data

PSD-based estimates of POC are validated against in situ measurements for the original algorithm (Fig. 10a) and the tuned algorithm (Fig. 10c). Both regressions have satisfactory R² values and also illustrate that in general higher POC values are associated with higher Chl (colormap). Notably, the original algorithm validation has a slope of ≈ 2 and exhibits substantial underestimates at low POC and overestimates at high POC. As intended, the tuning corrects this range exaggeration and significantly improves the slope, intercept, bias, RMS, and MAE. The regression with the N₀ tuning applied should not be considered a truly independent validation because the algorithm has been empirically tuned to retrieve POC well; however, the tuning was done with global POC imagery (using monthly images for 2004 and 2015) that uses the Stramski et al. (2008) empirical POC algorithm, not with these in situ POC data directly.

In addition to the validation with in situ POC, we performed a comparison of the matched satellite Chl and the corresponding PSD-based Chl estimate (Eq. 6) for the original (Fig. 10b) and the tuned algorithm (Fig. 10d). Both comparisons exhibit very high R² values, and similarly to POC, the original algorithm underestimated Chl at low values and overestimated Chl at high values. The tuning successfully addresses this, leading to excellent overall comparison of the tuned algorithm, with a slope near 1.0 and a low intercept. However, for the lowest Chl values (Chl <0.1 mg m⁻³), performance deteriorates. The tuned comparison is not a fully independent validation, as the algorithm was tuned to compare well with OC-CCIv5.0 satellite retrievals (using global monthly images for 2004 and 2015). Overall, the comparison with Chl is encouraging, indicating that the model is able to reasonably reproduce (with tuning) OC-CCI v5.0 standard satellite Chl values at the match-up points.

Figure 10 Validation of PSD-based POC retrievals vs. in situ POC measurements (a, c) and comparison of satellite retrievals of Chl using the standard OC-CCI v5.0 algorithm vs. Chl estimated from the PSD (b, d). Panels (a) and (b) have no tuning applied, whereas empirical tuning is applied to the N₀ parameter (Eq. 7) for panels (c) and (d). The tuning procedure applies a linear correction to N₀ in log10 space to ensure reasonable retrievals of POC and Chl, based on monthly satellite images from 2004 and 2015; for details, see Supplement Sect. S5. The data points in panels (a) and (c) are colored by their matched standard OC-CCI v5.0 satellite Chl values (colormap, in milligrams per cubic meter in log10 space).

Validation against in situ picophytoplankton carbon is presented in Fig. 11a (with no N₀ tuning applied) and in Fig. 11c with the tuning applied. The corresponding Chl comparisons between matched standard OC-CCIv5.0 Chl and Chl derived via the PSD model are shown in Fig. 11b and d. As with the POC match-ups (Fig. 10b and d), comparisons with Chl are better for the tuned version of the algorithm, indicating that the tuning is needed to reproduce more realistic Chl values globally. However, the tuning does not lead to any improvement in the validation results of picophytoplankton C (cf. Fig. 11a and c). The validation regression without tuning is statistically significant (p<0.05), albeit noisy (low R²=0.18); satellite retrievals and in situ data cover approximately the same ranges, and increasing Chl and in situ picophytoplankton C generally correspond to increasing satellite values as well, with some tendency for under- and overestimation as with the other variables. However, the tuned satellite retrievals have a very narrow range that does not cover the range of the in situ data, and validation statistics are generally worse than those of the original validation (the regression is not significant at the p=0.05 level). These validation results are generally consistent with the results of Martínez-Vicente et al. (2017), where the tuned version of the TK16 (Kostadinov et al., 2016a) algorithm was used.

Figure 11 Validation of picophytoplankton carbon derived from the PSD model using daily OC-CCI v5.0 satellite data vs. in situ measurements, as used in the POCO project (Martínez-Vicente et al., 2017) (a, c). Comparison of PSD-derived satellite Chl (y axes) with the matched satellite retrieval of Chl using the standard OC-CCI v5.0 algorithm at the locations of the in situ picophytoplankton carbon match-up points (x axes) (b, d). Panels (a) and (b) have no tuning applied, whereas empirical tuning is applied to the N₀ parameter (Eq. 7) for panels (c) and (d). The tuning procedure applies a linear correction to N₀ in log10 space to ensure reasonable retrievals of POC and Chl, based on monthly satellite images from 2004 and 2015; for details, see Supplement Sect. S5. The data points in panels (a) and (c) are colored by their matched standard OC-CCI v5.0 satellite Chl values (colormap, in milligrams per cubic meter in log10 space).

3.4 Further discussion, summary, and conclusions

The novel PSD–phyto C algorithm described here represents a major overhaul of the KSM09 algorithm (Kostadinov et al., 2009) (a comparison between KSM09 and the present algorithm is briefly discussed in Supplement Sect. S6). Unlike KSM09, two distinct particle populations are used: phytoplankton and NAPs. Phytoplankton backscattering is modeled using coated-sphere Mie calculations with inputs based on the Equivalent Algal Populations (EAP) approach (Bernard et al., 2009; Robertson Lain and Bernard, 2018). This model formulation allows assessment of the percent contribution of phytoplankton and NAPs to total b_bp as well as the estimation of Chl from the retrieved PSD. Underlying b_bp forward modeling is hyperspectral, facilitating adaptation of the algorithm to upcoming hyperspectral sensors like PACE (Werdell et al., 2019). PSD retrieval is achieved via spectral angle mapping (SAM), and no spectral shape is imposed on b_bp; operational end-members for current and past multispectral sensors and the OC-CCI v5.0 merged ocean color data set are created via band-averaging from the underlying hyperspectral modeled b_bp.

The algorithm has been used to create an accompanying data set based on the OC-CCI v5.0 data set (Kostadinov et al., 2022b; see “Data availability”). We emphasize that the PSD parameters and derived retrievals presented here and in the accompanying data set (Kostadinov et al., 2022b) are an experimental research satellite product with relatively large uncertainties. We do not claim that it is akin in validity and accuracy to the more established (and much more empirical!) algorithms for canonical products such as Chl and POC. As emphasized elsewhere in this text, the goal is to build an operational algorithm based on first principles as much as feasible, even at the expense of accuracy, in order to push the boundaries of what is retrievable from space and move the science of bio-optical algorithm development forward. Potential users of these PSD and derived data (Kostadinov et al., 2022b) need to be aware of their limitations, uncertainties, and validation status before using them, for example, in building or validating/constraining biogeochemical models. The choice of IOP algorithm to retrieve b_bp(λ) is key for the PSD–phyto C algorithm, as the spectral shape of b_bp is what the PSD slope retrieval is based upon (Eq. 4). The Loisel and Stramski (2000) IOP algorithm is chosen here, as in KSM09, because it allows spectral b_bp retrievals that are not constrained by a specific spectral function or parameterization of b_bp as is done, for example, in the QAA (quasi-analytical algorithm; Lee et al., 2002) and the GSM (Garver–Siegel–Maritorena) algorithm (Maritorena et al., 2002, 2010). For the wavelengths used in the PSD slope retrieval, modeled and satellite-derived b_bp spectral shapes compare well when the Loisel and Stramski (2000) algorithm is used, and global patterns of the retrieved PSD parameters appear reasonable. Preliminary tests with Loisel et al. (2018) indicate that this algorithm is not as suitable for PSD retrieval in this regard. Use of Loisel et al. (2018), Jorge et al. (2021), and other IOP algorithms will be further investigated in future development of the PSD algorithm.

An important assumption of the model is that N₀ for NAPs is twice that for phytoplankton so that the phyto-C-to-POC ratio is a constant 1:3. This ratio is expected to vary in the real ocean, and the value used here is a reasonable average choice (e.g., Behrenfeld et al., 2005; Jackson et al., 2017; Thomalla et al., 2017; and references therein). Graff et al. (2015) employed the cell sorting and chemical analysis methods of Graff et al. (2012) to measure phyto C in the equatorial Pacific and along the Atlantic Meridional Transect (AMT). Their results indicate that a phyto C : POC value of 1:3 is reasonable, falling within their observed ranges; however, they do observe many higher values, particularly in the oligotrophic gyres. The ratio of Roy et al. (2017) phyto C to Stramski et al. (2008) POC (applied to the May 2015 monthly image of the OC-CCI v5.0 data) indicates generally lower values of this ratio (with some high-latitude and coastal exceptions), and even lower values occur in the gyres, with values mostly below 0.1 in the low-latitude open ocean (data not shown, but consistent with work in progress by Shovonlal Roy). In light of this observation, note the difference between the Graff et al. (2015) and Roy et al. (2017) phyto C retrievals (Fig. 9 and Supplement Fig. S8). Further direct analytical observations of phyto C and the reconciliation and better understanding of the spatiotemporal variability in the phyto-C-to-POC ratio should be a high priority in order to improve understanding of carbon pools and their relationships in the ocean (Brewin et al., 2021) and to retrieve phyto C reliably from space.

The relatively poor PSD parameter validation results should be interpreted with caution, as there are multiple reasons for discrepancies between the in situ and satellite data and for the observed poor regression statistics, and the in situ data have their own limitations. Importantly, the in situ data PSD parameters are fit over a much narrower diameter range than the size range optically contributing the bulk of b_bp (see, e.g., Supplement Fig. S4), at least according to the modeled spectra. It is recognized that in the real world the particle assemblage is very complex, and its sources of backscattering are still not fully resolved (e.g., Stramski et al., 2004; Organelli et al., 2018). In particular, the composition and PSD of small sub-micron particles appear to be of importance and are not well known; here we assume the same PSD and NAP composition across all size classes and globally. There is also a mismatch in temporal and spatial scales of sampling between the satellite and in situ data. For example, the matched in situ PSD data do not exhibit the same negative correlation between ξ and N₀ that the satellite data do (Supplement Fig. S9). We note that this negative correlation in the satellite data has a theoretical underpinning because of what we know about global ocean ecosystems, namely that oligotrophic areas exhibit relative dominance of smaller phytoplankton (and smaller overall concentrations of particles/biomass), as opposed to increased importance of larger phytoplankton and increased biomass in more eutrophic areas. We thus expect backscattering in the ocean to become “bluer”, i.e., to have a steeper spectral slope, in oligotrophic areas. This is indeed observed in satellite data (e.g., Loisel et al., 2006) and is the basis for our algorithm. Therefore, in the ocean, we expect N₀ to decrease with increasing ξ globally and on average. This is not necessarily going to be captured by in situ data of limited spatiotemporal coverage and fit over a narrower size range.

We further note that the number of match-up points in the validation regression is different among PSD, POC, and picophytoplankton C, and their geographic distribution is different as well (Supplement Fig. S5). Namely, there are about an order of magnitude more POC match-ups than picophytoplankton carbon ones. Thus the different validation results presented here do not necessarily represent the same oceanographic conditions; e.g., the picophytoplankton C in situ data have less representation of eutrophic areas and span a smaller range of Chl than the POC validation, with very few points exceeding Chl = 1.0 mg m⁻³ (cf. Figs. 10b and 11b).

The picophytoplankton C data in Martínez-Vicente et al. (2017) are derived from cell counts (abundance) converted to carbon using specific conversion factors for different species/groups. Namely, 60 fg C per cell was used for Prochlorococcus, 154 fg per cell for Synechococcus, and 1319 fg per cell for pico-eukaryotes. This differs from the PSD-based phyto C retrieval algorithm in which the conversion is a function of cell volume and is continuous. For the allometric coefficients of Roy et al. (2017) used here, the equivalent conversion factor is ≈ 53 fg per cell for cells of the smallest diameter within the picophytoplankton range (0.5 µm) and is ≈ 1825 fg per cell for the largest-diameter cells within the picophytoplankton range (2.0 µm), indicating first-order consistency, but not full equivalency, with the methods of Martínez-Vicente et al. (2017).

The global relationships of the PSD parameters, derived PSD-based phyto C, and POC versus Chl for the May 2015 monthly image are illustrated in Supplement Fig. S10. Globally, as expected, there is a strong correlation of Chl with these variables, with increasing Chl associated with decreasing PSD slope and increasing N₀, phyto C, and POC. While the relationship is strong, there is significant spread of the PSD parameters and phyto C data for a given Chl value, suggesting that there is added value in retrieving them separately and that they should not all be treated as simply correlates of Chl. We note that there is a need for further investigation to avoid uniqueness of retrieval issues and degrees of freedom/independence issues, as well as more comprehensive and complete error propagation, since a lot of ecosystem properties are indeed correlated with Chl, and all these retrievals come from the same multispectral data.

The power law (Eq. 1) is a parameterization of real-world PSDs, and while there are theoretical underpinnings (e.g., West et al., 1997; Brown et al., 2004; Hatton et al., 2021) and observations (e.g., Quinones et al., 2003; Buonassissi and Dierssen, 2010; see also references in Boss et al., 2001) supporting its applicability, particularly over large size ranges, real-world PSDs may deviate from the power law, especially in coastal zones (e.g., Reynolds et al., 2010; Koestner et al., 2020; Runyan et al., 2020; Reynolds and Stramski, 2021). There is less information on living phytoplankton only and their specific PSDs because it has been historically difficult to separate living phytoplankton and measure, say, their PSD or carbon (e.g., Graff et al., 2012, 2015). A recent study by Haëntjens et al. (2022) investigates phytoplankton-specific PSDs; their observations support the conclusion that the phytoplankton-specific PSD shape is consistent with a power law to first order. We note that phytoplankton share their size domain with other organisms (bacteria on the low end and zooplankton at the high end), and we note that a drop-off in the phytoplankton PSD will be expected at the limits of the size range of autotrophs (see, e.g., Hatton et al., 2021); hence a phytoplankton-specific power law will have upper and lower range limits of applicability, and it is not expected to apply equally well over the same size range everywhere and always in the global ocean. Hatton et al. (2021) offer an assessment of the PSD of marine life over a huge range of sizes (body mass), demonstrating that a specific power law applies in the context of the Sheldon et al. (1972) hypothesis that equal biomass tends to occur in each logarithmically spaced size bin; their work offers support for use of the power law for modeling phytoplankton (over their size range) globally.

The power law is not a converging PSD model; i.e., it is sensitive to the chosen limits of integration (for a sensitivity analysis to the integration limits, see Kostadinov et al., 2016a). Gamma functions may be a better choice to represent marine PSDs (Risović, 1993, 2002). However, we choose to use the power law because of its theoretical underpinnings and because the goal is to build an operational algorithm (based on first principles as much as possible) for existing multispectral data with limited degrees of freedom. We additionally assume that the PSD slope for both phytoplankton and NAPs is the same, limiting the number of parameters to be retrieved. Hyperspectral data and observations of phytoplankton- and NAP-specific PSDs and IOPs will be needed to relax these assumptions in the future. Organelli et al. (2020) observed that the PSD slope steepened for small particles, deviating from a power law. This could partially explain the putative underestimates of the original algorithm in oligotrophic gyres. Moreover, the absolute number of particles retrieved is sensitive to uncertainties in the real index of refraction assumed. In this context, we note that the algorithm is able to pick up the concentration of particles, to first order, according to the N₀ validation (Fig. 8b). We find this to be impressive and consider it a success, given that the algorithm makes no a priori prescriptions about particle concentrations; they are solved for from the magnitude and shape of satellite-observed b_bp. While the goal here is to create a global algorithm which uses one set of end-members, we recognize that future implementations can be improved by assessing the impact of using regionally variable subsets of index-of-refraction distributions. The PSD parameterization and choices of Mie inputs, in particular complex indices of refraction, represent important sources of uncertainty and can also affect the need for tuning and the degree of suitability of estimating POC with our generic NAP population. Further algorithm analysis of performance and improvements need to focus on the index-of-refraction choices for the particle populations. For further discussion of algorithm uncertainties, see Kostadinov et al. (2009), Kostadinov et al. (2010), and Kostadinov et al. (2016a).

Graff et al. (2015) observe a relationship between phyto C and b_bp that is stronger than that for other proxies. This is encouraging for the use of backscattering as a proxy for phytoplankton carbon biomass. However, the link between the PSD and b_bp spectral shape is a second-order effect that is not easily observed in in situ observations (Kostadinov et al., 2009; Slade and Boss, 2015; Boss et al., 2018; Organelli et al., 2020), even though theoretical modeling demonstrates a clear link (Kostadinov et al., 2009; this study). Kostadinov et al. (2012) discuss some reasons why it may be difficult to observe this relationship in current in situ data, e.g., the fact that the PSD is fit over a narrow range of diameters compared to the size range theoretically affecting b_bp. Nevertheless, these considerations and the overall performance of the KSM09 homogeneous algorithm as compared to the algorithm presented here leads to the conclusion that there are four primary directions that should be priorities for moving forward. First, investigate the effect of choices of index-of-refraction distributions, as discussed above. Second, rather than relying only on b_bp for PSD and phyto C retrieval, a blended approach should be developed that also uses absorption; i.e., combine the approach here with that of Roy et al. (2017). Third, investigate the ability of hyperspectral data to provide more degrees of freedom for retrieval of more variables simultaneously, allowing relaxation of some key assumptions and perhaps a third particle population to represent POC and mineral particles separately; this is important in light of the upcoming PACE mission (Werdell et al., 2019). Hyperspectral absorption data in particular have the potential to increase information content and allow group-specific retrievals (e.g., Kramer et al., 2022; but see also Cael et al., 2020). Finally, collect more global, comprehensive in situ data sets of all relevant variables, including and especially phyto C (Graff et al., 2015), for further model development and validation. With regard to the latter, agencies and investigators should focus on building quality-controlled, one-stop-shop data sets.

Appendix A: Details on the OC-CCI v5.0 data set

Processing and analysis were done using the sinusoidal projection of OC-CCI v5.0. For user convenience, once the final products were generated, they were re-projected to equidistant cylindrical projection (unprojected latitude/longitude) before publication in the data repository linked above (“Data availability”). The empirical tuning (Sect. 3.3) is not applied to the variables in the published data set (“Data availability”). Instead, the spatially explicit linear-space multiplicative tuning factor (Supplement Fig. S7b) is given. The choice to provide an optional tuning to be applied at the user's discretion is dictated by the validation and comparison results discussed in the paper. Monthly and overall climatologies with propagated uncertainties are also provided, and for these climatologies, both tuned and original variables are included.