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

Abstract

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., ). Ocean ecosystems play a key role in Earth’s carbon cycle and climate by affecting atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ via the biological carbon pump, which sequesters some of the fixed carbon to the deeper ocean for longer timescales (e.g., ; ; ; ; ). The biological pump is influenced by the structure and function of oceanic ecosystems (e.g., ; ); therefore, mechanistic, predictive understanding of ocean ecosystems is of high priority to Earth systems and climate research (e.g., ). Satellite remote sensing of ocean color is a key tool for the global characterization of ocean ecology (e.g., ). This has led to large efforts to elucidate biological pump mechanisms using multiple platforms, including satellites, e.g., the EXPORTS program .

Phytoplankton cell size (diameters varying from $\approx$ 0.5 to $>\mathrm{50}$ $\mathrm{\mu }\mathrm{m}$ (e.g., ) is a key trait that affects multiple phytoplankton characteristics (), as well as sinking rates (e.g.,  ; ; ; ). Phytoplankton size classes (PSCs) thus tend to closely correspond to phytoplankton functional types (PFTs; e.g., ). 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., ; ; ; ). 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., ). 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.,  ; ; ; ; ; ; ):

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 }},$$ where $D$ is particle diameter; $N$ (${\mathrm{m}}^{-\mathrm{4}}$) is the differential number concentration of particles per unit volume seawater and per bin width of particle diameter; $N_{\mathrm{0}}=N\left(\mathrm{2}\mathrm{\mu }\mathrm{m}\right)$ is the particle number concentration at a reference diameter, here $D_{\mathrm{0}}=\mathrm{2}$ $\mathrm{\mu }\mathrm{m}$; and $\mathit{\xi }$ is the power-law slope of the PSD. Equation () has to be integrated over a given diameter range to get the total particle number concentration in that range, $N_{\mathrm{T}}$ (${\mathrm{m}}^{-\mathrm{3}}$).

Ocean color is quantified by the spectral shape and magnitude of the remote-sensing reflectance, $R_{\mathrm{rs}}\left(\mathit{\lambda }\right)$ (${\mathrm{sr}}^{-\mathrm{1}}$; also denoted simply $R_{rs}$ for brevity below), where $\mathit{\lambda }$ is the wavelength of light in vacuo. The Kostadinov–Siegel–Maritorena 2009 (KSM09, ) algorithm retrieves the parameters of an assumed power-law PSD ($\mathit{\xi }$ and $N_{\mathrm{0}}$ in Eq. ) from ocean color remote-sensing observations using the spectral shape and magnitude of the particulate backscattering coefficient, $b_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ (${\mathrm{m}}^{-\mathrm{1}}$). $b_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ can be retrieved using existing inherent optical property (IOP) inversion algorithms; KSM09 uses the IOP inversion. Subsequently, the retrieved PSD parameters allow the quantification of absolute and fractional PSCs: picophytoplankton, nanophytoplankton, and microphytoplankton based on bio-volume or phytoplankton carbon (henceforth TK16) via allometric relationships . Phytoplankton carbon (phyto C) is the key variable of interest for carbon cycle and climate studies and modeling, and TK16 (data set available: ) represents a relatively unique carbon-based approach among PSC/PFT algorithms as it is based on knowledge of the PSD and allometric relationships to get at size-partitioned phyto C. and 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 ; ). Recent results suggest that the structural complexity of oceanic particles enhances backscattering significantly and can explain the so-called “missing backscattering” in the ocean , i.e., the lack of optical closure between theoretically modeled and measured $b_{\mathrm{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., ), and they have significantly enhanced backscattering compared to their homogeneous equivalents .

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 . EAP explicitly models intracellular chlorophyll concentration, Chl${}_{\mathrm{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 , and in fact, sub-micron NAPs even smaller than the smallest autotroph ($\approx \mathrm{0.5}$ $\mathrm{\mu }\mathrm{m}$ in diameter) are critical for determining the spectral shape of $b_{\mathrm{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 in spite of evidence that phytoplankton/larger particles contribute more than Mie theory predicts, based on homogeneous spheres (e.g., ). 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_{\mathrm{bp}}$ of oceanic waters can be modeled (e.g.,  ; ; ). 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_{\mathrm{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 ; 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 for homogeneous spherical particles and the Aden–Kerker 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 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${}_{\mathrm{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  and .

As in the EAP model, the chloroplast is represented by the particle coat. Its relative volume, $V_{\mathrm{s}}$, is picked from a distribution as shown in Table . The chloroplast's imaginary refractive index (RI) (relative to seawater) at 675 $\mathrm{nm}$, $n^{\prime }\left(\mathrm{675}\right)$, is then computed as follows :

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)}},$$ where Chl${}^{*}=\mathrm{0.027}$ ${\mathrm{m}}^{\mathrm{2}}{\mathrm{mg}}^{-\mathrm{1}}$ is the theoretical maximum specific absorption coefficient of chlorophyll at 675 $\mathrm{nm}$ when dissolved in water , Chl${}_{\mathrm{i}}$ is the intracellular chlorophyll concentration in kilograms of chlorophyll per cubic meter of cellular material, and $n_{\mathrm{sw}}\left(\mathrm{675}\right)$ is seawater's absolute real RI at 675 $\mathrm{nm}$. A hyperspectral basis vector from the EAP model (based on measurements; for details see ; ) is then scaled using the value at 675 from Eq. (), obtaining a hyperspectral relative imaginary RI for the coat as chloroplast. In Eq. (), Chl${}_{\mathrm{i}}$ applies to the whole cell and is therefore scaled using $V_{\mathrm{s}}$ to obtain $n^{\prime }\left(\mathrm{675}\right)$ for the coat alone. The nominal chloroplast's relative real RI is then picked from a distribution as shown in Table  and modified as a function of its imaginary RI according to the Kramers–Kronig relations (implemented as a Hilbert transform) .

The cell cytoplasm is represented by the particle core. Its relative real RI is picked from a distribution given in Table , 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.

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 ) are given in the first three table rows. $\mathcal{N}(\mathit{\mu },\mathit{\sigma })$ stands for a normal distribution with mean $\mathit{\mu }$ and standard deviation $\mathit{\sigma }$. For the indices of refraction, apart from and , see , , , and .

Input parameter	Symbol, units, and notes
Pure seawater absolute real RI	$n_{\mathrm{sw}}=f\left(\mathit{\lambda }\right)$, after , usingtemperature $=$ 15 ${}^{\circ }$C and salinity = 33
PSD slope	$\mathit{\xi }\in$ [2.5, 6] in steps of 0.05 (the same value applies to both phytoplankton and NAPs; for the slope range,see, e.g., )
Wavelengths in vacuo	$\mathit{\lambda }\in$ [400, 700] $\mathrm{nm}$, hyperspectral – in steps of 1 $\mathrm{nm}$, band-averaging used for the nominal wavelengths of satellite sensors
Phytoplankton population inputs
Intracellular chlorophyll concentration	Chl${}_{\mathrm{i}}\in$ [0.5, 10], picked from, $\mathcal{N}$(2.5,2.5); $\mathit{\mu }\approx \mathrm{3.14}$ $\mathrm{kg}$ Chl ${\mathrm{m}}^{-\mathrm{3}}$ cellular material
Coat (chloroplast) relative volume	$V_{\mathrm{s}}\in$ [5, 35] %, picked from $\mathcal{N}$(20,5), resulting in mean coat relative thickness as fraction of cell radius $t_{\mathrm{coat}}=\mathrm{7.2}\mathit{\%}$ (cf. ); $t_{\mathrm{coat}}=\mathrm{1}-(\mathrm{1}-V_{\mathrm{s}}{)}^{\mathrm{1}/\mathrm{3}}$
Coat (chloroplast) relative real RI	$n_{\mathrm{coat}}\in$ [1.06, 1.22], picked from, $\mathcal{N}$(1.14,0.08); $\mathit{\mu }\approx \mathrm{1.14}$; wavelength-dependent via Kramers–Kronig relations
Core (cytoplasm) relative real RI	$n_{\mathrm{core}}\in$ [1.01, 1.03], picked from, $\mathcal{N}$(1.02,0.01); $\mathit{\mu }\approx \mathrm{1.02}$; wavelength-dependent via Kramers–Kronig relations
Coat (chloroplast) relative imaginary RI	$n_{\mathrm{coat}}^{\prime }\left(\mathit{\lambda }\right)$ computed from a hyperspectral basis vector (from ; ) that is scaled to the value at $n^{\prime }\left(\mathrm{675}\right)$ using Eq. ()
Core (cytoplasm) relative imaginary RI	$n_{\mathrm{core}}^{\prime }\left(\mathit{\lambda }\right)$ with a prescribed constant magnitude and detritus-like (exponential) spectral shape, with spectral slope $S=\mathrm{0.0123}$ ${\mathrm{nm}}^{-\mathrm{1}}$, resulting in $S_{\mathrm{a}}\approx \mathrm{0.014}$ ${\mathrm{nm}}^{-\mathrm{1}}$ for $a_{\mathrm{core}}\left(\mathit{\lambda }\right)$
Minimum outer particle diameter	$D_{{\min }_{\mathit{\varphi }}}=\mathrm{0.5}$ $\mathrm{\mu }\mathrm{m}$ (smallest autotroph; e.g.,  )
Maximum outer particle diameter	$D_{{\max }_{\mathit{\varphi }}}\in$ [20, 200] $\mathrm{\mu }\mathrm{m}$, picked from, $\mathcal{N}$(50,50); $\mathit{\mu }\approx \mathrm{72.3}$ $\mathrm{\mu }\mathrm{m}$
Differential number concentration at $D_{\mathrm{0}}=\mathrm{2}$ $\mathrm{\mu }\mathrm{m}$	$N_{\mathrm{0}\mathit{\varphi }}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{16}}$ ${\mathrm{m}}^{-\mathrm{4}}$, used in the forward modeling

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 . $\mathcal{N}(\mathit{\mu },\mathit{\sigma })$ stands for a normal distribution with mean $\mathit{\mu }$ and standard deviation $\mathit{\sigma }$.

NAP population inputs
Relative real RI	$n_{\mathrm{NAP}}\in$ [1.01, 1.2], picked from, $\mathcal{N}$(1.02,0.06); $\mathit{\mu }\approx \mathrm{1.06}$; wavelength-dependent via Kramers–Kronig relations
Relative imaginary RI	$n_{\mathrm{NAP}}^{\prime }\left(\mathit{\lambda }\right)$ with a prescribed constant magnitude and detritus-like (exponential) spectral shape , with spectral slope $S=\mathrm{0.0123}$ ${\mathrm{nm}}^{-\mathrm{1}}$, resulting in $S_{\mathrm{a}}\approx \mathrm{0.014}$ ${\mathrm{nm}}^{-\mathrm{1}}$ for $a_{\mathrm{NAP}}\left(\mathit{\lambda }\right)$
Minimum particle diameter	$D_{\min \mathrm{NAP}}=\mathrm{0.01}$ $\mathrm{\mu }\mathrm{m}$ (see Supplement Fig. S4 and ; )
Maximum particle diameter	$D_{\max \mathrm{NAP}}\in$ [200, 500] $\mathrm{\mu }\mathrm{m}$ (e.g.,  ), picked from, $\mathcal{N}$(400,10); $\mathit{\mu }\approx \mathrm{376.8}$ $\mathrm{\mu }\mathrm{m}$
Differential number concentration at $D_{\mathrm{0}}=\mathrm{2}$ $\mathrm{\mu }\mathrm{m}$	$N_{\mathrm{0}\mathrm{NAP}}=\mathrm{2}\times N_{\mathrm{0}\mathit{\varphi }}=\mathrm{1.0}\times {\mathrm{10}}^{\mathrm{17}}$ ${\mathrm{m}}^{-\mathrm{4}}$, resulting in total particle population $N_{\mathrm{0}}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{17}}$ ${\mathrm{m}}^{-\mathrm{4}}$, used in the forward modeling ($N_{\mathrm{0}\mathrm{NAP}}=\mathrm{2}\times N_{\mathrm{0}\mathit{\varphi }}$ constitutes an important model assumption, and it is discussed in the text)

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 $\approx$ 1.06. The input RIs and other input parameters for NAPs are summarized in Table .

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. ) with the same slope $\mathit{\xi }$, and (2) the scaling parameter $N_{\mathrm{0}}$ for NAPs is twice that of $N_{\mathrm{0}}$ for phytoplankton (the forward model uses default values as in Tables  and ). The latter assumption is chosen so that it results in a phyto C : POC ratio of $\mathrm{1}:\mathrm{3}$ (see ; ; and Sect.  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 $\mathit{\xi }$ value and one total $N_{\mathrm{0}}$ equal to the linear sum of the NAPs and phytoplankton $N_{\mathrm{0}}$ values).

2.2 Backscattering calculations

The backscattering efficiencies, $Q_{\mathrm{bb}}\left(\mathit{\lambda }\right)$, for a single phytoplankton cell and a single non-algal particle were computed using the inputs described above in Sect.  and Tables  and . The coated-sphere code of 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=\mathrm{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_{\mathrm{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_{\mathrm{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 $\mathrm{nm}$ wavelength in vacuo with a step of 1 $\mathrm{nm}$, allowing the resulting hyperspectral $Q_{\mathrm{bb}}\left(\mathit{\lambda }\right)$ 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 ) or existing hyperspectral sensors.

Before $b_{\mathrm{bp}}$ calculation, hyperspectral backscattering efficiencies, $Q_{\mathrm{bb}}$, for each Monte Carlo run were first pre-processed by applying quality control and band-averaging using a moving-average 11 $\mathrm{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 $\mathrm{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., ; ):

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,$$ where $m$ is the complex index of refraction (specified separately for coat and core in the case of phytoplankton). Equation () is applied separately to the modeled phytoplankton and NAP $Q_{\mathrm{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 $\mathit{\xi }$. First, individual total $b_{\mathrm{bp}}$ spectra from each Monte Carlo run ($N=\mathrm{3000}$) were normalized by the value at 555 $\mathrm{nm}$. The median of all normalized spectra at each waveband was used as the end-member for each PSD slope, from $\mathit{\xi }=\mathrm{2.5}$ to $\mathit{\xi }=\mathrm{6}$ in steps of 0.05 (see Table ). This approach allows the isolation of the $b_{\mathrm{bp}}$ spectral shape (dependent on $\mathit{\xi }$) and spectral magnitude (dependent on $N_{\mathrm{0}}$) (Eq. ). Using the hyperspectral underlying $Q_{\mathrm{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 $\mathit{\xi }$ and $N_{\mathrm{0}}$ are retrieved using the backscattering end-members, $E\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$, via the spectral angle mapping (SAM) technique (e.g., ). Briefly, the end-members and satellite-observed $b_{\mathrm{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

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).$$

Thus, spectral angle is an index of spectral shape similarity between two spectra, with more similar spectral shapes resulting in lower spectral angles. Equation () was used to calculate the spectral angle $\mathrm{\Theta }$ 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 $\mathit{\xi }$ corresponding to the smallest spectral angle is then assigned as the retrieved PSD slope. Three wavebands were used, namely 490, 510, and 550 $\mathrm{nm}$. For operational application to OC-CCI v5.0 remote-sensing reflectance ($R_{rs}\left(\mathit{\lambda }\right)$) data (which do not have the 550 $\mathrm{nm}$ band), band-shifting was applied to the input $R_{\mathrm{rs}}\left(\mathrm{560}\right)$ to estimate the corresponding $R_{\mathrm{rs}}\left(\mathrm{550}\right)$, which is used in the IOP inversion. The band-shifting was constructed using the band ratios between the respective original and target bands from a hyperspectral run of the (MM01) model. No other bands were shifted.

The $N_{\mathrm{0}}$ parameter is subsequently retrieved as the ratio of (1) the satellite-observed value of $b_{\mathrm{bp}}\left(\mathrm{443}\right)$ 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 $\mathit{\xi }$ 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 . Phytoplankton carbon in any size class spanning from cell diameter $D_{\min }$ to cell diameter $D_{\max }$ ($\mathrm{m}$) can be estimated as

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,$$ where $N_{\mathrm{0}\mathit{\varphi }}=\frac{\mathrm{1}}{\mathrm{3}}{N}_{\mathrm{0}}$, and $N_{\mathrm{0}}$ (${\mathrm{m}}^{-\mathrm{4}}$) for the total PSD is the satellite-retrieved parameter from total particulate backscattering; the other PSD parameters are as in Eq. (). Equation () was used to compute size-partitioned phyto C in three size classes – picophytoplankton (0.2 to 2 $\mathrm{\mu }\mathrm{m}$ in diameter), nanophytoplankton (2 to 20 $\mathrm{\mu }\mathrm{m}$ in diameter), and microphytoplankton (20 to 50 $\mathrm{\mu }\mathrm{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 . Given the first-order correspondence between PSCs and PFTs (e.g., ), these PSCs can also be interpreted as PFTs. The allometric coefficients of are used here, namely $a=\mathrm{0.54}$ and $b=\mathrm{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. ; see also ). Phyto C in Eq.  is given in milligrams per cubic meter; the conversion factors in Eq. () are used to convert from cubic meters to cubic micrometers and from picograms to milligrams of carbon . The factor of $\mathrm{1}/\mathrm{3}$ is an assumption of the model (Tables  and ). Thus, an estimate of POC (computed using the same size limits as total phyto C) was calculated as 3$\times$ phyto C.

Chlorophyll concentration was estimated from the PSD retrievals and the input intracellular chlorophyll concentration, Chl${}_{\mathrm{i}}$ (Table ; ), as follows: 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.$$

In Eq. (6), Chl${}_{\mathrm{i}}$, $D$, $D_{\mathrm{0}}$, and $N_{\mathrm{0}\mathit{\varphi }}$ all have to be expressed in consistent units so that Chl is obtained in milligrams per cubic meter. Here we use the median Chl${}_{\mathrm{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, $\mathit{\xi }$ and $N_{\mathrm{0}}$ (Eq. ). 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 ; (2) LISST-100X (Sequoia Scientific©) measurements from the Plumes and Blooms project (e.g., ; ) in the Santa Barbara Channel, as used in ; (3) Coulter counter measurements from the Atlantic Meridional Transect cruise no. 26 (AMT26) , as compiled and used in ; (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; ) 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) (data were downloaded from the NASA SeaBASS database; ); and (5) PSDs obtained using a volume scattering function (VSF) inversion technique from VSFs measured during the NASA EXPORTS campaign in the North Pacific in 2018 .

The compiled PSD data set was used to fit for the PSD parameters of Eq. () using the 2 to 20 $\mathrm{\mu }\mathrm{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 $R_{\mathrm{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., ), details of which are given in the caption of Fig. .

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 for algorithm validation ($N=\mathrm{3891}$), this study also incorporated recent in situ POC data (2013–2020) from the SeaWiFS Bio-optical Archive and Storage System archive (https://seabass.gsfc.nasa.gov/, last access: 1 May 2021). The daily, 4 $\mathrm{km}$, sinusoidal projection OC-CCI v5.0 data (1997–2020) 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 $\mathrm{3}\times \mathrm{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 $\mathrm{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_{\mathrm{rs}}\left(\mathit{\lambda }\right)$ 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 was used here to generate match-ups with satellite OC-CCI v5.0 $R_{\mathrm{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 $\approx$9 $\mathrm{km}$ buffer around the in situ sample location was less than or equal to 200 $\mathrm{m}$, or any grid cell elevations in that buffer were 0 $\mathrm{m}$ or higher, using a downsampled, 4 $\mathrm{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 $\mathrm{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. ). 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 algorithm and (2) the 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_{\mathrm{bb}}(D,\mathit{\lambda })$. 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. ). Hyperspectral $b_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ spectra modeled using a single forward optical model run are shown in Fig. . The computations use the median values of inputs varied in the Monte Carlo simulations (Tables  and ). These normalized spectra illustrate the strong spectral shape dependence on the PSD slope $\mathit{\xi }$. Phytoplankton $b_{\mathrm{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 $\mathrm{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_{\mathrm{bp}}$, as is done in KSM09 (see also ). NAP backscattering exhibits smooth shapes due to the smooth shape of their absorption (Fig. b). Fundamentally, it is evident from Fig. a and b that the higher the PSD slope $\mathit{\xi }$, the steeper the $b_{\mathrm{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 $\mathrm{nm}$. See Sect.  for more details.

[Figure omitted. See PDF]

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. a. They represent the modeled $b_{\mathrm{bp}}\left(\stackrel{\mathrm{‾}}{\mathit{\lambda }}\right)$ spectra against which satellite-measured $b_{\mathrm{bp}}$ spectra are compared using the SAM method (Eq. ). The spectral shape dependence on $\mathit{\xi }$ demonstrates the theoretical ability to retrieve this parameter from space.

Figure 2

(a) Normalized spectral shapes of the $b_{\mathrm{bp}}$ end-members (EMs) developed for spectral angle mapping (SAM), shown at the 19 unique wavelengths used for band-averaging (see Sect. ), and for PSD slopes as in the color legend. (b) The fraction of $b_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ due to phytoplankton as a function of the PSD slope $\mathit{\xi }$. 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 $\mathit{\xi }$ 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. $\mathit{\mu }$ in the legend stands for the mean, and $\mathit{\sigma }$ is the standard deviation. The standard deviation of this parameter is used to estimate uncertainties in the $N_{\mathrm{0}}$ retrieval.

[Figure omitted. See PDF]

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 , though progress has been made (e.g., ; ; ). This issue is of central importance to the PSD model, as it assigns varying fractions of the $b_{\mathrm{bp}}$ signal to phytoplankton vs. NAPs under certain assumptions (Tables  and ). Since in the two-component PSD model presented here phytoplankton and NAP $b_{\mathrm{bp}}$ are modeled separately, the fraction of $b_{\mathrm{bp}}$ due to phytoplankton vs. NAPs can be calculated. For a given PSD slope $\mathit{\xi }$ and wavelength, the assumptions of the model dictate fixed fractional contributions by NAPs and phytoplankton to total $b_{\mathrm{bp}}$, which are given in Fig. b. There is variability by wavelength, but the first-order variability is driven by the PSD slope, namely at low $\mathit{\xi }$ values ($\mathit{\xi }<\mathrm{4.0}$), phytoplankton contribution to $b_{\mathrm{bp}}$ is on the order of 30 % to 50 %, and it drops to near 0 % for higher slopes as $\mathit{\xi }$ approaches 6.0. The curves are not monotonic, and peak phytoplankton contribution to $b_{\mathrm{bp}}$ occurs at $\mathit{\xi }\approx \mathrm{3.25}$.

The fractional contributions of Fig. b 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_{\mathrm{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. estimated global absolute $b_{\mathrm{bp}}$ due to NAPs and its fractional contribution to total $b_{\mathrm{bp}}$ using analysis of correlations with Chl. Qualitatively and to first order, their global pattern of percent $b_{\mathrm{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 b here). Note that the 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_{\mathrm{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. c. These estimates are not symmetric about the $\mathit{\xi }$ 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 $\mathit{\xi }$ 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_{\mathrm{0}}$ parameter is shown in Fig. d 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. d at each $\mathit{\xi }$ 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_{\mathrm{rs}}\left(\mathit{\lambda }\right)$ data set . Both PSD parameters ($\mathit{\xi }$ and $N_{\mathrm{0}}$; Eq. ) and derived products were generated (Sect.  and ). 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. a) 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., ; ; 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_{\mathrm{0}}$ PSD parameter (Eq. ) is, as expected, higher in coastal, high-latitude, and eutrophic areas (indicating higher particle loads) and lower in the oligotrophic subtropical gyres (Fig. b). $N_{\mathrm{0}}$ varies over a few orders of magnitude. Note $N_{\mathrm{0}}$'s units of ${\mathrm{m}}^{-\mathrm{4}}$ (Eq. ) and that care should be taken when comparing Eq. () and $N_{\mathrm{0}}$ to other formulations of the PSD, e.g., the $k$ parameter in , as these are related but not equivalent (see also ).

Figure 3

Example operational retrievals of the PSD parameters (Eq. ) and their uncertainties, using monthly OC-CCI v5.0 $R_{\mathrm{rs}}$ data for May 2015: (a) PSD slope $\mathit{\xi }$, (b) $N_{\mathrm{0}}$ parameter (${\mathrm{m}}^{-\mathrm{4}}$ in log10 space), (c) uncertainty range for $\mathit{\xi }$, and (d) standard deviation (SD) of log10 of $N_{\mathrm{0}}$.

[Figure omitted. See PDF]

Algorithm uncertainties are provided on a per-pixel basis. The uncertainty range estimates for $\mathit{\xi }$ (Fig. c) (not necessarily symmetric about the $\mathit{\xi }$ value) indicate that the gyres are characterized by lower uncertainties than the more eutrophic areas, as can be expected from Fig. c. 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_{\mathrm{rs}}$ values and those due to the IOP inversion algorithm used are not included in Fig. c and d and in subsequent propagated errors. Uncertainties in the $N_{\mathrm{0}}$ parameter are more uniform spatially but higher in the gyres (Fig. d). 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. a) exhibit the expected lower values in the oligotrophic gyres and higher values elsewhere. Similarly to the results of , 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 or the method of scaling backscattering, and it is also higher than the range in CMIP5 model ensembles (cf. Fig. 1 in ). This putative underestimation in the gyres and overestimation in eutrophic areas suggests the need for algorithm tuning, which is discussed in Sect.  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 and still scarce. Here, an initial validation effort is undertaken using several other variables; see Sect. .

Figure 4

Example operational PSD-based retrievals of total and size-partitioned phytoplankton carbon, using monthly OC-CCI v5.0 $R_{\mathrm{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 $\mathrm{\mu }\mathrm{m}$), nanophytoplankton (2 to 20 $\mathrm{\mu }\mathrm{m}$), and microphytoplankton (20 to 50 $\mathrm{\mu }\mathrm{m}$).

[Figure omitted. See PDF]

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. (). Absolute concentrations of pico-, nano-, and microphytoplankton are illustrated for May 2015 in Fig. b, c, and d, respectively. Picophytoplankton C is mapped on the same color scale as total phyto C (Fig. a), but pico- and nanophytoplankton C maps have differing scales, illustrating that while picophytoplankton C varies over $\approx$ 3 orders of magnitude spatially and globally, nanophytoplankton C varies over $\approx$ 4–5 orders of magnitude, and microphytoplankton varies over $\approx$ 7 orders of magnitude (see also ; ). Note that empirical tuning will affect these ranges of variability; see Sect. . Fractional contributions of each of the three PSCs used here to total phyto C are illustrated in Fig. . Picophytoplankton dominate much of the open-ocean, lower-latitude, oligotrophic areas, contributing nearly 100 % of the carbon biomass there (Fig. a); nanophytoplankton (Fig. 5b) contribute up to $\approx$ 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 $\approx$ 45–50${}^{\circ }$ N (May 2015 is shown). As previously noted , 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. ). Retrievals use monthly OC-CCI v5.0 $R_{\mathrm{rs}}$ data for May 2015. The PSCs are picophytoplankton (a), nanophytoplankton (b), and microphytoplankton (c).

[Figure omitted. See PDF]

The fractional carbon-based PSCs (Fig. ) are ratios of two integrals of Eq. (); thus they are analytical functions of the PSD slope $\mathit{\xi }$ 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. , together with the satellite-observed $\mathit{\xi }$ 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 $\mathit{\xi }$. 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.

[Figure omitted. See PDF]

As an illustration of uncertainty propagation to derived products, the propagated uncertainty in total phyto C (Fig. a) and fractional picophytoplankton C biomass (Fig. b) is shown. Comparison of Fig. a with Fig. a 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  and ) 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 $\mathit{\xi }$ 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 $\approx \mathrm{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.  (namely , and ) or the or methods; POC products (e.g., ) 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.

[Figure omitted. See PDF]

The formulation of the PSD algorithm allows for both POC and Chl (Eq. ) to be estimated from the retrieved PSD. Due to the assumptions used, POC is phyto C multiplied by 3 (Fig. c). 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 ). 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 ). 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., ). 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 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 . An important benefit of POC is that it is a widely observed variable, available for global validation efforts (Sect. ), as is Chl. Similarly to POC, there are benefits of the PSD-derived estimate of Chl (Fig. d); 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. ). 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. ). 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., ). The values in parentheses are the standard deviations of the slope and $y$ intercept, respectively. (b) Same as in panel (a) for the $N_{\mathrm{0}}$ parameter (Eq. ) (axes in log10 space).

[Figure omitted. See PDF]

3.3

Algorithm validation, comparisons, and empirical tuning for the $N_{\mathrm{0}}$ 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) 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 $\mathit{\xi }$ (Fig. a) 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_{\mathrm{0}}$ parameter (Fig. b) is statistically significant (somewhat higher $R^{\mathrm{2}}$ than the $\mathit{\xi }$ regression) but also quite noisy. Strong clustering of the in situ observations around ${\mathrm{10}}^{\mathrm{15.5}}$–${\mathrm{10}}^{\mathrm{16.0}}$ ${\mathrm{m}}^{-\mathrm{4}}$ is observed, and the majority of these observations are somewhat overestimated in the satellite retrievals, which cluster around ${\mathrm{10}}^{\mathrm{16.25}}$. Notably, a separate cluster associated with low Chl (and low values of both satellite and in situ $N_{\mathrm{0}}$) shows that the satellite retrievals substantially underestimate these field data. Since $N_{\mathrm{0}}$ 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. ); 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_{\mathrm{0}}$ 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_{\mathrm{0}}$ 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_{\mathrm{0}}$, as done in TK16 , 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: 7$$N_{\mathrm{0}}\mathit{\_}\mathrm{tuned}={\mathrm{10}}^{\mathrm{0.3859}\log \mathrm{10}\left(N_{\mathrm{0}}\right)+\mathrm{9.5531}},$$ where $N_{\mathrm{0}}$ is the original (untuned) PSD parameter. This tuning changes $N_{\mathrm{0}}$ retrievals in a similar fashion to the TK16 tuning and is consistent with a tuning suggested by the $N_{\mathrm{0}}$ in situ validation presented here (Fig. b); namely, low satellite $N_{\mathrm{0}}$ values are increased, and high $N_{\mathrm{0}}$ values are decreased, narrowing the overall range of variability of retrieved $N_{\mathrm{0}}$ 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. ).

The tuned $N_{\mathrm{0}}$ 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_{\mathrm{0}}$ 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 and with the algorithm in Supplement Fig. S8. The histograms of these four images are compared in Fig. . The tuned retrievals are similar to those of , whereas the original retrievals are similar to those of , and the latter two have wider ranges globally compared to the former two. Of these algorithms, the simplest is that of , as it is a direct scaling of $b_{\mathrm{bp}}$, and it is based on in situ chemical analytical measurements of phyto C . 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. a). 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.

[Figure omitted. See PDF]

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. a) and the tuned algorithm (Fig. c). Both regressions have satisfactory $R^{\mathrm{2}}$ 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 $\approx$ 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_{\mathrm{0}}$ 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 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. ) for the original (Fig. b) and the tuned algorithm (Fig. d). Both comparisons exhibit very high $R^{\mathrm{2}}$ 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 $<\mathrm{0.1}$ $\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$), 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_{\mathrm{0}}$ parameter (Eq. ) for panels (c) and (d). The tuning procedure applies a linear correction to $N_{\mathrm{0}}$ 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).

[Figure omitted. See PDF]

Validation against in situ picophytoplankton carbon is presented in Fig. a (with no $N_{\mathrm{0}}$ tuning applied) and in Fig. c 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. b and d. As with the POC match-ups (Fig. b 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. a and c). The validation regression without tuning is statistically significant ($p<\mathrm{0.05}$), albeit noisy (low $R^{\mathrm{2}}=\mathrm{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=\mathrm{0.05}$ level). These validation results are generally consistent with the results of , where the tuned version of the TK16 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 (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_{\mathrm{0}}$ parameter (Eq. ) for panels (c) and (d). The tuning procedure applies a linear correction to $N_{\mathrm{0}}$ 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).

[Figure omitted. See PDF]

3.4 Further discussion, summary, and conclusions

The novel PSD–phyto C algorithm described here represents a major overhaul of the KSM09 algorithm (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 . This model formulation allows assessment of the percent contribution of phytoplankton and NAPs to total $b_{\mathrm{bp}}$ as well as the estimation of Chl from the retrieved PSD. Underlying $b_{\mathrm{bp}}$ forward modeling is hyperspectral, facilitating adaptation of the algorithm to upcoming hyperspectral sensors like PACE . PSD retrieval is achieved via spectral angle mapping (SAM), and no spectral shape is imposed on $b_{\mathrm{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_{\mathrm{bp}}$.

The algorithm has been used to create an accompanying data set based on the OC-CCI v5.0 data set (; see “Data availability”). We emphasize that the PSD parameters and derived retrievals presented here and in the accompanying data set 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 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_{\mathrm{bp}}\left(\mathit{\lambda }\right)$ is key for the PSD–phyto C algorithm, as the spectral shape of $b_{\mathrm{bp}}$ is what the PSD slope retrieval is based upon (Eq. ). The IOP algorithm is chosen here, as in KSM09, because it allows spectral $b_{\mathrm{bp}}$ retrievals that are not constrained by a specific spectral function or parameterization of $b_{\mathrm{bp}}$ as is done, for example, in the QAA (quasi-analytical algorithm; ) and the GSM (Garver–Siegel–Maritorena) algorithm . For the wavelengths used in the PSD slope retrieval, modeled and satellite-derived $b_{\mathrm{bp}}$ spectral shapes compare well when the algorithm is used, and global patterns of the retrieved PSD parameters appear reasonable. Preliminary tests with indicate that this algorithm is not as suitable for PSD retrieval in this regard. Use of , , and other IOP algorithms will be further investigated in future development of the PSD algorithm.

An important assumption of the model is that $N_{\mathrm{0}}$ for NAPs is twice that for phytoplankton so that the phyto-C-to-POC ratio is a constant $\mathrm{1}:\mathrm{3}$. This ratio is expected to vary in the real ocean, and the value used here is a reasonable average choice (e.g., ; ; ; and references therein). employed the cell sorting and chemical analysis methods of 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 $\mathrm{1}:\mathrm{3}$ is reasonable, falling within their observed ranges; however, they do observe many higher values, particularly in the oligotrophic gyres. The ratio of phyto C to 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 and phyto C retrievals (Fig.  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 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_{\mathrm{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., ; ). 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 $\mathit{\xi }$ and $N_{\mathrm{0}}$ 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., ) and is the basis for our algorithm. Therefore, in the ocean, we expect $N_{\mathrm{0}}$ to decrease with increasing $\mathit{\xi }$ 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 $\mathrm{mg}{\mathrm{m}}^{-\mathrm{3}}$ (cf. Figs. b and b).

The picophytoplankton C data in are derived from cell counts (abundance) converted to carbon using specific conversion factors for different species/groups. Namely, 60 $\mathrm{fg}$ C per cell was used for Prochlorococcus, 154 $\mathrm{fg}$ per cell for Synechococcus, and 1319 $\mathrm{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 used here, the equivalent conversion factor is $\approx$ 53 $\mathrm{fg}$ per cell for cells of the smallest diameter within the picophytoplankton range (0.5 $\mathrm{\mu }\mathrm{m}$) and is $\approx$ 1825 $\mathrm{fg}$ per cell for the largest-diameter cells within the picophytoplankton range (2.0 $\mathrm{\mu }\mathrm{m}$), indicating first-order consistency, but not full equivalency, with the methods of .

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_{\mathrm{0}}$, 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. ) is a parameterization of real-world PSDs, and while there are theoretical underpinnings (e.g., ; ; ) and observations (e.g., ; ; see also references in ) supporting its applicability, particularly over large size ranges, real-world PSDs may deviate from the power law, especially in coastal zones (e.g., ; ; ; ). 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.,. A recent study by 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., ); 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. 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 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 ). Gamma functions may be a better choice to represent marine PSDs . 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. 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_{\mathrm{0}}$ validation (Fig. b). 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_{\mathrm{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 , , and .

observe a relationship between phyto C and $b_{\mathrm{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_{\mathrm{bp}}$ spectral shape is a second-order effect that is not easily observed in in situ observations , even though theoretical modeling demonstrates a clear link (; this study). 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_{\mathrm{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_{\mathrm{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 . 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 . Hyperspectral absorption data in particular have the potential to increase information content and allow group-specific retrievals (e.g., ; but see also ). Finally, collect more global, comprehensive in situ data sets of all relevant variables, including and especially phyto C , 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. ) 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.

Code and data availability

Code and data associated with algorithm development as well as operational application to OC-CCI v5.0 data are published on the Zenodo^® repository and are available at the following DOI: 10.5281/zenodo.6354654.

An OC-CCI v5.0-based satellite PSD–phyto C data set (monthly, 1997–2020, plus monthly and overall climatologies) has been published on the PANGAEA^® repository and is freely available in netCDF format and browse images at the following DOI: 10.1594/PANGAEA.939863.

The supplement related to this article is available online at: https://doi.org/10.5194/os-19-703-2023-supplement.

Author contributions

TSK designed the study, conducted the modeling and algorithm development and data analyses, and wrote the manuscript. SB and LRL provided the EAP model code and technical support for EAP modeling. SM helped with error propagation estimates and designed the band-shifting methodology. CEK, BJ, VMV, and SS extracted match-ups and/or provided match-up in situ and satellite data set compilations. XZ provided the coated-sphere code and technical support for it, as well as validation PSD data. HL and DSFJ provided technical assistance with IOP code testing. EK tested backscattering spectral shape sensitivity. SR provided algorithm output data. SB, LRL, XZ, SM, EK, SR, CEK, BJ, SS, HL, and DSFJ read the manuscript and provided comments/edits.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

The views and opinions expressed here are those of the authors and do not necessarily express those of NASA. Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

Funding for this project was provided by NASA grant no. 80NSSC19K0297 to Tihomir S. Kostadinov, Irina Marinov, and Stéphane Maritorena. Tihomir S. Kostadinov also acknowledges support from California State University San Marcos. This work is a contribution to the Simons Foundation project Computational Biogeochemical Modeling of Marine Ecosystems (CBIOMES, 549947, Shubha Sathyendranath). This paper is also a contribution to the ESA projects Ocean Colour Climate Change Initiative (OC-CCI) and the Biological Pump and Carbon Exchange Processes (BICEP). Additional support from the National Centre for Earth Observations (UK) is also gratefully acknowledged. Shubha Sathyendranath also acknowledges additional support from the UK's National Centre for Earth Observations. NASA grants 80NSSC17K0568 and NNX15AE67G are acknowledged for support for EXPORTS and NAAMES LISST data collection.

We acknowledge Olaf Hansen, Harish Vedantham, Marco Bellacicco, Salvatore Marullo, Irina Marinov, Ivona Cetinic, Giorgio Dall'Olmo, Emmanuel Boss, Nils Haëntjens, and David Desailly for various help/useful discussions. We acknowledge the ESA OC-CCI and BICEP project teams and contributors. We acknowledge in situ PSD validation data contributors as follows: those given in ; David Siegel and the UCSB ERI/Plumes and Blooms project team; Emmanuel Boss, Nils Haëntjens, and the NASA EXPORTS and NAAMES cruise teams; and Giorgio Dall’Olmo, Emmanuele Organelli , and the AMT26 cruise team. We acknowledge all in situ data contributors to the BICEP/POCO project compilations of POC and picophytoplankton carbon data sets. The OC-CCI reference is , and the v5.0 specific reference is . The modeling and processing are done using the sinusoidal projection (one of the projections provided by OC-CCI), whereas maps here are presented in equidistant cylindrical projection (unprojected lat/long). Erik Fields, the ESA, BEAM (Brockmann Consult GmbH), and NASA are acknowledged for the re-projection algorithm. Coastlines in maps shown here are from v2.3.7 of the GSHHS data set; see . The NOAA ETOPO1 data set (https://www.ngdc.noaa.gov/mgg/global/, last access: 23 October 2015) was used in validation for bathymetry masking. Modeling and data processing were done in MATLAB^®.

The cividis colormap used in most visualizations is a variant of the viridis colormap optimized for color vision deficiency perception and is from and the function to implement it in MATLAB^® is due to Ed Hawkins.

We acknowledge Emmanuel Boss and one anonymous reviewer for their very useful comments, which helped improve the manuscript.

Financial support

This research has been supported by the National Aeronautics and Space Administration (grant nos. 80NSSC19K0297, 80NSSC17K0568, and NNX15AE67G) and the Simons Foundation (grant no. 549947).

Review statement

This paper was edited by Aida Alvera-Azcárate and reviewed by Emmanuel Boss and one anonymous referee.