Sensitivity in reflectance attributed to phytoplankton cell size : forward and inverse modelling approaches

Synoptic scale knowledge of the size structure of phytoplankton communities can offer insight in to primary ecosystem diversity and biogeochemical variability from operational to the decadal scales. Accordingly, obtaining estimates of size and other phytoplankton functional type descriptors within known confidence limits from remotely sensed data has become a major objective to extend the use of ocean colour data beyond chlorophyll a retrievals. Here, a new forward and inverse modelling structure is proposed to determine information about the cell size of phytoplankton communities using Standard size distributions of two layered spheres to derive a full suite of algal inherent optical properties for a coupled radiative transfer model. This new capability allows explicit quantification of the remote sensing reflectance signal attributable to changes in phytoplankton cell size. Inversion of this model reveals regions within the parameter space where ambiguity may limit potential of inversion algorithms. Validation of the algorithm within the Benguela upwelling system using independent data shows promise for ecosystem applications and further investigation of the interaction between phytoplankton functional types and optical signals. The results here suggest that the utility of assemblage related signals in spectral reflectance is highly sensitive to algal biomass, the presence of other absorbing and scattering constituents and the resultant constituent-specific inherent optical property budget. As such, optimal methods for determining phytoplankton size from (in situ or satellite) ocean colour data will likely rely on appropriately spectrally dense and optimised sensors, well characterised measurement errors including those from atmospheric correction, and an ability to appropriately limit ambiguity within the context of regional inherent optical properties. © 2014 Optical Society of America OCIS codes: (010.4450) Oceanic optics; (010.5620) Radiative transfer; (010.1690) Color. References and links 1. S. Bernard, R. Kudela, P. J. S. Franks, W. Fennel, A. Kemp, A. Fawcett, and G. C. Pitcher, “The requirements for Forecasting Harmful Algal Blooms in the Benguela,” Large Mar. Ecosyst. 14, 281–302 (2006). 2. E. Marañón, “Inter-specific scaling of phytoplankton production and cell size in the field,” J. Plankton Res. 30(2), 157–163 (2008). #203303 $15.00 USD Received 17 Dec 2013; revised 6 Feb 2014; accepted 6 Feb 2014; published 6 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.011536 | OPTICS EXPRESS 11536 3. J. Uitz, H. Claustre, A. Morel, and S. B. Hooker, “Vertical distribution of phytoplankton communities in open ocean: An assessment based on surface chlorophyll,” J. Geophys. Res. 111, C08005 (2006). 4. S. Alvain, C. Moulin, and Y. Dandonneau, “Remote sensing of phytoplankton groups in case 1 waters from global SeaWiFS imagery,” Deep Sea Res. I 52, 1989–2004 (2005). 5. A. Ciotti and M. R. Lewis, “Assessment of the relationships between dominant cell size in natural phytoplankton communities and the spectral shape of the absorption coefficient,” Limnol. Oceanogr. 47(2), 404–417 (2002). 6. T. Konstadinov, D. A. Siegel, and S. Maritorena, “Global variability of phytoplankton functional types from space: assessment via the particle size distribution,” Biogeosciences 7, 3239–3257 (2010). 7. R. J. W. Brewin, N. J. Hardman-Mountford, S. J. Lavender, D. Raitsos, T. Hirata, J. Uitz, E. Devred, A. Bricaud, and B. Gentili, “An intercomparison of bio-optical techniques for detecting dominant phytoplankton size class from satellite remote sensing,” Remote Sens. Environ. 115(2), 325–339 (2011). 8. J. R. V. Zaneveld, “A theoretical derivation of the dependence of the remotely-sensed reflectance of the ocean on the inherent optical properties,” J. Geophys. Res. 100(C7), 13135–13142 (1995). 9. A. Morel, D. Antoine, and B. Gentili, “Bidirectional reflectance of oceanic waters: accounting for Raman emission and varying particle scattering phase function,” Appl. Opt. 41(30), 6289–6306 (2002). 10. A. Bricaud, M. Babin, A. Morel, and H. Claustre, “Variability in the chlorophyll-specific absorption coefficients of natural phytoplankton: Analysis and parameterization,” J. Geophys. Res. 100(C7), 13321–13332 (1995). 11. A. Bricaud and A. Morel, “Light attenuation and scattering by phytoplanktonic cells: a theoretical modeling,” Appl. Opt. 25(4), 571–580 (1986). 12. A. Quirantes and S. Bernard, “Light scattering by marine algae: two-layer spherical and nonspherical models,” J. Quant. Spectrosc. Radiat. Transfer, 89(1–4), 311–321 (2004). 13. Y.-H. Ahn, A. Bricaud, and A. Morel, “Light backscattering efficiency and related properties of some phytoplankters,” Deep Sea Res. 39(11–12), 1835–1855 (1992). 14. A. L. Whitmire, W. S. Pegau, L. Karp-Boss, E. Boss, and T. J. Cowles, “Spectral backscattering properties of marine phytoplankton cultures,” Opt. Express 18(14), 15073–15093 (2010). 15. W. Zhou, G. Wang, Z. Sun, W. Cao, Z. Xu, and S. Hu, “Variations in the optical scattering properties of phytoplankton cultures,” Opt. Express 20, 11189–11206 (2012). 16. J. C. Kitchen and J. R. Zaneveld, “A three-layered sphere model of the optical properties of phytoplankton,” Limnol. Oceanogr. 37(8), 1680–1690 (1992). 17. S. Bernard, T. A. Probyn, and A. Quirantes, “Simulating the optical properties of phytoplankton cells using a two-layered spherical geometry,” Biogeosci. Discuss. 6, 1–67 (2009). 18. S. Bernard, F. A. Shillington, and T. A. Probyn, “The use of equivalent size distributions of natural phytoplankton assemblages for optical modeling,” Opt. Express 15(5), 1995–2007 (2007). 19. M. Matthews and S. Bernard, “Using a two-layered sphere model to investigate the impact of gas vacuoles on the inherent optical properties of M. aeruginosa,” Biogeosciences 7, 3239–3257 (2013). 20. M. Defoin-Platel and M. Chami, “How ambiguous is the inverse problem of ocean color in coastal waters?,” J. Geophys. Res. 112, C03004 (2007). 21. A. Morel and L. Prieur, “Analysis of variations in ocean color,” Limnol. Oceanogr. 22(4), 709–722 (1977). 22. R. G. Barlow, H. Sessions, N. Siliulwane, H. Engel, S. Hooker, J. Aiken, J. Fishwick, V. Vicente, A. Morel, M. Chami, J. Ras, S. Bernard, M. Pfaff, J. W. Brown, and A. Fawcett, “2003: BENCAL Cruise Report, NASA/TM 2003-206892,” Tech. rep. (2003). 23. C. S. Yentsch, “Measurement of visible light absorption by particulate matter in the ocean,” Limnol. Oceanogr. 7, 207–217 (1962). 24. C. S. Roesler, “Theoretical and experimental approaches to improve the accuracy of particulate absorption coefficients derived from the quantitative filter technique,” Limnol. Oceanogr. 43(7), 1649–1660 (1998). 25. T. R. Parsons, Y. Maita, and C. M. Lalli, A Manual of Chemical and Biological Methods for Seawater Analysis (Pergamon, 1984). 26. J. E. Hansen and L. D. Travis, “ Light scattering in planetary atmospheres,” Space Sci. Rev. 16, 527–610 (1974). 27. C. S. Roesler and M. J. Perry, “ In situ phytoplankton absorption, fluorescence emission, and particulate backscattering spectra determined from reflectance,” J. Geophys. Res. 100(C7), 13279–13294 (1995). 28. S. Bernard, T. A. Probyn, and F. A. Shillington, “ Towards the validation of SeaWiFS in southern African waters: the effects of gelbstoff,” S. Afr. J. Mar. Sci. 19(1), 15–25 (1998). 29. C. D. Mobley and L. K. Sundman, HydroLight 5.0, Technical Documentation (Sequoia Scientific, Inc., 2008). 30. H. Dierssen, R. Kudela, and J. Ryan, “Red and black tides: Quantitative analysis of water-leaving radiance and perceived color for phytoplankton, colored dissolved organic matter, and suspended sediments,” Limnol. Oceanogr. 51(6), 2646–2659 (2006). 31. C. D. Mobley, “ Fast light calculations for ocean ecosystem and inverse models,” Opt. Express 19(20), 18927– 18944 (2011). 32. A. Morel and S. Maritorena, “Bio-optical properties of oceanic waters: A reappraisal,” J. Geophys. Res. 106, 7163–7180 (2001). 33. A. Albert and C. D. Mobley, “An analytical model for subsurface irradiance and remote sensing reflectance in deep and shallow case-2 waters,” Opt. Express 11(22), 2873–2890 (2003). #203303 $15.00 USD Received 17 Dec 2013; revised 6 Feb 2014; accepted 6 Feb 2014; published 6 May 2014 (C) 2014 OSA 19 May 2014 | Vol. 22, No. 10 | DOI:10.1364/OE.22.011536 | OPTICS EXPRESS 11537 34. I. Reda and A. Andreas, “Solar position algorithm for solar radiation applications,” Sol. Energy 76(5), 577–589 (2004). 35. P. Werdell and S. Bailey, “An improved bio-optical data set for ocean color algorithm development and satellite data product validation,” Remote Sens. Environ. 98, 122–140 (2005). 36. J. A. Nelder and R. Mead, “A Simplex Method for Function Minimization,” Comput. J. 7(4), 308–313 (1965). 37. J. O’Reilly, S. Maritorena, B. G. Mitchell, D. A. Siegel, K. L. Carder, and S. A. Garver, “Ocean color chlorophyll algorithms for SeaWiFS,” J. Geophys. Res. 103, 24937–24953 (1998). 38. G. Zibordi, F. Mélin, S. Hooker, D. D’Alimonte, and B. Holben, “An autonomous above-water system for the validation of ocean color radiance data,” IEEE Trans. Geosci. Remote Sens. 42(2), 401–415 (2004). 39. M. Matthews, S. Bernard, and L. Robertson, “An algorithm for detecting trophic status (chlorophyll-a), cyanobacterial-dominance, surface scums and floating vegetation in inland and coastal waters,” Remote Sens. Environ. 124, 637–652 (2012). 40. S. Agusti, C. M. Duarte, and J. Kalff, “Algal cell size and the maximum density and biomass of phytoplankton,” Limnol. Oceanogr. 32(4), 983–986 (1987). 41. M. J. Sauer, C. S. Roesler, P. J. Werdell, and A. Barnard, “Under the hood of satellite empirical chlorophyll a algorithms: revealing the dependencies of maximum band ratio algorithms on inherent optical properties,” Opt. Express 20(19), 20920–20933 (2012). 42. E. Rehm and C. D. Mobley, “Estimation of hyperspectral inherent optical properties from in-water radiometry: error analysis and application to in situ da


Introduction
The spectral nature of the light emerging from the world's oceans is intrinsically linked to the biogeochemical constituents of ocean waters.Ocean colour data (primarily remote sensing reflectance (R rs )) thus represents a vast resource of information for scientists through the use of schemes which relate optical measurements to constituents of interest.Developments in satellite radiometry over the last twenty years have enabled measurement of ocean colour on unprecedented spatial and temporal scales, whilst bio-optical modelling and in situ data have established a detailed understanding of the interactions between ocean constituents and the light field.Although retrieval of chlorophyll-a concentrations ([Chl a]) from ocean colour data has been achieved with reasonable success in the open ocean, confident derivation of phytoplankton functional type indices remains a major challenge, particularly in the coastal ocean.
The southern Benguela offers substantial opportunities to develop and test new bio-optical methods related to phytoplankton characteristics, thanks to the regular occurrence of high biomass blooms.These blooms, which can often be harmful [1], are of great interest to the aquaculture and fisheries industries of the region and additionally provide distinctive signals and monospecific case studies for algorithm development and testing.For Harmful Algal Bloom (HAB) monitoring, an ability to distinguish between species on a size or ideally taxonomic level is vital [1].Beyond the Benguela, phytoplankton size structure is essential for biogeochemical model development and validation towards ecosystem and carbon cycling studies [2].
Several key approaches have been adopted to address the challenge of deriving information about phytoplankton communities from ocean colour.Abundance approaches draw on ecologically understood, empirical relationships between ocean colour derived [Chl a] and abundance of various phytoplankton sizes or functional types [3].Other approaches have used various facets of the R rs signal and related these to observed phytoplankton taxonomy [4].A final set of approaches begin to use the impacts of the algal and non-algal particle size distribution and various phytoplankton characteristics on absorption and scattering e.g.[5,6].Differences in approach have complicated direct comparison of algorithm performance [7].
Semi-analytical algorithms apply radiative transfer theory to R rs inversion.However most of these do not yet include the causative effects of cell size on ocean optics.Many semi-analytical reflectance inversions algorithms are based around a simplified reflectance approximation (Eq. (1), [8]).
Where R rs is the remote sensing reflectance (sr −1 ), λ denotes wavelength dependence, L u (0 + ) (μW cm −2 nm −1 sr −1 ) is the upwelling radiance leaving the water surface, E d (0 + ) (μW cm −2 nm −1 ) is the downwelling irradiance just above the water surface, f /Q represents the bi-directional character of the upwelling radiance [9], b b is the total backscattering coefficient and a is the total absorption coefficient.In semi-analytical inversion algorithms, R rs is parameterised with respect to inherent optical properties (IOPs) including the total absorption (a) and backscattering (b b ) coefficients which can be further parameterised using additive models of ocean constituents typically including pure water, algal particles, their detrital break down products (chromophoric dissolved organic matter (CDOM) or gelbstoff), suspended sediments and other small particles (e.g.bacteria).Phytoplankton absorption has historically been expressed as a function of [Chl a], using a wavelength specific phytoplankton absorption coefficient (a * ph ) [5,10] However, to understand how phytoplankton characteristics influence optical signals, there is a need for IOP models which can parameterise second order variability related to broad characteristics of functional types (e.g.size, internal structures, accessory pigments), in a coupled, causative way.
Size related influences on algal absorption have been realistically simulated using Mie theory [11,12].However, algal backscattering remains comparatively poorly understood and simulated [12,13,14,15].Defining parameters which accurately reflect both the diversity of phytoplankton communities and represent the effect of this on optical properties is difficult for two main reasons.Firstly, the particle size distribution can be complex and difficult to represent simplistically with one variable.Secondly, phytoplankton are not homogenous spheres as defined in Mie theory.Studies show that this assumption may be too simplistic to account for the influences of internal cellular structure on the light field [12,16,17].Standard size distributions, summarised through an effective diameter (D e f f ), can be used to address the first issue [18], whilst the use of multi-layered spheres with different refractive indices has shown promise in simulating the optical influence of intracellular structures [12,16,17,19].
Ambiguity can be a significant problem in ocean colour inversion [20], in that various combinations of ocean constituents (and their IOPs) can result in indistinguishable R rs spectra.The varying proportion and co-variance of different IOPs across the global ocean has lead to the definition of "case 1" and "case 2" waters [21].Coastal ocean (typically "case 2") waters, have significantly different IOP budgets to open ocean ("case 1") waters and in these instances multiple competing sources of optical variability, i.e. gelbstoff/non-algal particulates can enhance the ambiguity in inversion techniques [20].A coupled forward model/inversion algorithm framework offers the ability to quantity the extent of these problems in a variety of optical water types, as the sensitivity of parameter retrieval can be investigated through inversion of forward model simulated [20] and/or well resolved in situ data.To develop an operationally usable product using these principles requires consideration of computational efficiency, error estimates and a comprehensive in situ validation.Presented here is a new algorithm framework -the Equivalent Algal Population (EAP) inversion algorithm; incorporating inherent optical properties derived from a two-layered sphere model of Standard size distributions of algal cells.
An outstanding question in the field of ocean colour research is: how much of the variance in ocean colour signals is attributable to changes in various characteristics of the phytoplankton assemblage?To address this question with regards to cell size, the EAP approach is firstly used in a forward mode to generate a set of simulated data to encompass a range of optical water types from phytoplankton dominated "case one" to "case two" where substances other than phytoplankton are dominant.This simulated data is then evaluated to establish where spectral sensitivities in remote sensing reflectance (R rs ) can be attributed to cell size related effects and how this sensitivity varies in different optical water types, under competing sources of bio-optical variability (i.e.biomass levels, presence of gelbstoff/CDOM).The simulated data is then inverted to assess to what extent this sensitivity can be exploited in a semi-analytical inversion algorithm to retrieve accurate estimates of cell size.The likely errors under different bio-optical regimes as a result of inversion ambiguity are quantified.The EAP approach is then used for two applications to give perspective on issues of interest to the ocean colour community.Firstly, a further set of simulated data is used to evaluate the extent to which size related variability in R rs , can account for the second order variability commonly observed in empirical band ratio algorithms applied to satellite ocean colour data.Secondly, the approach is applied in the southern Benguela with an in situ validation conducted using data gathered from field campaigns in the St Helena Bay region.Finally, the future use and limitations of this approach in algorithm development and understanding how algal characteristics influence ocean colour is discussed.

Field measurements
A validation data set was compiled from sampling campaigns undertaken in the southern Benguela between 2002 and 2005, including shore based field campaigns at Lamberts Bay, Elands Bay, Saldanha Bay (for full methodological details of data collection see Bernard et al. [18]) and the Benguela Calibration (BENCAL) cruise report (Barlow et al. [22].Radiometric measurements were made using a Hyperspectral Tethered Surface Radiometer Buoy (Hyper-TSRB S/N 018) manufactured by Satlantic, Inc. Measurements of upwelling spectral radiance (L u (z) (z = 0.66 m depth) and above surface downwelling spectral irradiance (E d (0 + )) were made.Particulate absorption data from the shore-based field campaigns were measured using the quantitative filter pad technique [23,24].Samples included from the BENCAL cruise were similarly analysed [22]).Chlorophyll-a concentrations were analysed using fluorometric analysis [25].Measurements of particle size were made using two particle sizers: a 128 channel Coulter Multisizer II with either a 50 μm or 140 μm aperture.In addition to this, samples from 2005 were analysed using an analogous method for a Beckmann Coulter Z2 cell and particle counter.This methodology allows particles in the range 1 -70 μm to be detected, an adequate range within the ecological context of the southern Benguela, where larger cells typically dominate.The effective diameters (D e f f ) of the algal particle size distributions were calculated as shown in Eq. (2) [13,26] where d is the particle diameter and F(d)d(d) is the number of particles per unit volume in the size range d ± 1/2d(d).

Optical theory and forward model development
Any analytical derivation of R rs or L u (z) requires at a minimum, estimates of both total absorption and backscattering (if not the full volume scattering function).In additive models these are sub-divided in to the respective absorption and backscattering characteristics of the ocean constituents of interest, parameterised through modelled or empirically derived relationships.In the Equivalent Algal Population (EAP) inversion algorithm, four major components are considered: water, phytoplankton, combined gelbstoff and detritus, and non-algal particles.Of these four components, water and phytoplankton contribute to both absorption and backscattering, whilst combined gelbstoff and detritus and non-algal particles contribute only to absorption and backscattering respectively.
To account for absorption by gelbstoff/detritus (a gd ) in the algorithm, an exponential shape function is used to represent their combined effect (Eq. ( 3), [27]).The exponential slope factor S is defined as 0.012, a suitable value for the southern Benguela based on observed field data and modelling studies [28].
Similarly, the effect of backscattering by non-algal particles including detritus, bacterial or lithogenic particles is represent by a simple power low following a λ −1.2 relationship [27].Absorption and backscattering of seawater were accounted for using the data from a number of sources which is offered as default within Hydrolight 5.0 [29].
To couple the optical influence of second order variability such as cell size and pigment related features to estimation of R rs through IOPs, the forward model of Bernard et al., 2009 [17], is used here to generate [Chl a] and size related spectral phytoplankton absorption and backscattering vectors.This model makes two assumptions.Firstly, that complex phytoplankton communities can be feasibly represented using Standard size distributions, providing the effective diameter and total particle surface area are equivalent [18].Secondly, that a two-layered spherical geometry can reasonably represent the absorbing and backscattering properties of diverse phytoplankton populations [17,19].Diatoms and dinoflagellates typically dominate algal assemblages in the Benguela [1].Absorption spectra are typically relatively coherent across this group of species [30].As such, the EAP basis vectors are generated based on imaginary refractive indexes representative of the intracellular absorption characteristics of this group [17].Full details of both the real and imaginary refractive indices, their derivation and use in the two-layered sphere model can be found in Bernard et al., 2009 [17].Phytoplankton absorption (a φ ) and backscattering (b bφ ) are calculated from the [Chl a] and D e f f specific basis vectors as per Eqs.( 4) and ( 5) where, F*(d) is the Chl a specific size distribution (Eq.( 8)) and a * φ and b * bφ are the [Chl a] specific absorption and backscattering respectively.
The particle size distribution (F(d)) associated with F*(d) in Eqs. ( 4) and ( 5) is given by a Standard distribution (Eq.( 6)) ( [17,18,26] where ASF is an area scaling factor [18] which adjusts the magnitude of the equivalent distribution to the total projected surface area in the particle size distribution.V e f f is the effective variance which describes the width of the Standard distribution.A constant of 0.6 is used for V e f f , determined as an average of in situ measured distributions and deemed an appropriate parameterisation for multiple distributions through experiment [18].To scale this to a [Chl a] specific size distribution as required for generation of the basis vectors, the total relative particle volume of the Standard distribution ( V ) is calculated: And the Chl a specific size distribution is calculated from this: where c i is the intracellular [Chl a].The phytoplankton absorption and backscattering basis vectors derived from the two-layered sphere model therefore have several dependencies.Firstly, the spectral refractive index data, discussed in Bernard et al. [17].Secondly, as inferred by Eq. ( 8), F * d is dependent on the c i .A value of 2.5 kg m −3 was used for the diatom/dinoflagellate group here.In turn, the phytoplankton absorption and backscattering coefficients are also dependent on the [Chl a] and D e f f , as solved for via the inversion process.Two methods of radiative transfer calculation were tested; a reflectance approximation (REFA) approach and the EcoLight-S radiative transfer code [31] (ES).The reflectance approximation (Eq.( 1)) represents the bidirectional character of the upwelling light field using an f /Q parameter [9].This parameter is dependent on the inherent optical properties of surface waters, in particular the phase function and backscattering efficiency [32] and other factors influencing the angular structure of the light field (e.g.sun zenith and surface roughness).Look-up tables from Morel et al. [9] are used here.This spectral formulation is explicitly dependent on the [Chl a] and the solar zenith angle.To adapt these tables for use with the Satlantic radiometry, the original data from seven wavelengths were linearly interpolated at 5 nm resolution.Further, as the maximum [Chl a] value of 10 mg m −3 used by Morel et al. [9] is often exceeded in the Benguela, the maximum f /Q data is used at higher [Chl a].
As in-water measurements must in practice be made at some depth below the surface, a propagation scheme is required to determine water leaving radiance (L u (0 + )) and remote sensing reflectance (R rs ) for modelling and satellite validation purposes.In the REFA approach, the diffuse attenuation coefficient for upwelling radiance (K u ) is parameterised as per the formulation of Albert and Mobley [33] (Eq.( 9)), using total absorption and backscattering and the solar zenith angle (θ s ), derived from the algorithm of Reda and Andreas [34].
The final formulation for the REFA method is shown in Eq. (10), where η 2 τ is an assumed constant for the transmittance of upwelling radiance across the air-sea interface.
The above approach was altered to include EcoLight-S [31] to replace the reflectance approximation and associated parameterisations ( f /Q , K u etc).As with the REFA approach, L u (-0.66m) is simulated and used for inversion.EcoLight-S provides an up to 1000 fold increase in speed compared to Hydrolight within 10% accuracy for Photosynthetic Active Radiation (PAR) values [31], making it an optimal choice for inversion problems requiring compromise between accuracy of radiative transfer calculations and computational efficiency.Increased speed is achieved through several simplifications of the full Hydrolight implementation, including an azimuthally integrated version of the radiative transfer equation without inelastic scattering effects [31].Computational speed can be optimised in a number of ways (discussed fully by Mobley [31]), dependent on water type and accuracy requirements.

Generation of simulated data
A parameter space investigation and estimation of minimum inversion method error was made by forcing the forward model with the ranges of parameters listed in table 1.Values were selected to encompass both open ocean conditions and the extremes which may be experienced in highly turbid waters and Harmful Algal Blooms.Values that are ecologically unlikely with regard to known allometric scaling laws were removed from analysis e.g. a bloom dominated by 2 μm cells at [Chl a] of 100 mg m −3 [40].To contextualise the sensitivity to size seen here, with respect to current ocean colour algorithms and realistic co-variance of IOPs, a second set of forward model runs (henceforth referred to as FWD N ) were conducted to approximate the Nasa bio-Optical Marine Algorithm Development (NOMAD) data set [35].The NOMAD data set has been used in the generation and evaluation of many ocean colour algorithms.Hence, creating a similar data set here, allows for a preliminary investigation of the role of size in second order variability in the context of both global optical measurements and empirical band ratio algorithms.Compared to the original forward model runs above (table 1), more conservative ranges of a gd and b bs were used, where a gd = 50% (N-) and 400% (N+) of a φ (445) and b bs (550) = 0.005 (N-) and 0.01 (N+).
To solve for the four unknowns, a Nelder-Mead simplex, non-linear optimisation method was used [36].Input in all cases is L u (-0.66m), simulated from the forward model or from the Satlantic H-TSRB.The inversion approach was initialised with a constant set of initial conditions as follows: [Chl a] = 10 mg m −3 , D e f f = 5 μm, a gd (400) = 0.197 m −1 , b bs (550) = 0.0051 m −1 .A reduced convergence weighting was applied between 675 and 700 nm to account for the lack of fluorescence term.Residuals are calculated between the simulated or measured L u (-0.66m) and that modelled by the EAP approach.Modelled L u (-0.66m) from the reflectance approximation are then converted to R rs using derived K u , for ease of interpretation.EcoLight-S provides R rs as an automatic output.For the in situ validation in St Helena Bay, the inversion approach was applied to 75 casts from the Satlantic H-TSRB.Outputs were then compared to the coincident measured data from the validation set.Quantifying the range of R rs associated with a range of sizes indicates how much signal is potentially available to differentiate between these communities regardless of algorithm technique.Figures 1(a) and 1(b) shows the range of R rs , at specific wavelengths and [Chl a] associated with the modelled range of D e f f .In this case, both absorption by gelbstoff/detritus (a gd ) and backscattering by non-algal particles (b bs ) are low, which is typical in the southern Benguela, where there is little terrigenous input from river outflow.It is apparent that R rs sensitivity to size related effects, is spectrally variant and biomass dependent.As biomass levels increase (i.e.above 1-3 mg m −3 ), the effects of changes in size become most significant in the green (550 nm) when the effects of a minima in chlorophyll a and size related absorption are combined with a decrease in backscattering from pure water at these wavelengths.Smallest cells do not always produce the highest R rs as spectral shape imparted by the combination of changes in biomass and size in both the absorption and backscattering of phytoplankton create the non-linear response at specific [Chl a] and D e f f at the selected wavelengths.As [Chl a] increases, the position of the phytoplankton absorption minima with respect to that of pure water absorption shifts from shorter to longer wavelengths [30], resulting in varying locations of peak R rs .In the two-layered sphere model, the magnitude of phytoplankton backscattering is enhanced with increasing biomass and smaller cell sizes.Differences between the REFA and ES approaches shown in Figs.1(a) and 1(b) respectively, are primarily attributable to the different f /Q values used to generate each set of forward data (not shown).The f /Q values derived from the EcoLight-S show elevated values in the red (particularly around 709 nm), explaining the features observed at this wavelength in Fig. 1.

Forward model sensitivity to cell size
By looking at the R rs response to size in the forward model under "case two" like conditions, it becomes clear that there are high levels of non-uniqueness in R rs spectra when highly backscattering, non-algal particles (at least as characterised here) are present (Fig. 2).The IOP budget is dominated by b bs under these conditions, and size related effects only become apparent when [Chl a] becomes substantial and the higher absorption by smaller cells becomes dominant against non-algal backscattering.Similarly, additional absorption in the presence of high a gd , suppresses size related variability in the blue.In these "case two" type waters, particularly the highly scattering case, the differences in the two radiative transfer schemes becomes most pronounced.The R rs values generated by ES under high b bs (Fig. 2(b)), are significantly higher than those from the REFA approach.Also, different trends in size related variability with biomass are observed, particularly in the red, as a result of the propagation of the high b bs influence in to the bidirectional effects, which in the REFA approach are constrained by the f /Q parameterisation both in terms of magnitude and spectral shape.
In interpreting these results in the context of ocean colour data, it should be noted that the substantial sensitivity observed in Fig. 1 likely represents a "best case" scenario.The simplified size parameterisation may overstate the available signal (see discussion below (section 3.5).It is likely that IOP covariance is often more extreme (i.e.high biomass is often associated with case 2 coastal waters (Fig. 2) and measurement and/or atmospheric correction errors are likely to add substantial uncertainty.

Inversion errors in cell size retrievals
Comprehensive discussion and estimation of the ambiguity and minimum errors associated with the physical problem of ocean colour inversion was conducted by Defoin-Platel and Chami [20].An approximation of this approach, using the simulated data, can provide insight in to both the role of phytoplankton size in this ambiguity, and the minimum errors one may expect from the inversion approach used here prior to application to in situ data where measurement error will also factor.In addition to the Nelder-Mead simplex, several mathematical techniques including Levenberg-Marquardt optimisation and an evolutionary algorithm were explored, however the Nelder-Mead simplex provided the most consistent and optimal results for both the in situ and simulated data inversion.
Figure 3 shows the root mean squared errors in D e f f prediction over the simulated data set, for four generalised conditions: low a gd and low b bs ("case one"/the "Benguela type" waters above), high a gd and low b bs , low a gd and high b bs , and high a gd and high b bs (case two/gelbstoff and sediment influenced waters).Lowest errors occur in the context of Benguela type waters, with low a gd and low b bs .High error and substantial scatter in RMSE values across biomass levels for the high b bs scenarios suggests significant ambiguity may be introduced under highly scattering conditions.As a result of this presumably more accurate handling of bidirectional- ity, using the ES approach results in more consistent trends in error reduction with increased biomass under high b bs scenarios (Fig. 3(b)).Figure 4 summarises the errors in D e f f estimates, against errors in [Chl a] under low a gd and low b bs conditions.Figure 4 suggests that obtaining an accurate size estimate is less likely when the [Chl a] is estimated incorrectly and/or [Chl a] is low (<10 mg m −3 , red dots).

Applications: Cell size related variability in R rs in the context of satellite ocean colour products
To put these results in to the context of current ocean colour products, Fig. 5 shows an approximation of the maximum band ratio (MBR) approach used in the OC4 algorithm [37] using forward model output (ES) analogous to the data ranges within the NOMAD data.Firstly, Fig. 5 suggests that changes in size could be directly responsible for some of the scatter associated with MBRs, however the extent of this is influenced by coincident variability in gelbstoff/detritus and non-algal particles, as shown by the overlapping ranges in Fig. 5(a).Secondly, Fig. 5. a) Size dependent ranges of OC4 maximum band reflectance ratios, versus [Chl a] from the forward model using EcoLight-S to simulate the range of water types covered by the NOMAD data set.In all cases larger cells are associated with the higher R rs ratio.b) Example ranges of spectra from samples with similar (within 10%) R rs ratios to highlight the significant ambiguity within a reasonable error that could be associated with R rs measurements.
similar changes in ratio can result from changes in [Chl a], size, gelbstoff/detritus and non-algal particles independently or coincidently, hence the retrieval of variables as distinct from each other is highly ambiguous (Figs.5(a) and 5(b)).These results support those of Sauer et al. [41], suggesting that variability in a * φ (in our case, coincident with changes in size) may be obscured by a gd , particularly at lower biomass, where the majority of the size related signal occurs in the blue and MBR approaches are typically applied (Fig. 1).Sauer et al., [41] indicate little influence of total particle backscattering (b bp ) on MBR approaches in the context of NOMAD.The data here supports this, and closer examination of the IOP budgets under both this and the previous forward model scenarios, suggests that b bφ (and the associated size influence inferred by the two-layered sphere model) only becomes a dominant contributing factor (versus scattering by water) in the absence of significant contribution from non-algal particles and at [Chl a] greater than 3-100 depending on the wavelength and size of the assemblage; whilst a φ at all sizes can represent the dominant contribution to total absorption at [Chl a] as low as 1, at some wavelengths and when there is minimal influence from a gd [43] .

Applications: Validation in the southern Benguela
A first assessment of the efficacy of the non-linear optimisation approach used in the EAP inversion algorithm is its ability to closely match modelled and measured radiance spectra.Figure 6 shows the spectral correlation and root mean squared errors (RMSE) between the output modelled L u (-0.66 m), and the input H-TSRB derived L u (-0.66 m).Across the spectra there is a high degree of correlation, accompanied by low RMSE, with several notable exceptions.Firstly, lower convergence and highest RMSE are seen in the blue, between ≈ 400 and 430 nm.As the Benguela is typically dominated by high biomass, R rs values in the blue may be comparatively weak.The sensitivity study by Defoin-Platel and Chami [20] suggested that highest ambiguities in ocean colour inversion for total absorption may be found when this is the case.The algorithm also applies a relatively simple, constant parameterisation for the absorption of gelbstoff and detritus (a gd ).Although a gd is a relatively minor constituent in the St Helena Bay region, any uncertainty in accounting for this would likely manifest in this region of the Fig. 6.Correlation between H-TSRB derived (range given by blue fill) and modelled L u (-0.66 m), with grey shading denoting standard error and root mean squared errors (RMSE) represented by the dot colour (n=75) for a) REFA and b) ES methods.spectrum.Secondly, convergence decreases in the red, beyond 710 nm.As L u measurements are low, typically with relatively high noise in the red region of the spectrum, this may not be unexpected.Finally, a lower level of convergence is seen around 675-700 nm.This is due to the lack of fluorescence term in the model and low convergence weighting set in this range for the simplex.For the EcoLight-S variant, similar trends are observed with slightly higher RMSE.Overall the high levels of correlation and low RMSE and standard error suggest that modelled L u (-0.66 m) is converging well across the validation data set using the simplex method.Figure 7 compares the spectral phytoplankton absorption from the in situ validation data set against that selected by the EAP inversion approach.Again, we see high levels of correlation and low RMS and standard errors between the measured and modelled spectra, which are in this case entirely independent.Highest RMSE are observed in the blue (≈ 400 -475 nm), lending weight to the hypothesis above that high levels of ambiguity in the prediction of total absorption create larger error in this region of the spectrum.Lower correlation is also observed at red wavelengths, contributing to an explanation of the lower convergence observed in Fig. 6.The results from the EcoLight-S algorithm variant are almost identical, suggesting that the increased accuracy in calculation of modelled L u (-0.66 m) is not substantially influencing the optimisation process.Although no suitable data exists for validation of the phytoplankton backscattering for the Benguela, the high level of convergence in the inversion, combined with the accuracy of phytoplankton absorption measurements and low levels of contribution from other IOP sources (i.e.gelbstoff/ non-algal particles) in the Benguela, suggests that there is a high level of validity in the phytoplankton backscattering estimates.
Comparing [Chl a] from in situ samples to those predicted using the EAP algorithm yields r 2 values of 0.86 and 0.80 (p < 0.001) (REFA and ES respectively), suggesting high predictive capacity across a range of concentrations (0.86 to 309 mg m −3 ) (Fig. 8  The results above show that the EAP inversion approach has substantial utility for the prediction of [Chl a] associated with both lower biomass waters and Harmful Algal Blooms in the southern Benguela.Although the coupling of IOP choice to the radiative transfer methods and ultimate calculation of L u /R rs prevents application of other algorithms to this data in a fully independent way, the range successfully retrieved is substantially larger than any offered by current satellite ocean colour products (see Matthews et al., 2012 for further discussion of these products in the southern Benguela).The EAP results compare favourably with those of the Maximum Peak Height (MPH) algorithm applied in the southern Benguela [39].Though direct comparison is not possible in this instance as the MPH is applied to level 1 satellite data, it appears that the EAP is less prone to overestimation, particularly in high biomass conditions.The EAP method can also give an indication of the broad size structure of a community through the effective diameter (D e f f ) estimates, which are likely most useful for discriminating between assemblages dominated by smaller or larger cells under bloom scenarios.

Uncertainties in the IOP model and inversion method
Whilst the behaviour of the EAP has shown that some second order variability in radiometric measurements can be attributed to size, other sources of variability may cause additional scatter.The inclusion of additional parameterisations must be balanced against the potential for increased ambiguity (particularly for inverse modelling), so for simplicity, the EAP approach currently only includes basis vectors for one group of algal species.Additional basis vectors could be included in an admixture approach, and may allow for great accuracy of modelled variables, if these basis vectors impart significant optical variability.For example, whilst distinguishing between dinoflagellates and diatoms would be very useful for understanding HAB development, many of these species share very similar absorption spectra [30].However, it maybe be possible (and in terms of inversion ambiguity, suitable) to include basis vectors characterising unique pigments, such as those associated with cryptophytes and cryptophyte endosymbionts (e.g.Myrionecta rubra [17]), especially at elevated [Chl a] when these spectral signatures can become unique from those associated with other species [30].Additionally, several assumptions in the two-layered model may also introduce scatter and explain some sources of ambiguity.The use of a constant intracellular chlorophyll concentration (c i ), whilst a necessary simplification given the current model and data available, may result in under/over expression of size related, optical characteristics.The use of a variable c i , informed by observations, could be applied at the basis vector level, with no need for additional parameterisation.Secondly the assumption of constant Standard size distributions, whilst appropriate for mono-specific blooms, may be less relevant for samples where a background of smaller cells dominate.Alternative size distributions could be incorporated on the basis of known ecological context e.g. the dominance of smaller cells at lower [Chl a].Finally, the optical properties of living algae and the extent of their degradation products are known to be influenced by physiology, life stage, and response to growth conditions, including the light environment [30].Though there is currently no accommodation for this in the two-layered model, other bio-optical measurements (e.g.fluorescence products linked to physiology) and ancillary environmental data (e.g.bloom phenology, nutrient levels) could indicate when these factors may play a role in inversion error and ambiguity [20].Further constraint of the inversion method with regards to the above factors and IOP covariance may reduce the likelihood of finding local as opposed to global minima [20,42].

Radiative transfer considerations, application to satellite data and further uses for the EAP approach in ocean colour research
Assumptions made when solving or approximating solutions to the radiative transfer equation can introduce significant errors in the radiometric quantities estimated [31].However, for coupled modelling and inversion approaches, compromises in accuracy must be optimised against gains in computational efficiency.Here, using EcoLight-S to invert in situ measurements, results in a substantial increase in computational time (≈ 300 fold), versus using the REFA method.Whilst for the limited number of measurements in the in situ validation data set here, this is not unpractical, scaling up the approach for eventual application to satellite derived time series of ocean colour data requires further consideration of the accuracy gained.The sensitivity study conducted with the forward model, reveals a few significant differences between the two methods, which may become increasingly important in other water types.It has been suggested that substantial errors in radiative calculations can arise when assumptions made associated with single scattering are invalid and/or the single scattering albedo (ratio of scattering to total attenuation) is greater than 0.6 [44], as in the highly scattering cases here.Comparison of f /Q values from EcoLight-S to those from Morel et al. [9] suggest that as scattering increases, the f /Q factors used in the REFA approach, become increasingly inadequate.In highly scattering waters (e.g. the high b bs scenarios above e.g.Fig. 2), these differences are most substantial, resulting in significantly different R rs responses to size related sensitivity.The sensitivity to size remains highly ambiguous in both cases.Extension of the existing f /Q tables may further improve accuracy at higher biomass and small cell sizes in Benguela waters, whilst the use of radiative transfer models such as EcoLight-S appears particularly prudent in highly scattering waters such as those dominated by sediments.

Conclusions
The EAP approach has substantial utility for detecting the large range of chlorophyll a concentrations associated with algal blooms in the southern Benguela.Additionally, the approach allows for the derivation of size parameter (the effective diameter, D e f f ) that is causally linked to the inherent optical properties (IOPs) selected by the optimisation approach.However, substantial errors can occur in the prediction of this parameter in situations where phytoplankton absorption and scattering do not represent the dominant component of the IOP budget.The sensitivity of R rs measurements to size related changes in phytoplankton IOPs has also be investigated using the EAP approach in both forward and inverse formats.This allows for a preliminary understanding of the role that size related optical signals play in ocean colour ambiguity.The initial sensitivity study conducted here suggests a number of key points: • Size related sensitivity in R rs is highly dependent on algal biomass, as determined by the relative algal contribution to the IOP budget.
• Even under optimal Case 1 conditions, the error associated with radiometrically derived size retrievals increases substantially below [Chl a] of around 10 mg m −3 .This effect is exacerbated in waters with high non-algal scattering and in the presence of high absorption from gelbstoff and detritus.
• Radiative transfer techniques able to quantitatively account for phytoplankton spectral scattering and bidirectionality are needed to propagate assemblage related information in IOPs through to R rs , especially in optically complex waters.However, application of these techniques does not necessarily decrease ambiguity in the inverse problem.

Fig. 1 .
Fig. 1.Ranges of modelled R rs with variations in size (effective diameter) and biomass ([Chl a]), under low b bs and low a gd conditions for a) REFA and b) ES methods.Dots indicate R rs associated with smallest cells.c) Shows example ranges of spectral R rs at selected [Chl a] across the modelled size range using ES.

Fig. 2 .
Fig. 2. Ranges of modelled R rs to variations in D e f f and [Chl a], under high b bs and high a gd conditions.Note the differences in scale, where (b -ES) shows much higher R rs values than (a -REFA).Dots indicate R rs associated with smallest cells.c) Shows example ranges of spectral R rs at selected [Chl a] across the modelled size range using ES.

bFig. 3 . 3 Fig. 4 .
Fig. 3. Root mean squared errors (RMSE) for D e f f retrievals at different [Chl a] under various parameter combinations for a) the REFA and b) ES methods.

Fig. 7 .
Fig. 7. Correlation between spectrophotometer derived (range given by blue fill) and modelled a φ .Gray shading denotes standard error and root mean squared errors (RMSE) represented by the dot colour (n=49) for a) REFA and b) ES methods.

Fig. 8 .
Fig. 8. Correlations between measured and algorithm predicted a) [Chl a] (n= 73) and b) D e f f (n=44).Shaded areas show 95% confidence intervals based on linear regression for each data set.Red and blue dots represent values derived from the REFA and ES approaches respectively.r 2 values for [Chl a] estimation were 0.86 (REFA) and 0.80 (ES), and 0.45 (REFA) and 0.25 (ES)for D e f f estimates, with p<0.001 in all cases.Estimates of absolute percentage error (|ψ|) and bias (ψ) are given as per the method used in[38].

Table 1 .
Summary of parameter ranges used in forward model data simulation