HYDROPT: An Open-Source Framework for Fast Inverse Modelling of Multi- and Hyperspectral Observations from Oceans, Coastal and Inland Waters

: Biomass estimation of multiple phytoplankton groups from remote sensing reﬂectance spectra requires inversion models that go beyond the traditional band-ratio techniques. To achieve this objective retrieval models are needed that are rooted in radiative transfer (RT) theory and exploit the full spectral information for the inversion. HydroLight numerical solutions of the radiative transfer equation are well suited to support this inversion. We present a fast and ﬂexible Python framework for forward and inverse modelling of multi- and hyperspectral observations, by further extending the formerly developed HydroLight Optimization (HYDROPT) algorithm. Computation time of the inversion is greatly reduced using polynomial interpolation of the radiative transfer solutions, while at the same time maintaining high accuracy. Additional features of HYDROPT are speciﬁcation of sensor viewing geometries, solar zenith angle and multiple optical components with distinct inherent optical properties (IOP). Uncertainty estimates and goodness-of-ﬁt metrics are simul-taneously derived for the inversion routines. The pursuit to retrieve multiple phytoplankton groups from remotely sensed observations illustrates the need for such ﬂexible retrieval algorithms that allow for the conﬁguration of IOP models characteristic for the region of interest. The updated HYDROPT framework allows for more than three components to be ﬁtted, such as multiple phytoplankton types with distinct absorption and backscatter characteristics. We showcase our model by evaluating the performance of retrievals from simulated R rs spectra to obtain estimates of 3 phytoplankton size classes in addition to CDOM and detrital matter. Moreover, we demonstrate HYDROPTs capability for the inter-comparison of retrievals using different sensor band settings including coupling to full spectral coverage, as would be needed for NASA’s PACE mission. The HYDROPT framework is now made available as an open-source Python package.


Introduction
Ocean-color remote sensing has opened up the opportunity to monitor the biological and chemical processes of the ocean on an unprecedented scale. Satellite sensors continue to provide a synoptic view of ocean biogeochemistry at high spatial and temporal resolution that would be impossible to acquire through in-situ sampling campaigns [1].
More than four decades ago the first ocean-color instrument, the Coastal Zone Color Scanner (CZCS; a list of abbreviations is given under section "Abbreviations") was launched and for the first time provided a detailed picture of phytoplankton dynamics in the upper layers of the ocean. The main product derived from CZCS observations were estimates of the primary pigment found in almost all phytoplankton, chlorophyll-a [1]. The first order variation in remote sensing reflectance (R rs ) in open ocean waters is due to chlorophyll-a [2] characterized by two absorption peaks in the blue and in the red, causing a greening of the waters with increasing phytoplankton concentrations.
However, chlorophyll-a alone does not give the full picture of phytoplankton diversity [3]. The diversity in phytoplankton is characterized by different physiological and morphological traits that affect biogeochemical processes and the ecological niches that these species inhabit [4,5]. In turn, these physiological and morphological differences between phytoplankton groups affect the optical properties of the water column and in theory could be detected in the remotely sensed signal. Exploiting this second order variation in the remotely sensed signal to derive a more detailed description of phytoplankton community structure has become a top priority in the ocean-color community [3,6].
Many physiological traits of phytoplankton, which in turn affect ecological and biogeochemical processes, are correlated with cell size and pigment composition [4,5,7]. The smallest phytoplankton such as Prochlorococcus and Synechococcus are numerically the most abundant phytoplankton in the ocean and are important contributors to global primary production [8]. The largest phytoplankton such as diatoms dominate in cold, nutrient rich waters and are important contributors to the biological pump through efficient draw down of carbon to the sea floor [9]. In addition to chlorophyll-a, Prochlorococcus contains (divinyl) chlorophyll-b, Synechococcus contains a diversity of phycobilin pigments, and diatoms contain chlorophyll-c and fucoxanthin, with which they exploit different parts of the light spectrum [7,10].
Cell size and pigment composition affect absorption and backscatter characteristics of phytoplankton [9][10][11][12]. The packaging effect is an important driver of the variability in phytoplankton absorption which is depended on phytoplankton cell size and intracellular pigment concentration [10]. Numerous studies have investigated the link between community size structure and chlorophyll-a specific absorption: the phytoplankton absorption normalized to chlorophyll-a concentration derived from HPLC-based measurements. The chlorophyll-a specific absorption, especially around the blue absorption peak (440-490 nm), decreases with increasing cell size ( Figure 1) [11,13]. Small-celled phytoplankton communities are characterized by pronounced absorption peaks whereas communities dominated by larger cells exhibit less pronounced peaks. Brewin et al. [14] empirically derived backscatter coefficients for phytoplankton communities of different dominating cell sizes. Phytoplankton communities dominated by the largest cells show a low and spectrally flat backscatter, whereas smaller celled communities exhibit elevated backscatter at shorter wavelengths and exponentially decreasing scatter with increasing wavelength (Figure 2). The differences in absorption and backscatter characteristics could in theory be exploited to detect and discriminate the size structure of phytoplankton groups from the remotely sensed signal [6].
However, efforts to exploit the full spectral reflectance to obtain phytoplankton community structure have been challenging and have led to the identification of certain phytoplankton groups without estimates of biomass [2,3,15,16]. Most of the semi-analytical retrieval models only decompose either the retrieved absorption or backscatter coefficients to obtain information on community structure and ignore the inter-dependency between the IOPs and measured R rs signal. The empirical models suffer from noise in the detected pigment concentrations of representative species and are biased to the data collected. Inversion methods based on radiative transfer principles are needed to accurately separate the contribution of various phytoplankton groups to R rs [15,17]. HydroLight [18] is the state-ofthe-art radiative transfer model for ocean-color applications; however, it is computational costly and can therefore not be used to invert R rs spectra in near real time.
The HYDROPT model [19] was developed to overcome the time intensive computations of RT simulations and at the same time sustain highly accurate R rs calculations over several wavelengths. HYDROPT is based on HydroLight RT simulations and speed and accuracy are realized through interpolation of the radiative transfer solutions. The inversion is achieved by finding the best fit to the measured R rs by means of nonlinear optimization techniques. The original HYDROPT algorithm was adapted to coastal and inland waters and was only capable of fitting three optical components to a limited set of 9 wavebands from the Medium Resolution Imaging Spectrometer (MERIS). Coastal and inland waters are characterized by high levels of absorption and scattering by CDOM and particulate matter, leaving only one component to represent the absorption and scattering budget of the phytoplankton population. To retrieve multiple phytoplankton groups, HY-DROPT needs to be updated to accommodate several optical spectral models for different phytoplankton groups in addition to CDOM and detrital matter. The increasing number of ocean-observing sensors also necessitate the need for flexible waveband configuration beyond MERIS band settings. A final impetus to this research has been the publication of new absorption values of pure water in the 250-550 nm wavelength range [20], a part of the spectrum important for the discrimination of phytoplankton groups [4,7]. Chlorophyll-a specific absorption for (a) micro-(b) nano-and (c) pico-phytoplankton. A hundred absorption spectra are randomly sampled to visualize variability within and across size classes (blue lines). Mean absorption spectra are used as for the inversion with HYDROPT (red dashed lines). Absorption coefficients and standard errors are obtained from Uitz et al. [11]. Chlorophyll-a specific backscatter for (a) micro-(b) nano-and pico-phytoplankton. A hundred backscatter spectra are randomly sampled to visualize variability within and across size classes (blue lines). Mean backscatter spectra are used for the inversion with HYDROPT (red dashed lines). Backscatter coefficients and standard errors are obtained from Table  3-database D in Brewin et al. [12].
The objective of this research is three-fold. First, we present the update to the formerly developed HYDROPT algorithm [19]. HYDROPT has been completely rewritten in Python and is made open-source [21]. The most notable update is the coupling to full spectral coverage as would be needed for hyperspectral missions like NASA's Plankton, Aerosol, Cloud, ocean Ecosystem (PACE) mission [22]. HYDROPT is sensor agnostic, allowing flexible band settings in the 400-710 nm range while accounting for the bidirectional nature of R rs [23][24][25] . Second, we test the feasibility of the framework to retrieve three phytoplankton size classes in addition to CDOM and detrital matter from simulated R rs spectra. The dataset incorporates changes in phytoplankton community structure in addition to variable mass specific IOPs for phytoplankton, detrital matter and CDOM to best represent the optical variation encountered in natural water bodies [5,11,26,27]. Third, we demonstrate the applicability of HYDROPTs flexible waveband settings to both multiand hyperspectral sensors and we compare the performance of retrievals for Sea-viewing Wide Field-of-view Sensor (SeaWiFS), Ocean and Land Color Instrument (OLCI) and the PACE ocean-color sensor.

HydroLight Simulations
HydroLight forward simulation are run to obtain the remote sensing reflectance, R rs (units: sr −1 ), as a function of IOPs, viewing geometry (θ v , φ v ) and solar zenith angle (θ s ): where λ is the wavelength, L w (units: Wm −2 sr −1 nm −1 ) the water-leaving radiance and E d (units: Wm −2 nm −1 ) the downwelling irradiance just above the water surface. The remote sensing reflectance is calculated for a viewing geometry consistent with the standard HydroLight quad layout with a 10°resolution in nadir viewing angle (θ v ) and 15°resolution in azimuthal angle (φ v ). The RT simulations are run at a 5 nm resolution between 400 and 710 nm and for solar zenith angles between 0-80°at a 10°resolution. For the forward simulations a bio-optical model is chosen that includes the IOP of water, CDOM, mineral particles and phytoplankton [18]. The goal of the forward radiative transfer simulations is to obtain R rs for a realistic set of absorption and backscatter coefficients. Details on the IOP spectral models, wind and sky conditions can be found in the supplementary information (Table A1). Fluorescence and Raman scattering are ignored. A lookup table (LUT) is constructed that contains all permutations of CDOM absorption at 440 nm (ρ CDOM : 0.005-10 m −1 ) and concentrations of minerals (ρ min : 0.01-100 g m −3 ) and phytoplankton (ρ chl : 0.01-32 mg m −3 ). Every component is varied over 10 logarithmic spaced intervals between their respective minimum and maximum values, yielding a total of a 1000 simulations per waveband for every viewing geometry.

Polynomial Forward Model
HYDROPT uses a polynomial function to describe the relationship between R rs and the IOPs of the water column, the total absorption (a) and backscatter (b b ) coefficients (units: m −1 ). Since the last version of the framework [19], the use of the total scattering coefficient has been replaced by the backscatter coefficient. Through polynomial interpolation of the RT solutions a fast and accurate forward model is constructed that takes into account the bidirectional nature of R rs . The rationale behind a polynomial description of the RT solution space is that it provides an analytical form of the relationship between R rs and IOPs, the total spectral absorption and backscatter [19]. Since polynomial expressions are easily differentiable, they provide a fast way to invert reflectance spectra in conjunction with gradient-based optimization routines. The following polynomial form is chosen to relate R rs to the IOPs: where R rs is the natural log-transformed reflectance at wavelength λ, n is the degree of the polynomial expression, c ij (λ) is the fitted coefficient of the corresponding polynomial term and a and b b are the log-transformed absorption and backscatter coefficients. The powers of the polynomial terms are indicated by i and j. The (log-transformed) reflectance values, R rs , together with the absorption and backscatter coefficients, a and b b , are obtained from HydroLight simulations. It should be noted that R rs and the coefficients c ij are also dependent on the viewing geometry and sky conditions but are omitted from Equation (2) for brevity. The polynomial coefficients c are linearly interpolated to obtain R rs for any waveband, and viewing geometry. Using Equation (2) the remote sensing reflectance can now be calculated for any combination of water constituents given that: where a and b b are the total absorption and backscatter coefficients respectively, a w and b b,w are the absorption and backscatter due to water and a * i and b * b,i are the mass specific absorption and backscatter coefficients for constituent i and ρ is the concentration of the constituent. Multiple optical constituents (m) can be used to calculate R rs such as different phytoplankton groups as long as the total absorption and backscatter values are within the same range as the IOP values used in the radiative transfer simulations.
Given the matrix C of all coefficients c ij and the matrix X of the transformed polynomial terms of a and b b , Equation (2) can be rewritten in vector notation: To prevent over-fitting, the degree of the polynomial model, n, is determined through 10-fold cross-validation. Since reflectance values differ in order of magnitude between wavebands, the root mean squared relative error (RMSRE) is used as a measure of accuracy [28,29]: where k is the number of wavebands (k = 63), R rs the reflectance given by HydroLight and R rs the reflectance predicted by our forward model (Equations (2) and (4)). The root mean squared relative error scales the errors across all wavebands. The 1-SE rule is used to select the most parsimonious model within one standard error of the best model [30].

Optimization
The estimation of IOPs and concentration of aquatic optical constituents from R rs relies on minimizing the difference between the observed (R rs ) and predicted (R rs ) reflectance spectra. The loss function to be minimized is therefore expressed as: where ρ m is the absorption or concentration of water constituent m (e.g., CDOM, non-algal particles, chlorophyll), k the number of wavebands and σ i the weight assigned to waveband i. The concentration of constituents that minimize the loss function are assumed to be the best estimate for that respective reflectance spectrum. Gradient-based optimization routines, such as Levenberg-Marquardt, use the gradient of the loss function to converge to a local or global solution. Several Python packages implement these optimization routines [31,32] and either approximate the gradient numerically or allow the user to supply an analytical expression of the gradient. Since numerical approximations are computationally more costly and time consuming, providing an analytical expression of the gradient is preferred. Here we derive the gradient of R rs w.r.t the constituent concentrations ρ: where ∇R rs is the gradient for waveband i w.r.t constituents 1 to m and ∂R rs /∂ρ m is the partial derivative of R rs w.r.t constituent m. The partial derivatives can be further decomposed as follows: and ∂R rs /∂a and ∂R rs /∂b b are the first order derivative of the forward model (Equation (2)) given by: The partial derivatives ∂a /∂a and ∂b /∂b are simply: The two terms ∂a/∂ρ m and ∂b b /∂ρ m denote the first order derivative of the absorption and backscatter model for constituent m. In the case of a linear model as in Equation (3) the derivatives reduce to a * m and b * b,m respectively. Evaluating the gradient in Equation (7) for every waveband λ and constituent ρ gives the Jacobian of our forward model to be used in gradient-based optimization routines: The Levenberg-Marquardt implementation also estimates the co-variance matrix for every retrieval [32]. Given the co-variance matrix C , the relative retrieval error for constituent i can be approximated as follows: where C i,i is the variance in the estimation of constituent i andρ i is the estimated concentration of constituent i.

IOCCG Dataset
To test the generalizability of the polynomial forward model we first validate the model against an independent dataset that assumes different optical conditions under which the simulations have been conducted. Here we assess HYDROPT's capability to accurately predict R rs with changes in the particle volume scattering function (VSF). The VSF, and thus the backscatter ratio, is extremely difficult to determine in the field and laboratory. Most studies have relied on approximations of the VSF by either interpolating the VSF between angles that can be measured in the field [33] or through theoretical calculations [34]. Because of the uncertainty in the VSF of marine particles, several assumptions are made when performing in-water radiative transfer calculations. For the HydroLight simulations in this study we assumed a Fournier-Forand (FF) phase function with a 1.4% backscatter ratio for phytoplankton as well as for minerals. The well-established IOCCG dataset [27] comprises forward simulations with different phase functions for phytoplankton and particulate matter: a FF phase function with a 1% backscatter ratio for phytoplankon and a Petzold phase function with a 1.8% backscatter ratio for minerals. To establish generalizability or our forward model, the absorption and backscatter coefficients from the IOCCG simulations are used as input to the polynomial model to predict R rs . The predicted R rs spectra by HYDROPT are validated against the IOCCG R rs values.

Hyperspectral Phytoplankton Size Class Dataset
To establish the potential of HYDROPT to retrieve multiple phytoplankton size classes from R rs , we create a hyperspectral dataset containing simulated spectra [35]. The dataset takes into account the natural variability in chlorophyll-a and mass specific IOPs for three phytoplankton size classes as well as detritus and CDOM absorption. In addition, the dataset also incorporates the co-variation between the optical components as would be encountered in the case-I waters [27].
By combining HPLC samples from open ocean waters and diagnostic pigment analysis, Uitz et al. [11] estimated chlorophyll-a specific absorption spectra for three different phytoplankton size classes: pico-(<2 µm), nano-(2-20 µm) and micro-phytoplankton (>20 µm) (Figure 1). One hundred spectra for each size class are sampled from the distribution of absorption coefficients in Uitz et al. [11] (Figure 1). These spectra are randomly selected and used for forward modelling of our dataset. For the inversion part of the exercise, we only supply the average chlorophyll-a specific absorption spectra (red dashed lines in Figure 1) to HYDROPT.
Brewin et al. [12] estimated the chlorophyll-a specific backscatter coefficients for these three phytoplankton size classes using several in-situ databases that contain HPLC measurements and particulate backscatter coefficients. By applying diagnostic pigment analysis to the HPLC measurements, biomass for each size class could be estimated and subsequently their fractional contribution to the total particulate backscatter [12]. The mass specific backscatter is described by a power-law according to: where b * bp (units: m 2 mg −1 ) is the mass specific backscatter coefficient, λ 0 is the reference wavelength at 470 nm and γ 1 is the spectral slope. The mass specific backscatter coefficient, spectral slope and variability in the estimates of these coefficients were determined by nonlinear least squares fitting and bootstrapping. Here, we model phytoplankton backscatter using the parameters described in Table 3-database D in [12]. The given 95% confidence intervals are converted to standard errors assuming normally distributed errors, yielding the hundred spectral backscatter curves in Figure 2. Estimates of spectral backscatter of pico-and nano-phytoplankton did not yield significantly different spectral backscatter curves and are therefore grouped together (Figure 2b) [12].
Several studies have derived empirical relationships between the phytoplankton community size structure and total chlorophyll-a [36,37]. These population models assume that small-celled phytoplankton dominate at low chlorophyll concentrations and large cells at high chlorophyll concentrations [38]. A third group, the nano-phytoplankton, was later added to the population model and shown to dominate at intermediate chloro-phyll concentrations [5,11,39]. Here, we parameterize a three-component phytoplankton model according to Brewin et al. [5] (Table 2; parameters below first optical depth, τ < 0.6). Noise is added to the dataset by transforming the chlorophyll concentrations associated with the three size classes (c i ) to a random variable (C i ) with a log-normal distribution: log(C i ) = N (log(c i ), 1). The chlorophyll concentration c i is varied between 0.01-30 mg m −3 .
The absorption of CDOM (a cdom ) and detrital matter (a dm ) are modelled as an exponential function: with reference wavelength λ 0 at 440 nm and a spectral slope s was randomly drawn from a normal distribution N (0.0176, 0.002) for CDOM and N (0.0123, 0.0013) for detrital matter ( Figure A1) [26]. Detrital backscatter is described by a power-law according to: For the phase function, indicated byb dm , we assume a 1.4% FF phase function and reference wavelength, λ 0 , at 550 nm. The spectral slope coefficient, γ 2 , randomly varies between −2.2 and 0.2 as a function of the chlorophyll concentration [27]. Higher spectral slopes are indicative of oligotrophic waters whereas lower slope coefficients are usually found in more eutrophic areas [40]. Variability in mass specific detrital backscatter is shown in Figure A2. Hereafter we adopt IOCCG [27] to model the co-variation between chlorophyll, CDOM and detrital matter. The IOPs are forward modelled to obtain R rs at 5 nm intervals using our polynomial approximation (Equations (2) and (4)). The hyperspectral dataset yields a total of 430 simulated R rs spectra that are used for the inversion. Retrievals with a relative error (Equation (13)) higher than 200% are discarded. All accuracy metrics for the successful retrievals are calculated on log 10 transformed variables and mean absolute error (MAE) and bias are backtransformed out of log 10 space according to Seegers et al. [41].

Ocean-Color Instruments
To demonstrate the flexibility of HYDROPTs waveband settings, we have conducted an inter-comparison of chlorophyll-a retrievals between three important ocean-color instruments: SeaWiFS, OLCI and PACE. The number of wavebands for these instruments in the 400-710 nm range are 6, 11 and 63 respectively. For the PACE sensor we have assumed equally spaced wavebands at a 5 nm resolution. Together these instruments cover the global observations of the oceans since 1997 [1] and for years to come. Of interest is the variation in retrieval accuracy that can potentially arise from the differences in band setting between these instruments. However, this exercise is not meant to be a complete analysis of the differences, since this would involve a perfect error budget of each band for each instrument (σ i in Equation (6)). Even with complete knowledge of instrument measurement errors, the accuracy of atmospheric correction is instrument-dependent, introducing errors and bias in the derived R rs [42].
The retrieval accuracy between sensors for total chlorophyll-a are evaluated using 5 metrics: slope of the linear model, bias, MAE, R 2 , and fraction of successful retrievals (f-score; see Table A2). The slope, bias and MAE statistics are projected on a 0-1 scale for ease of comparison according to: with x being the metric to be transformed. A value of one indicates a perfect score whereas a score of zero reflects a deviation of ≥100% from a perfect score.

Forward Model Validation
The speed and accuracy of HYDROPT are achieved by the adoption of polynomial approximations to the exact solutions of HydroLight RT simulations. In this section, we first assess and validate the accuracy of these approximations.
By varying the degree of the polynomial in Equation (2) a measure of accuracy is obtained as a function of the complexity of the model (Figure 3). The RMSRE accuracy metric is calculated over the R rs values of all 63 wavebands with the sun at a zenith angle of 30°. Increasing the model complexity reduces the prediction error in both the validation and training set up till the 4th degree. Further increasing the polynomial degree increases the prediction error and the variance, indicating model over-fitting. A 4th degree polynomial is chosen as the most parsimonious model (Figure 3). The expected average error across all wavebands is ≈ 1%. Figure 3. Cross-validation of the model in Equation (2) for different polynomial degrees (n). The model is fitted on R rs data at nadir with the sun at zenith angle of 30°. The accuracy metric used is the root mean squared relative error (RMSRE). The blue line and purple line show the validation and training score respectively. The 68% confidence envelope is shown. The dashed line indicates the most parsimonious model chosen in this study.

RMSRE
The average error is further decomposed to visualize the variability across wavebands. Instead of relying on the RMSRE, the relative error is calculated ([R rs − R rs ]/R rs ; withR rs the reflectance predicted by HYDROPT and R rs the simulated reflectance calculated by HydroLight). The distribution of the relative error for 8 out of a total of 63 wavebands is visualized in Figure 4. The errors are calculated for a nadir viewing direction and solar zenith angle of 30°. The relative error per waveband shows a slight variation; however, most of the errors are well below 1% consistent with Van Der Woerd and Pasterkamp [19].
The polynomial model is validated for all possible viewing geometries. Figure 5 shows the mean relative error for all HydroLight quads (10 · 24 = 240) at a nominal wavelength of 440 nm. The errors at nadir (see also Figure 4) are consistent with other sensor viewing angles ( Figure 5) with only marginal deviations. Mean relative errors vary between 0.4% and 0.6% with the lowest values found in a sun-facing direction (φ = 0°), with a slight increase to 0.6% when the sensor is directed away from the sun (φ = 180°).  To test the generazibility of our forward model, we tested the predictions of HYDROPT against the independent IOCCG dataset that contains RT simulations with different phase functions for both NAP and phytoplankton. The total absorption (a) and backscatter (b b ) coefficients of the IOCCG datasets are used as inputs to the polynomial model to predict R rs . Good agreement is reached between R rs from the IOCCG dataset and the predicted values from the polynomial model ( Figure 6). The MAE indicates a relative average error of 1-2% for the three OLCI bands compared to the IOCCG R rs (Figure 6a-c). The maximum deviation in R rs are 4% for the 442.5 and 560 nm band, and 6% for the 708.75 nm band (Figure 6d-f). Together with a negligible bias and perfect R 2 score these results indicate robust predictions by the HYDROPT forward model.

Hyperspectral Inversion
HYDROPT can use the full spectral information for the inversion of R rs to IOPs and concentrations of optical constituents. The inversion of R rs to IOPs and the comparison with the forward modelled spectra are shown in Figure 7. HYDROPT aims to minimize the difference between the forward-and inverse-modelled R rs , which results in a near perfect fit for the R rs example spectrum in Figure 7a. Since HYDROPT models R rs as a function of the total spectral absorption and backscatter (Equations (1) and (2)-given the viewing geometry and solar angle) it is not surprising that the total inverse-modelled absorption and backscatter are in close agreement with the forward modelled spectra (Figure 7b,c). Figure 7d-f illustrates how HYDROPT further decomposes the total absorption into the absorption by pico-, nano-and micro-phytoplankton respectively. Of the total absorption budget (at 440 nm), roughly 50% is accounted for by phytoplankton. Pico-phytoplankton dominate the phytoplankton absorption, followed by nano-and micro-phytoplankton. The quality of the match between forward-and inverse-modelled spectra follows the same order: the retrieved pico-phytoplankton absorption shows the best agreement to the forward modelled spectrum followed by nano-and micro-phytoplankton (Figure 7d-f).
Similarly, the spectral backscatter of the different optical components can be disentangled from the total backscatter budget. As such, HYDROPT can be used as a diagnostic tool to investigate the contribution of the individual optical components to R rs .

Retrieval of Phytoplankton Size Classes
In the previous section one hyperspectral R rs spectrum was inverted to illustrate the decomposition into the spectral absorption by three phytoplankton size classes. Here we show the retrieval results for all 430 simulated R rs spectra that cover a diverse range of optical conditions and phytoplankton community composition. Pico-phytoplankton concentrations ranged from approximately 4 · 10 −3 to 2 mg m −3 (Figure 8a). Of the 430 simulated spectra, HYDROPT was able to successfully invert 78% of the spectra to chlorophyll-a concentrations contained within the pico-phytoplankton size class. The average relative error expected for pico-phytoplankton is 120% (MAE = 2.2), taken over the entire concentration range. These deviations from ground truth values are also reflected in the negative R 2 and bias of 1.81, indicating that on average HYDROPT overestimates the concentration of pico-phytoplankton by roughly 80%.  Table A2. The f-score indicates the fraction of reflectance spectra that could be successfully inverted. Data point density is indicated by color (yellow = high, blue = low).
The concentration in nano-phytoplankton varied between 4 · 10 −4 and 8 mg m −3 . Retrieval results for nano-phytoplankton were consistent with the observed concentrations (Figure 8b). HYDROPT was able to retrieve concentrations for 90% of the spectra with an R 2 statistic of 0.93 and the MAE indicating an expected relative error of 34% over this concentration range. An average overestimation is observed of 18%.
Micro-phytoplankton retrievals were highly depended on the concentration (Figure 8c). At high chlorophyll concentrations (i.e., >10 mg m −3 ) , retrievals reflected the ground truth observations more closely. However, with decreasing concentration the accuracy of the retrievals deteriorated leading HYDROPT to overestimate micro-phytoplankton more consistently (by 38%, bias = 1.38). Close to 76% of the spectra could be inverted to obtain estimates of the concentration of micro-phytoplankton with an R 2 of 0.82 and MAE of 125%.
The retrievals for total chlorophyll concentration, the sum of chlorophyll contained in the three size classes, are shown in Figure 8d. The general trend in total chlorophyll is captured well by HYDROPT with an R 2 of 0.97 and MAE of 30%. As with the retrievals for micro-phytoplankton, estimates of total chlorophyll are closer to the 1:1 line at concentrations above 10 mg m −3 . HYDROPT shows no bias in the retrievals of total chlorophyll-a (bias = 1).
The absorption and backscatter coefficients for detrital matter and CDOM could be successfully retrieved for most of the spectra (f-score 1 and 0.993 respectively) ( Figure 9). Retrieval of the absorption coefficients for detrital matter showed a relatively large spread around the 1:1 line, which is reflected in an R 2 value of 0.67 and a MAE of 121%. Detrital backscatter and CDOM absorption both yielded a R 2 of 0.98. The MAE for detrital matter and CDOM were 13% and 23% respectively. Detrital backscatter and CDOM absorption retrievals were in good agreement with the observed values -only at higher values a slight but consistent overestimation is observed (bias of 6% and 14% respectively).  Table A2. The f-score indicates the fraction of reflectance spectra that could be successfully inverted. Data point density is indicated by color (yellow = high, blue = low).

Comparison of Multi-and Hyperspectral Retrievals
The HYDROPT framework allows for the comparison of retrievals using different band settings. Here we assess the performance of the inversion using band settings of three different sensors to retrieve the total chlorophyll-a concentration contained in pico-, nano, and micro-phytoplankton ( Figure 10). The multispectral SeaWifs and OLCI and hyperspectral PACE sensor are compared to emphasize how band placement and number of spectral bands potentially influence the capability of these sensors to retrieve phytoplankton size class information. It should be kept in mind that differences in instrument measurement errors and atmospheric correction are not evaluated. Figure 10. Inter-comparison of retrieval statistics for total chlorophyll-a (sum chlorophyll contained in pico-, nano-and micro-phytoplankton) between sensors. Mean bias (MB), mean absolute error (MAE) and slope are projected on a [0-1] interval according to Equation (17). A transformed statistic of 1 indicates a perfect score whereas a score of 0 indicates ≥100% deviation from a perfect score. For calculation of statistics see Table A2. Three sensors are compared: the hyperspectral PACE sensor and multispectral OLCI and SeaWiFs sensors.

MB MAE
The hyperspectral PACE instrument outperformed the two other sensors on all the 5 accuracy statistics. Out of the 5 accuracy statistics, MAE shows the largest variation among the sensors (Figure 10). For the PACE sensor we observe a reduction in MAE performance of roughly 30% (i.e., an average relative error of 30% compared to ground truth), followed by 44% for OLCI and 73% for SeaWiFS. Although PACE will operate with a vastly greater number of wavebands (k ≈ 63), the accuracy for the other 4 metrics (R 2 , f-score, slope and mean bias) closely resembles that for the retrievals with band settings adopted by OLCI (k = 11). This indicates that with a limited number of bands, OLCI (and heritage sensors) achieves optimal waveband placement to retrieve total chlorophyll-a from R rs for water with different phytoplankton community structure. Even with 6 bands placed at SeaWifs waveband center it is possible to retrieve size class information, albeit with further reduction in accuracy and in absence of errors in sensor measurements and atmospheric correction.

Discussion
Previous attempts have been made to derive phytoplankton community structure from optical measurements, but these approaches either relied on empirical models [2,5,[43][44][45][46] or the decomposition of (remotely estimated) a ph or b bp [9,13,14,39,47,48]. The advantage of empirical models are speed and ease of use [49]. In addition, empirical models, such as abundance-based models, also only require remotely retrieved chlorophyll estimates to derive functional types [43] and/or size classes [5,46]. Other empirical models directly derive community structure from R rs [44,45]. However, empirical models require reparameterization for different biogeographic regions [5,[50][51][52] and are biased toward the training data used to build the model.
Retrieval models that adopt absorption and scattering-based approaches usually require less empiricism than the previously discussed models and adopt radiative transfer principles [3,6]. However, most methods use only part of the spectral information for the inversion and only exploit knowledge of the relationship between cell size and either the absorption or backscatter coefficients. The inter-dependency between R rs and the absorption and backscatter budget is not fully captured.
One of the few attempts to directly retrieve community composition from R rs using radiative transfer principles is Werdell et al. [15]. Werdell et al. [15] used the Gordon approximation for R rs for waters with phytoplankton community consisting of two different species. The IOPs are forward modelled and the best fit to the measured R rs is determined through nonlinear optimization (Levenberg-Marquardt). In this case, full spectral information is used for the inversion as well as knowledge of both the absorption and backscatter characteristics of the optical components. However, using their method, Werdell et al. [15] note that further decomposing the contribution of each phytoplankton group to the total IOP budget significantly degrades the accuracy of the retrievals. Their effort serves as an example of the difficulty of directly retrieving multiple phytoplankton groups from R rs , without relying on empirical methods.
Here we applied the updated HYDROPT inversion framework to R rs spectra and showed that under ideal circumstances it is possible to directly retrieve the concentration of multiple phytoplankton groups in addition to CDOM and detrital matter (Figures 8 and 9). Possible explanations for the discrepancy between our results and those of Werdell et al. [15] are (1) chlorophyll-a specific IOPs used by Werdell et al. [15] are too similar to estimate the separate contribution of the two phytoplankton groups to the total IOP budget, and/or (2) in contrast to Werdell et al. [15], we incorporated different backscatter spectra for the phytoplankton groups ( Figure 2). So detailed knowledge of both the absorption and backscatter spectra could potentially aid in improved retrievals of different phytoplankton groups from R rs .
To test the feasibility of HYDROPT to invert multi-and hyperspectral reflectance measurements to obtain phytoplankton size structure information we created a dataset with simulated R rs spectra. In this research we have shown that it is possible to obtain concentration estimates of the three phytoplankton size classes from reflectance spectra under ideal conditions (e.g., no sensor noise, perfect atmospheric correction). By supplying detailed absorption and backscatter spectra of 3 phytoplankton size classes, CDOM and detrital matter (Figures 1, 2, A1 and A2), HYDROPT can optimize the spectral fit and estimate the contribution of the optical constituents to R rs . HYDROPT additionally calculates error statistics for the retrievals. A retrieval was considered valid when the relative error is within 200%. Nano-phytoplankton could be retrieved for 90% of the spectra, followed by pico-phytoplankton (78%) and micro-phytoplankton (76%) (Figure 8). Ambiguity in retrievals could be due to several factors, including chlorophyll biomass [17] but also the absorption and scattering properties [34]. However, future research on the sensitivity of R rs to changes in phytoplankton community composition should elucidate the limits on the retrieval accuracy using semi-analytical approaches such as HYDROPT.
Many empirical and semi-analytical retrieval algorithms exist [3,6]. However, most algorithms are configured for specific wavebands making it difficult to apply these algo-rithms to different sensors to derive the same bio-optical parameters-hence obstructing data fusion and inter-comparison experiments [53,54]. The retrieval inter-comparison experiment ( Figure 10) shows how HYDROPT can easily be applied to different ocean-color sensors to invert R rs to the IOPs and concentration of optical components of interest. The retrieval performance of OLCI, and to a lesser extent SeaWiFS show that with a limited number of wavebands and correct band placement total chlorophyll-a and phytoplankton size structure can be obtained from remotely measured R rs . However, given the number of wavebands and thus the available degrees of freedom for the inversion, hyperspectral measurements as would be obtained from PACE allow many more parameters to be estimated from R rs (e.g., larger set of phytoplankton groups, aerosols, atmospheric composition etc. [22]) or target distinct optical characteristics such as absorption features of accessory pigments to distinguish different phytoplankton taxa [4,7].
For the accurate inversion of R rs spectra, HYDROPT relies on (1) a forward model rooted in radiative transfer theory, (2) specification of IOP spectral models for the region of interest, and (3) fast optimization routines that find the best fit to the measured R rs spectrum. We shortly discuss these three characteristic features of the HYDROPT framework below.
In most reflectance models the bidirectional nature of R rs is accounted for by the so-called f /Q factor [23,55]. Interestingly, the f /Q factor exhibits spectral dependence due to the relative importance of molecular scattering over particle scattering and therefore changes in the total VSF with wavelength [24]. The dependency of the f /Q factor on wavelength is complex and is further influenced by the effects of CDOM absorption and number of scattering events in the water column [56,57]. HYDROPT accounts for changes in the bidirectional effects with wavelength by fitting the forward model on HydroLight simulations at every 5 nm in the visible domain (400-710 nm). Furthermore, variability in R rs due to sensor-and sun position [24,25,58] is resolved by refitting the polynomial model (Equation (2)) for different sensor viewing geometries and solar zenith angles.
Our validation benchmark shows that HYDROPT accurately predicts R rs with an average error of 1% compared to the RT simulations (Figure 3), consistent with Van Der Woerd and Pasterkamp [19]. The relative errors in R rs show a weak dependency on sensor position, ranging between 0.4-0.6% for the 440 nm band ( Figure 5). The importance of the VSF in predicting R rs has been noted in several studies, especially for phytoplankton blooms and highly scattering environments [33,56,59]. Validation against the independent IOCCG dataset shows that HYDROPT is capable of accurately predicting R rs with a maximum deviation of 4-6%, even for waters that assume a different VSF ( Figure 6).
The validation of HYDROPT shows that the forward calculations can be applied to waters for a wide range of optical conditions. Inverting R rs for these waters is achieved by finding the best spectral fit through nonlinear optimization routines in conjunction with our forward model. The updated HYDROPT framework allows the configuration of multiple optical components and various optimization routines to invert R rs such as Levenberg-Marquardt (LM) [15,19,60], Nelder-Mead [17] and simulated annealing [61]. In addition, when signal to noise (SNR) levels are well characterized for the spectral bands, the importance of the individual bands during the optimization can be accounted for by assigning appropriate weights to every waveband (σ i in Equation (6)) [19].
HYDROPT was originally developed with the aim to retrieve chlorophyll-a concentrations from optically complex waters where absorption and scattering by CDOM and sediment dominate [19]. We like to stress that HYDROPT is now also able to give a more in-depth analysis of complex coastal waters by including multiple spectral models for CDOM and sediments originating from different sources. Given the natural variability in absorption and backscatter characteristics for phytoplankton, CDOM and particulate matter, potential future improvements could be the addition of Bayesian optimization routines to the HYDROPT framework [62]. Instead of relying on average spectral IOP models for the inversion (red lines in Figures 1, 2, A1 and A2) a probabilistic approach would allow to account for variability in the mass specific IOP reference spectra (blue lines).
The R rs predicted by HYDROPT showed a higher relative error for the IOCCG data ( Figure 6) than for our own simulated spectra using HydroLight (Figure 4). Furthermore, the prediction errors of the IOCCG data are more pronounced at longer wavelengths, reaching up to +6% at 708 nm ( Figure 6). These results indicate that accounting for changes in the VSF is an important step toward improvement of highly accurate R rs predictions. Park and Ruddick [58] included an extra parameter in their forward model to account for changes in the VSF when transitioning between case-I and case-II waters. In the future, HYDROPT could be adopted in a similar way. The inclusion of shape details of the VSF into the forward model [63,64] would be a significant improvement and allow HYDROPT to accurately predict R rs for the large optical diversity encountered in the world's oceans, coastal waters, rivers and lakes.

Conclusions
In this paper, we demonstrated the ability of the HYDROPT framework to retrieve 3 phytoplankton size classes from simulated reflectance spectra using the principles of radiative transfer theory ( Figure 8). Moreover, HYDROPT can decompose the contribution of all optically active components to R rs including CDOM and detrital matter ( Figure 9). Our framework exhibits fast and robust prediction of R rs with an expected average error of ≈1% compared to HydroLight's radiative transfer solutions ( Figure 3) while accounting for bidirectional effects across wavebands ( Figure 4) and different solar-viewing geometries ( Figure 5).
The flexible configuration of wavebands makes it possible to apply the framework to past, current and future multi-and hyperspectral missions such as PACE ( Figure 10). In addition to satellite measurements, HYDROPT is also able to process in-situ above-water R rs . In our simulated hyperspectral R rs dataset we were not able to account for all the variability that can be encountered in the real world such as errors in the measurements, fluorescence, Raman scattering and imperfect atmospheric correction to name a few. We have made HYDROPT available as an open-source Python package [21] and invite the aquatic remote sensing community to apply the framework to existing [65,66] and future multi-and hyperspectral measurements [22] and benchmark the performance under various real-world scenarios.

Conflicts of Interest:
The authors declare no conflict of interest. The funder had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: Backscatter coefficient b

RT
Backscatter ratio a * , b * b chlorophyll-a-or mass specific absorption and backscatter coefficient   Figure A2. Mass specific detrital backscatter. Spectral slopes varies randomly between −0.2 and 2.2 [27]. Red dashed line indicates spectral backscatter with averaged spectral slope (γ 2 = 1) used for the inversion.  Table A2. Definition of retrieval statistics. x are the log 10 -transformed forward modelled values, y are the estimated log 10 -transformed values from the inversion. Total number of successful retrievals are indicated by n (with δ ≤ 2; see Equation (13)) and total number of R rs spectra in the dataset by n * (=430).

MAE
Mean absolute error 10ˆ fraction of successful retrievals n/n *