Algorithmic procedure for retrieving calorific contents of marine phytoplankton from space

Graphical abstract


Specifications
Environmental Science More specific subject area: Ocean Sciences Method name: Semi-analytical ocean-colour algorithm for phytoplankton calorific contents Name and reference of original method:

Method details
Calorific value of microscopic phytoplankton in the ocean is determined by their cellular macromolecular composition. The elemental composition of carbon, nitrogen and phosphorous in phytoplankton cells i.e., the stoichiometric ratio, varies across oceanographic regions [4,9] , and this variation alters the nutritional quality of phytoplankton as food to the marine grazers [5,16] . To estimate and monitor the calorific contents of marine phytoplankton on a global scale, satellite remote-sensing is invaluable. A comprehensive method to derive calorific contents from ocean colour, captured by the Earth-Observation satellites, is described here. This method was originally introduced and applied in a study by Roy [11] . Some essential components of this bio-optical method was independently developed in previous studies by Roy el al. [12,13] , and some of the steps were customised for ad-hoc applications to state-of-the-art satellite observations of ocean. The method is generic, in the sense that it is applicable to any satellite ocean-colour database, for retrieving the concentrations of phytoplankton carbohydrate, protein and lipid, and hence the total energy value of phytoplankton, on local, regional or global scales.

Method overview
The method relies on the optical fingerprints of the living phytoplankton cells in the ocean, which is captured implicitly in ocean colour by multi-spectral satellite sensors. In a nutshell, from the raw ocean-colour data, calorific contents of phytoplankton (concentrations of carbohydrate, protein, lipid) can be retrieved following four major steps as shown in Fig. 1 . Multi-spectral ocean-colour data are generally freely available from satellite repositories managed by e.g. NASA and European Space Agency. In the first step, the light-absorption spectra of living phytoplankton in the visible light (i.e., wavelength 40 0-70 0 nm) is retrieved from remote-sensing reflectances using the so-called Inherent Optical Properties (IOP) algorithms. In the second step, the absorption coefficient of the main light-absorbing pigment inside phytoplankton cells i.e., chlorophyll-a (Chl-a) is computed from the phytoplankton absorption spectra. Following this, in the third step, information on cell-size distribution within a phytoplankton community is computed from the absorption coefficients of chlorophyll. Finally, in the fourth step, the macromolecular concentrations i.e. the concentrations of carbohydrate, protein and lipid are computed from the phytoplankton size spectra. The method further allows partitioning of macromolecular concentrations into generic, user-defined, phytoplankton size classes (PSCs).

Technical details
In the following, the major steps of the method are described with technical details. All notations used in the method are described at first use, and for the benefit of the users, those are also compiled and described on Table 1 .  Step 1: Method to retrieve total phytoplankton absorption from satellite ocean colour The ocean-observing satellite sensors capture radiances that originate from the ocean surface and pass through the atmosphere. On the ocean surface, light is reflected, absorbed and scattered by the constituents of the ocean water. The inherent optical properties (IOPs) of the dissolved and suspended constituents in the ocean water determine the processes of scattering and absorption, thereby the water-leaving radiance at different wavelengths. Phytoplankton absorption coefficient is an IOP, obtaining which, from the water-leaving radiances, is not straightforward. To derive these coefficients from remote-sensing reflectance, different IOP algorithms have been developed (details of these algorithms are documented in [6] ). A suitable IOP algorithm needs to be implemented first to obtain phytoplankton absorption coefficients, which are then used in Step 2. To illustrate (as in [12] ), one can apply a semi-analytical inversion method developed by Carder et al. [1] (further described in [6] , Chapter 9). In this method the concentration of chlorophyll-a and IOPs, including phytoplankton absorption coefficients, sum of absorption coefficients of non-algal particles plus yellow substances, and backscattering coefficient of particles at any wavelength can be computed from remote-sensing reflectance spectrum (described in [6] ). The algorithm is semi-analytical, i.e, a combination of algebraic relationships and empirical relationships. This algorithm can capture the large global variability observed in specific absorption coefficients of phytoplankton, and it can be used to invert the absorption coefficients of phytoplankton at the red peak [6] . The algorithm utilises remote-sensing reflectance at four wavebands (412 nm, 443 nm, 488 nm and 547 nm), and computes specific absorption coefficients of phytoplankton in the red peak ( a * ph (676) ) (further details of this algorithm is beyond the scope of this paper but can be found in [6] ). The coefficient a * ph (676) is an essential input to the Step 2.
It should, however, be noted that alternative IOP algorithms can be used to retrieve phytoplankton absorption coefficients. The performance of the IOP algorithms are generally comparable to each other [6] . Alternative IOP algorithms that have been used include generalized ocean color inversion model for retrieving marine inherent optical properties (GIOP) [18] , and ad-hoc empirical algorithms particularly for phytoplankton absorption coefficient at 676 nm. One constraint the users should note is that, often one may encounter 'bad pixels' (e.g. due to cloud or ice effect), at which semi-analytical inversion algorithms may not work. In such scenarios, a combination of semi-analytical and empirical inversion algorithms need to be used to generate a * ph (676) maps over large area. For illustration, global maps of the derived a * ph (676) values for a month is shown in Fig. 2 a.
Step 2: Method to compute of chlorophyll-specific absorption coefficient of phytoplankton from phytoplankton absorption spectra Chlorophyll-a is considered to be the principal light-absorbing pigment within phytoplankton cells, and so the absorption properties of phytoplankton are generally reported relative to the concentration of chlorophyll-a. The optical thickness ρ c (dimensionless) of a phytoplankton cell (assuming a homogeneous sphere with chlorophyll-a as the main absorbing pigment e.g. [8,10] ), for a ray of light of wavelength λ passing through the centre of cell, can be expressed as: and m = 0 . 06 [13] , so that ρ c Based on Duysens (1956) and van de Hulst (1957) the absorption efficiency Q a (dimensionless) of a cell, which is defined as the ratio of the light absorbed by the cell to the incident light, can be expressed as: Now, the flattening effect F (ρ c (λ)) (dimensionless), defined as the ratio of the absorption coefficient of a substance in discrete particles form suspended in the medium ( a p (λ) ), to the absorption coefficient of the same substance dissolved uniformly in a medium ( a sol (λ) ), can be expressed as: So, the specific absorption coefficient of chlorophyll-a a * chl (λ) (in m 2 ( mg Chl-a ) −1 ) of intact phytoplankton cells in water can be expressed as: Although chlorophyll-a is considered as the main absorbing pigment inside a phytoplankton cell, auxiliary pigments may also contribute to the total phytoplankton absorption a ph (λ) , especially at low chlorophyll concentrations. Decoupling the absorption by auxiliary pigments (say, a u (λ) ) and that by chlorophyll-a ( a chl (λ) ), the total phytoplankton absorption a ph (λ) at any wavelength λ can be expressed as: The coefficients of phytoplankton absorption normalized to chlorophyll-a ( a * ph (λ) ), can then be obtained by dividing both sides of Eq. (5) by B , the concentration of chlorophyll-a: where, a * chl (λ) is the specific absorption of chlorophyll-a, which is eventually a function of cell diameter as in Eq. (4) . Usually the second term on the right-hand side of the Eq. (6) is ignored on the assumption that a * chl (λ) ≈ a * ph (λ) , e.g., in the red absorption peak of chlorophyll-a (at λ ≈ 676 nm). But Eq. (6) suggests that in low chlorophyll-a conditions, this assumption may overestimate the magnitude of a * chl (λ) significantly. To overcome this uncertainty, Roy et al. [13] derived the following non-linear relationship between a * chl (λ) and a * ph (λ) : with magnitude of the parameter σ c (λ) derived as: Here, a m (λ) is the maximum value that a * ph (λ) approaches to in a limiting condition, as the total absorption decreases. Using phytoplankton absorption spectra based on in situ measurements, Roy et al. [13] estimated the magnitude a m (λ) at λ = 676 nm as 0 . 0412 m 2 ( mgChl-a ) −1 . Additionally, for λ = 676 nm, a * ci (676) is taken as 0.028 m 2 ( mg Chl-a ) −1 , based on the values reported for laboratory cultures. Eq. (7) suggests that at any given wavelength (say, 676 nm) chlorophyll-a specific absorption ( a * chl ) increases with phytoplankton absorption normalized to chlorophyll-a ( a * ph ) to a saturation level, and the parameter σ c (λ) determines how slowly the saturation level is reached.
Step 3: Method to compute phytoplankton size spectra from chlorophyll-specific absorption coefficient In seawater, phytoplankton cells are assumed to follow a particle-size distribution defined by a Jung-type power law [7,14] . Under this assumption, the differential number concentration N of phytoplankton cells per unit volume of seawater with a cell diameter of D is given by: where ξ represents the exponent of the phytoplankton size spectrum and k is the constant of power law related to the abundance of the phytoplankton population. So, within the diameter range [ D min , D max ] , the concentration of phytoplankton chlorophyll-a ( B ) can be integrated as a product of the volume of the cell ( π 6 D 3 ), the intracellular chlorophyll-a concentration ( c i ), and the total number of cells N T per unit volume of seawater within that diameter range, as follows: Similarly, at a wavelength λ = 676 nm, the total absorption by chlorophyll-a due to the phytoplankton cells in the diameter range [ D min , D max ] can be obtained from Eq. (4) by integration, as follows: Therefore, the specific absorption of chlorophyll-a at λ = 676 nm, due to phytoplankton cells in the diameter range [ D min , D max ] , can be express as: with the limiting condition: The coefficients of phytoplankton absorption normalized to chlorophyll-a at 676 nm ( a * ph (676) ) is obtained from satellite-derived chlorophyll and a ph (676) . Eq. (7) is then used to convert a * ph (676) into a * chl (676) values, which are then used on the left-hand side of Eq. (12) , to numerically estimate ξ over a given range of diameter [ D min , D max ] . The magnitude of ξ is estimated iteratively within an optimisation method, where the integral in the right-hand side of Eq. (12) is calculated numerically for an initial value of ξ , which is then optimized to minimize the differences between the two sides of Eq. (12) . Corresponding to a given value of a * ph (λ) , the estimate of ξ can be obtained through an optimisation method chosen by the user. Given the deterministic relationship, the estimates of ξ would be unique and independent of the choice of the optimisation method. Therefore, the user is left to choose the optimisation method and coding platform/package. For example, Roy et al. [12] used a single-variable bounded nonlinear-function minimization method based on golden section search and parabolic interpolation implemented on MATLAB(R) [3] . For illustration, a global map of the derived ξ values for a month is shown in Fig. 2 Step 4: Method to compute of phytoplankton macromolecular concentration from phytoplankton size spectra The allometric relationship between the cellular concentration of a phytoplankton macromolecule (e.g. carbohydrate, protein or lipid, with the notation  [11] : By computing the ratio, the parameter k is removed from the expression of macromolecule-tochlorophyll ratio χ M in Eq. (14) . This way χ M is first computed from Eq. (14) with inputs of ξ from Step 3, and then M total is computed from the observed value of B total as follows: For illustration, global maps of the derived annually-averaged concentrations of phytoplankton carbohydrate, protein and lipid are shown in Fig. 2 with χ M i j and B i j , respectively, are the macromolecule-to-chlorophyll ratio and the concentration of chlorophyll B i j in the size class [ D i , D j ] , where χ M i j follows from Eq. (14) : and based on Roy et al. [12] , B i j /B total is given by: , so that, Therefore, [ M] total can also be expressed as: Any number of PSCs can be implemented within the equations above. In particular, to obtain macromolecular concentrations in three major PSCs, i.e., picoplankton, nanoplankton and microplankton, the diameter bounds need to be specified as: D 0 = 0 . 25 μm, D 1 = 2 μm, D 2 = 20 μm, and D 3 = 50 μm (e.g. [12,15,17] ). For illustration, a MATLAB script to compute macromolecular concentrations in total and those partitioned into different size classes can be found in the supplementary materials.

Computation of uncertainty
The estimates of macromolecular concentrations based on the method described above ideally should be validated with direct in situ measurements. In the absence of primary in situ measurements (which is often the case), uncertainty in the estimates can be computed as a function of the satellite inputs to the algorithm. These uncertainty levels are valuable for the users, given that the most prominent sources of uncertainties are associated with the satellite products used as inputs to the model e.g. chlorophyll-a and phytoplankton absorption coefficients. To quantify the overall uncertainty levels in the satellite-derived estimates of macromolecular concentrations, a theoretical sensitivity analysis with respect to most important input parameters can be carried out.