Multi-Spectral Remote Sensing of Phytoplankton Pigment Absorption Properties in Cyanobacteria Bloom Waters: A Regional Example in the Western Basin of Lake Erie

Abstract

Phytoplankton pigments absorb sunlight for photosynthesis, protect the chloroplast from damage caused by excess light energy, and influence the color of the water. Some pigments act as bio-markers and are important for separation of phytoplankton functional types. Among many efforts that have been made to obtain information on phytoplankton pigments from bio-optical properties, Gaussian curves decomposed from phytoplankton absorption spectrum have been used to represent the light absorption of different pigments. We incorporated the Gaussian scheme into a semi-analytical model and obtained the Gaussian curves from remote sensing reflectance. In this study, a series of sensitivity tests were conducted to explore the potential of obtaining the Gaussian curves from multi-spectral satellite remote sensing. Results showed that the Gaussian curves can be retrieved with 35% or less mean unbiased absolute percentage differences from MEdium Resolution Imaging Spectrometer (MERIS) and Moderate Resolution Imaging Spectroradiometer (MODIS)-like sensors. Further, using Lake Erie as an example, the spatial distribution of chlorophyll a and phycocyanin concentrations were obtained from the Gaussian curves and used as metrics for the spatial extent of an intense cyanobacterial bloom occurred in Lake Erie in 2014. The seasonal variations of Gaussian absorption properties in 2011 were further obtained from MERIS imagery. This study shows that it is feasible to obtain Gaussian curves from multi-spectral satellite remote sensing data, and the obtained chlorophyll a and phycocyanin concentrations from these Gaussian peak heights demonstrated potential application to monitor harmful algal blooms (HABs) and identification of phytoplankton groups from satellite ocean color remote sensing semi-analytically.

1. Introduction

Spectral light absorption and backscattering are the two inherent optical properties directly controlling the light field in water and further influencing water color. Optical properties of phytoplankton, specifically the absorption coefficients of the pigments inside them, play a key role in determining both the penetration of radiant energy in water and the use of this radiant energy for photosynthesis. These pigment absorption coefficients and their concentrations are important for understanding photosynthetic rate [1,2], identifying and quantifying phytoplankton functional groups [3] and determining size class distributions ([4,5] and references therein). These properties of phytoplankton and their associated backscattering, along with colored dissolved organic matter absorption and non-algal particle absorption and scattering directly control the light field of water.

With the objective of expanding the ability to detect more pigments than just chlorophyll a, various methods have been proposed [6,7,8]. Most of these methods could obtain one or more pigments in addition to chlorophyll a from R_rs(λ), but the accuracy was influenced by the presence of non-algal particles and dissolved organic matters in the water column. To reduce the influence from components other than phytoplankton, Wang et al. [9] developed a multi-pigment inversion model (MuPI) to obtain information of multiple pigments from hyperspectral R_rs(λ). This model incorporated a Gaussian decomposition scheme into a semi-analytical inversion model and demonstrated that 13 Gaussian curves that contain important phytoplankton pigment information can be retrieved from hyperspectral remote sensing reflectance. The Gaussian scheme was first described by Hoepffner and Sathyendranath [10,11], and was one of the methods proposed to obtain multiple pigments from phytoplankton absorption coefficients [12,13,14,15,16,17]. Hoepffner and Sathyendranath [10] used these Gaussian curves to represent light absorption of different pigments and demonstrated that some of the Gaussian peak heights showed good relationships with phytoplankton pigments, such as peaks at 435 and 675 nm for chlorophyll a, peaks at 643 nm for chlorophyll c, peaks at 460 and 655 nm for chlorophyll b, and peaks at 489 and 582 nm for carotenoids. This Gaussian decomposition scheme sheds light on the potential of obtaining information for multiple phytoplankton pigments beyond chlorophyll a semi-analytically using bio-optical techniques.

Most of the past and current operational ocean color satellite missions are multispectral. It is thus necessary to evaluate the viability and associated uncertainties of obtaining these Gaussian curves from multi-spectral remote sensing data. Using measurements from cyanobacteria bloom waters in different regions, we validated the MuPI performance in obtaining multiple independent Gaussian curves, and the sensitivity of MuPI to the specific spectral bands of existing ocean color satellite sensors. Further, the inversion scheme was applied to Hyperspectral Imager for the Coastal Ocean (HICO), and Moderate Resolution Imaging Spectroradiometer aboard the Aqua satellite (MODIS-Aqua), and MEdium Resolution Imaging Spectrometer (MERIS) imagery over Lake Erie to obtain the spatial distributions of pigment absorption coefficients, chlorophyll a (Chl-a) and phycocyanin (PC) for a cyanobacteria bloom event in the western basin of Lake Erie.

2. Materials and Methods

2.1. In Situ Measurements

2.1.1. Study Areas

The in situ measurements were obtained from different areas at different seasons around the world, which include ponds in the state of Mississippi, U.S.A. (MS), Lake Taihu, China (LT), and Lake Erie, U.S.A. (LE).

Mississippi Pond dataset (MS): This dataset includes 41 samples of remote sensing reflectance (R_rs(λ), sr⁻¹) and absorption coefficients of phytoplankton (a_ph(λ), m⁻¹), non-algal particles (a_d(λ)), and gelbstoff (a_g(λ)). This dataset (see Mishra et al. [18] for details) was collected from a series of highly turbid and productive aquaculture ponds at various bloom stages (initiation, peak, and senescence) of cyanobacteria (Chl-a concentration varied from 59 to 1377 mg∙m⁻³), located in northwestern Mississippi, U.S.A., where cyanobacteria blooms with Planktothrix agardhii as the most abundant species were occurring when the in situ data were collected

Lake Taihu dataset (LT): This dataset includes 45 R_rs(λ) spectra, corresponding spectra of a_d(λ), a_g(λ), and a_ph(λ) (350–750 nm) collected from Lake Taihu, China under cyanobacteria bloom conditions (Chl-a: 10–222 mg·m⁻³) where the dominant species was Microcystis aeruginosa. Water samples and optical data of surface water were collected during two surveys in January–August 2011 and November 2011, respectively.

Lake Erie dataset (LE): This dataset is composed of 36 R_rs(λ), a_d(λ), a_g(λ), and a_ph(λ) samples, and 20 Chl-a and PC concentrations. Cyanobacteria blooms dominated by Microcystis aeruginosa have been seasonally reoccurring in the western basin of Lake Erie since the mid-1990s [19]. The field measurements were collected during cruises from 20–23 August 2013 and 18–21 August 2014 in western Lake Erie at 16 and 20 stations, respectively, during the peak cyanobacteria bloom period [20,21], with chlorophyll concentrations ranging from >1 to 100 s of mg·m⁻³.

Since the ponds are too small for the MODIS and MERIS 1 km spatial resolution imagery, Lake Erie was selected as a local example for application of the method to HICO, MODIS and MERIS imageries to obtain pigment concentrations and Gaussian absorption distribution. The cyanobacteria bloom which occurred in Lake Erie in August 2014 was one of the most serious, causing disruption to the drinking water supply in Toledo, U.S.A. The 2014 cruise was conducted during this event and the sampling locations along with the Chl-a and PC concentrations at each station, are shown in Figure 1, where PC concentrations as high as ~200 mg∙m⁻³ were noted in the southern portion of Lake Erie.

2.1.2. Remote Sensing Reflectance

The remote sensing reflectance of the MS dataset was acquired in the range of 400–900 nm with a sampling interval of 0.3 nm by deploying a dual sensor system with two inter-calibrated ocean optics spectroradiometers (Ocean Optics Inc., Dunedin, FL, USA). The in situ R_rs(λ) from LT and LE were measured with a hand-held ASD spectroradiometer (Analytical Spectral Device, Inc. Boulder, Colorado, U.S.A.), from 350 to 1050 nm at 1.4 nm intervals.

2.1.3. Absorption Coefficients

Surface water samples were collected and filtered immediately after the R_rs(λ) measurements and analyzed on the same day in the laboratory. The particulate absorption coefficient was quantified utilizing the transmittance or transmittance/reflectance method of Tassan and Ferrari [22] along with the NASA Ocean Optics Protocols, Revision 4, Volume IV protocol [23]. Percent transmittance and reflectance were measured for each sample filter to calculate the particulate absorption spectrum. Filters were then bleached, and the transmittance and reflectance measured to calculate the absorption coefficients of non-algal particles. The phytoplankton absorption coefficients (a_ph(λ)) at each station were calculated by subtraction of non-algal particle absorption coefficients from total particle absorption coefficients. The Gaussian peak heights (a_Gau(λ)) were obtained from a_ph(λ) using the Gaussian decomposition scheme described in Wang et al. [9], which used a least-square curve fitting technique in Matlab to obtain the a_Gau(λ). These a_Gau(λ) were used as ground truth to validate the outputs from R_rs(λ) inversions.

For the MS dataset, a Perkin Elmer Lambda 850 Spectrophotometer (Perkin Elmer Inc., Waltham, MA, USA) was used to measure the absorption coefficient of phytoplankton, detrital matter, and gelbstoff in the 380–750 nm range at 1 nm spectral resolution. Detailed information regarding the environmental characteristics and measurement methods can be found in Mishra et al. [18]. A Shimadzu UV-2401 Spectrophotometer was used for LT. More information about LT and details about water sample collection, measurement protocols, and processing methods can be found in Duan et al. [24] and Ma et al. [25]. Samples from LE were measured with a Perkin-Elmer Lambda 35 UV/Vis Dual-beam Spectrophotometer from 300–800 nm at 1 nm resolution as described in Mouw et al. [26].

2.1.4. Pigment Concentrations and Group Composition

Water samples for Chl-a and PC concentration from Lake Erie were filtered through 0.7 and 0.2 μm glass fiber filters (Whatman GF/F, 25-mm and 47-mm diameter), respectively. The Chl-a was estimated using a Turner Designs 10-AU Fluorometer in the laboratory following the NASA Ocean Optics Protocols, Revision 5, Volume V [23]. Following NASA protocols [27], filters for phycocyanin concentration determination were extracted in phosphate buffer (Ricca Chemical, Arlington, Texas, U.S.A., pH 6.8) using two freeze-thaw cycles, followed by sonication. The concentration of PC was estimated from the relative fluorescence measured using a Turner Aquafluor fluorometer [21]. Phytoplankton populations for the Lake Erie stations in 2013 were also used in this study which were identified and counted using standard light microscopy [21].

2.2. Satellite Imagery

The Hyperspectral Imager for the Coastal Ocean (HICO) was the first space-borne hyperspectral imaging spectrometer designed to sample the coastal ocean. HICO covered selected coastal regions at 90 m spatial resolution with full spectral coverage (380 to 960 nm sampled at 5.73 nm intervals). HICO imagery of the same time frame as in situ measurements in Lake Erie in 2014 was downloaded from the Oregon State University website (http://hico.coas.oregonstate.edu/), and following the steps provided on the website the level 1B (L1B) imagery was atmospherically corrected to level 2 (L2) data using the online version of Tafkaa_6s model, followed by georectification. The HICO R_rs(λ) data have about the same spectral resolution as the in situ measured reflectance (5.73 nm for HICO and 5 nm for original R_rs(λ)). These data were used as an example for MuPI application in hyperspectral satellite imagery.

The Moderate Resolution Imaging Spectroradiometer onboard the Aqua satellite (MODIS-Aqua) launched in 2002 is one of the contemporary satellite ocean color missions. The L1B imagery with 1 km spatial resolution was downloaded from the National Aeronautics and Space Administration (NASA) website (https://ladsweb.nascom.nasa.gov/). The L2 imagery (with 748 nm included in R_rs(λ)) from the same time frame as in situ measurements in Lake Erie in 2014 was obtained using the data processing L2gen from the SeaWiFS Data Analysis System (SeaDAS) software using the standard near infrared (NIR) scheme [28,29,30] for atmospheric correction. MODIS imagery was used as an example of MuPI application to multi-spectral satellite remote sensing data.

MEdium Resolution Imaging Spectrometer (MERIS) on board the European Space Agency’s Envisat platform, was launched in 2002 and ended its mission in May 2012. The L3 seasonal composed imagery and L2 reduced resolution imagery (1 km spatial resolution) was obtained from NASA ocean color website (https://oceancolor.gsfc.nasa.gov/). Lacking coincidence with in situ measurements, the same day imagery of MODIS and MERIS on 3 September 2011 was used for the comparison of the retrieving results. The MERIS seasonal composed imagery in Lake Erie was used to obtain a general seasonal variation of Gaussian absorption coefficients.

2.3. The a_Gau(λ) Estimation

2.3.1. MuPI Model

To obtain a_Gau(λ) from remote sensing reflectance, a multi-pigment inversion method [9] was used. The corresponding functions used in this method are included in

Table 1. Summary of the functions and symbols used in multi-pigment inversion model (MuPI).

Index	Function	Description	References
1	$R_{\mathrm{rs}}\left(\lambda \right)=F(\frac{b_{\mathrm{b}}\left(\lambda \right)}{a\left(\lambda \right)+b_{\mathrm{b}}\left(\lambda \right)})=\frac{0.52\left({\mathrm{G}}_{1}u+{\mathrm{G}}_2{u}^2\right)}{1-1.7\left({\mathrm{G}}_{1}u+{\mathrm{G}}_2{u}^2\right)}$ with $u=\frac{b_{\mathrm{b}}\left(\lambda \right)}{a\left(\lambda \right)+b_{\mathrm{b}}\left(\lambda \right)}$	Remote sensing reflectance as a function of a(λ) and bb(λ) G1 = 0.089 sr−1; G2 = 0.125 sr−1	[31,32,33,34,35]
2	$a\left(\lambda \right)=a_{\mathrm{ph}}\left(\lambda \right)+a_{\mathrm{dg}}\left(\lambda \right)+a_{\mathrm{w}}\left(\lambda \right)$ $b_{\mathrm{b}}\left(\lambda \right)=b_{\mathrm{bw}}\left(\lambda \right)+b_{\mathrm{bp}}\left(\lambda \right)$	Absorption and backscattering coefficients as the total of their components
	aph(λ): phytoplankton absorption coefficient; adg(λ): absorption coefficient of non-algal particles and gelbstoff; aw(λ): absorption coefficients of pure seawaters; bbp(λ): beam attenuation coefficients of suspended particles; bbw(λ): beam attenuation coefficients of water molecules
3	$b_{\mathrm{bp}}\left(\lambda \right)=b_{\mathrm{bp}}{\left(\frac{{\lambda }_0}{\lambda }\right)}^{\eta }$ with $\eta =2.0\left(1-1.2\exp \left(-0.9\frac{r_{\mathrm{rs}}\left(443\right)}{r_{\mathrm{rs}}\left(55x\right)}\right)\right)$	λ0: 440 nm	[33]
4	$a_{\mathrm{ph}}\left(\lambda \right)={\displaystyle \sum _{i=1}^n{a}_{\mathrm{Gau}}\left({\lambda }_i\right)\exp \left[-0.5{\left(\frac{\lambda -{\lambda }_i}{{\mathsf{\sigma }}_i}\right)}^2\right]}$	σi: width of the ith Gaussian curve (FWHM = 2.35 × σi, FWHM as full width at half maximum) aGau(λi): the peak height of peak centered at λi	[10]
5	$a_{\mathrm{dg}}\left(\lambda \right)=a_{\mathrm{dg}}\left({\lambda }_0\right)\exp \left(-S\left(\lambda -{\lambda }_0\right)\right)$	λ0 = 440 nm; S: spectral slope of adg(λ)	[36]
6	$\delta =\frac{\sqrt{\frac{1}{{\mathrm{N}}_{\lambda }}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\lambda }}{\left({\widehat{R}}_{\mathrm{rs}}\left({\lambda }_i\right)-R_{\mathrm{rs}}\left({\lambda }_i\right)\right)}^2}}}{\frac{1}{{\mathrm{N}}_{\lambda }}{\displaystyle \sum _{i=1}^{{\mathrm{N}}_{\lambda }}{R}_{\mathrm{rs}}\left({\lambda }_i\right)}}$	Cost function for spectral optimization Ȓrs(λ) and Rrs(λ): Estimated and measured remote sensing reflectance	[36,37,38,39]

. In this method, the Gaussian scheme, which uses 13 Gaussian curves to reconstruct the spectrum of phytoplankton absorption coefficient, was incorporated into the relationship in which R_rs(λ) is a function of total spectral absorption (a(λ)) and backscattering (b_b(λ)) coefficients [31,32]. With an R_rs(λ) spectrum as the only input, the function was solved with a spectral optimization method which returns the values of unknowns (a_Gau(λ), absorption of non-algal particles and gelbstoff: a_dg(λ), and backscattering coefficient of particles, b_bp(λ)), that minimizes the difference between the estimated and input R_rs(λ) spectra. A summarized flowchart of the method is shown in Figure 2. The parameters used for the 13 Gaussian curves in this study are shown in

Table 2. The Gaussian curve peak locations, widths (σ = FWHM/2.35) and the corresponding power-law relationships among the peak heights. The peak heights of the Gaussian curves located at 386.6, 414, 451.7, 484, 515.6, and 677 nm can be estimated from the peak height of the Gaussian curve centered at 435 nm (annotated as x₁), those of 548.8, 584.4, 653, and 693.5 nm can be estimated from the peak height at 617.6 nm (x₂).

Gaussian Bands and Pigment Relationships
Peaks	Pigment	Peak (nm)	Width (nm)	Peak Height	R2
1	Chl-a	386.6	18.8	y = 1.52x1	0.99
2	Chl-a	414	10.7	y = 0.97x1	0.998
3	Chl-a	435	12	x1	1
4	Chl-c	451.7	18.5	y = 0.90x1	0.995
5	Carot	484	19.6	y = 0.95x1	0.99
6	Carot	515.6	18	y = 0.53x1	0.99
7	PE	548.8	15.7	y = 0.76x20.92	0.99
8	Chl-c	584.4	17	y = 90x20.94	0.997
9	PC	617.6	16	x2	1
10	Chl-c	636	11.6	y = 0.35x21.1	0.99
11	Chl-b	653	14	y = 0.82x20.87	0.99
12	Chl-a	677	10.6	y = 0.69x1	0.99
13	others	693.5	20	y = 0.37x20.92	0.99

, which includes the peak locations, widths and empirical relationships between the peak heights.

2.3.2. Gaussian Parameters

As with other semi-analytical algorithms (e.g., HOPE [38]; GSM [40], GIOP [36]), a_Gau(λ) in MuPI were retrieved using a spectral optimization method. The Gaussian scheme used in MuPI plays a critical role in a_Gau(λ) estimation, and the Gaussian peak heights, locations, and widths are very important factors in the Gaussian scheme. To determine the peak locations and widths of the Gaussian curves, the peaks and troughs of a_ph(λ) spectra were used as the initial inputs in the optimization process. As shown in Table 2 and Figure 3, a total of 13 Gaussian curves, corresponding to peaks or shoulders of measured a_ph(λ) spectra, were used to represent the absorption of phytoplankton pigments at different wavelengths: chlorophyll a: 386.6, 414, 435, 677, and 693.5 nm; chlorophyll b: 653 nm; chlorophyll c: 451.7, 584.4, and 636 nm; carotenoids: 484 and 515.6 nm; phycoerythrin: 548.8 nm; and phycocyanin: 617.6 nm. These peak locations and widths (the second and third columns of Table 2) were potentially influenced by pigment composition. Upon further analysis, we found that using fixed values for peak locations and widths in Table 2 cause no more than 5% differences in a_ph(λ) estimation versus when these peak locations and widths were set as free variables. Comparatively, the peak height of Gaussian curves influenced by pigment concentrations play a critical role in a_ph(λ) estimation. Thus, in the MuPI method, the peak locations and widths presented in Table 2 were fixed as constants, and the peak heights of these 13 Gaussian curves, containing bio-optical information of phytoplankton pigments, are the parameters to be obtained from remote sensing reflectance.

It is impractical to retrieve 13 peaks simultaneously from an R_rs(λ) spectrum via spectral optimization because more variables in the function usually means much higher uncertainties in the outputs. In determining the number of independent a_Gau(λ) to be retrieved from the MuPI model, a series of tests were conducted, and the mean, median, maximum and minimum unbiased absolute percentage difference (UAPDs, Equation (1)) were compared for different datasets when the number of independent a_Gau(λ) vary among 13, 5 (Peak 2, 5, 7, 9 and 12), 3 (Peak 3, 5 and 9) and 2 (Peak 3 and Peak 9). The combinations of independent a_Gau(λ) were determined based on the highest correlation coefficients (R²) with other a_Gau(λ). As expected, an obvious decrease in the percentage differences are noticed when the number of independent a_Gau(λ) are reduced from 13 to 2 for all three in situ datasets, as shown in

Table 3. The mean (Mea.), median (Med.), maximum (Max.) and Minimum (Min.) UAPD of different dataset with different number of unknowns (N) to be retrieved from MuPI.

Unknowns (N)	Datasets	UAPD (%)
		Mea.	Med.	Max.	Min.
13	All	47	43	73	24
5	All	37	38	46	27
3	All	33	34	40	25
2	All	32	34	38	24
13	MS	38	41	81	13
5	MS	29	28	50	15
3	MS	25	29	35	15
2	MS	26	28	37	16
13	LT	42	36	71	23
5	LT	37	38	47	25
3	LT	35	34	47	24
2	LT	32	31	49	22
13	LE	64	60	111	39
5	LE	45	45	56	36
3	LE	41	40	49	33
2	LE	39	39	50	33
13	LE-2014	62	65	108	30
5	LE-2014	31	28	45	17
3	LE-2014	27	28	35	16
2	LE-2014	22	21	33	16

, although the retrieval accuracy varies from dataset to dataset because of the different contributions of phytoplankton to the total absorption (100% × a_ph(λ)/a(λ), Figure 4). Thus the height of two peaks, here a_Gau(435) and a_Gau(617.6), instead of 13, were retrieved as free variables from an R_rs(λ) spectrum to achieve the highest retrieval accuracy from MuPI. This also controls the number of unknowns in the model to maximize its application to different satellite remote sensing data and improve computation efficiencies.

$$\mathrm{UAPD}=\frac{\left|\widehat{S}-S\right|}{0.5(\widehat{S}+S)}\times 100\%$$

(1)

where Ŝ is the estimated value and S is the measured value.

Once a_Gau(435) and a_Gau(617.6) were obtained from R_rs(λ), the other peaks could be estimated from these two values through a set of relationships shown in Table 2. Whether these relationships can represent data from different regions over different seasons presents uncertainty and a potential error source for a_Gau(λ) retrieval from MuPI. In our analysis of three datasets, we found that the empirical relationships among a_Gau(λ) to be very stable. The data from different regions over different seasons follow the trend predicted by the regression line based on one dataset. Figure 5 shows two examples in which all the data points from different datasets scattered tightly around the regression lines.

In contrast with the parameters determined in Wang et al. [9], to reduce potential errors updates were made in this study regarding the independent peak heights and their empirical relationships with other peaks. These modifications were made for the following reasons: (1) when further checking the co-variances among a_Gau(λ) from the MS and LT datasets, it was found that, sometimes, power-law relationships among the peak heights were more robust and representative than the previously adopted linear relationships (Figure 3); and (2) to make the retrievals more meaningful, a_Gau(435) and a_Gau(617.6) instead of a_Gau(515.6) and a_Gau(584.4) were used as the two independent variables. This change was made because a_Gau(435) showed a robust relationship with chlorophyll a concentration, and a_Gau(617.6) with phycocyanin. Further, their absorption coefficients or concentrations directly retrieved from R_rs(λ) can be more widely used in the estimation of phytoplankton biomass, productivity and cyanobacteria bloom detection and monitoring. Further, because of the high co-variation among the peak heights (R² > 0.92), the selection of different peaks as free variables caused no more than 10% difference in the retrieved results.

2.4. a_Gau(λ) Spectra Shape

In exploring the a_Gau(λ) application in HABs detection, a spectral shape analysis was conducted, in which a_Gau(435), a_Gau(584.4), and a_Gau(617.6) were normalized by their respective root of sum of squares,

$$n{a}_{\mathrm{Gau}}\left({\lambda }_i\right)=\frac{a_{\mathrm{Gau}}\left({\lambda }_i\right)}{\sqrt{{\displaystyle \sum _{i=1}^N{a}_{\mathrm{Gau}}{\left({\lambda }_i\right)}^2}}}$$

(2)

where the index N represents the total number of wavelengths, and λ_i corresponds to the wavelengths of 435, 584.4 and 617.6 nm. The na_Gau(λ) spectra vary over the range between 0 and 1, while it retains the “shapes” pertaining to the original a_Gau(λ) spectra, i.e., the band ratios of na_Gau(λ) retain the same as a_Gau(λ).

2.5. Multi-Spectral R_rs(λ)

To characterize the model performance with multi-spectral R_rs(λ) data as inputs, a set of different spectral band configurations were used (

Table 4. The multi-spectral bands used in this study: a composed All-band, OLCI, MERIS, MODIS, VIIRS, MSI and OLI bands.

Index	All-Band	OLCI	MERIS	MODIS	VIIRS	MSI	OLI
1	400	400
2	412	413	413	412	410
3	443	443	443	443	443	444	443
4	490	490	490	488	486	497	483
5	510	510	510
6	530			531
7	555	560	560	547	551	560	563
8	620	620	620
9	645			645
10	655						655
11	665	665	665	667		665
12	675	674		678	671
13	681	681	681
14	709	709	709			704
15	745	754	754	748	745	740

). This included the existing sensors of Ocean and Land Color Instrument (OLCI), MERIS, MODIS, Visible Infrared Imaging Radiometer Suite (VIIRS), Sentinal-2 Multispectral Instrument (MSI), Landsat-8 Operational Land Imager (OLI) and an All-band configuration that combined all of the existing spectral bands of the different sensors. The inclusion of the All-band configuration was intended to test the sensitivity of the algorithm to specific spectral bands of the existing sensors. The R_rs(λ) of these multiple spectral bands were obtained by applying their spectral response functions (RSF) (https://oceancolor.gsfc.nasa.gov/docs/rsr/rsr_tables/, https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library/-/asset_publisher/Wk0TKajiISaR/content/sentinel-2a-spectral-responses) to the in situ hyperspectral data.

3. Results

3.1. a_Gau(λ) from Multi-Spectral R_rs(λ)

The MuPI performance was validated for multi-spectral R_rs(λ) at six different band configurations.

Table 5. The mean (Mea.), median (Med.), maximum (Max.), and minimum (Min.) of unbiased absolute percentage differences (UAPD, %) of a_Gau(λ) estimated from R_rs(λ) with the spectral bands of different sensors for the MS, LT and LE datasets. MERIS+745: MERIS with the 745 nm band; MODIS-645/748: MODIS bands without 645 or 748 nm band.

Bands	UAPD (%)
	Mea.	Med.	Max.	Min.
ALL-band	30	32	36	23
OLCI	35	36	41	28
MERIS + 754	32	33	39	27
MERIS	35	36	41	28
MODIS	34	33	40	30
MODIS-645	34	33	43	29
MODIS-748	45	45	50	40
VIIRS	36	34	48	31
MSI	35	34	45	32
OLI	48	46	63	42

shows the mean UAPD of 13 a_Gau(λ) for all three in situ datasets with the spectral bands of different sensors, respectively.

Table 6. For LE 2014: The unbiased absolute percentage differences (UAPD, %) of a_Gau(λ) estimated from R_rs(λ) with the original hyperspectral bands and those at different spectral configuration. P1–P13: 13 Gaussian peaks. MERIS+745: MERIS with the 745 nm band; MODIS-645/748: MODIS bands without 645 or 748 nm band.

UAPD (%)	P1	P2	P3	P4	P5	P6	P7	P8	P9	P10	P11	P12	P13
ALL-band	29	23	23	19	20	20	14	22	19	18	36	21	18
OLCI	29	23	23	19	19	21	17	21	21	22	36	24	22
MERIS + 754	26	21	20	18	20	22	16	23	23	23	37	26	21
MERIS	27	21	22	17	19	21	17	23	23	24	37	25	22
MODIS	33	27	28	23	21	22	17	22	21	18	35	24	29
MODIS-645	32	26	24	24	19	22	17	18	17	17	31	24	27
MODIS-748	45	39	40	35	31	30	27	29	25	28	32	29	40
VIIRS	33	29	28	25	25	25	24	29	26	24	40	27	28
MSI	30	26	27	21	27	27	25	31	29	28	46	33	25
OLI	31	30	32	27	33	30	26	33	35	41	46	43	29

provides detailed information for LE 2014 data regarding the variation of the mean UAPD in the a_Gau(λ) retrieval with R_rs(λ) at different spectral configurations. OLCI and MERIS produced similar mean UAPD, which is around 35% for all three datasets and 28% for the LE 2014 data. Compared with MODIS, VIIRS showed slightly higher mean UAPD (All data: VIIRS: 36%, MODIS: 34%; LE 2014: VIIRS: ≤40%, MODIS: ≤37%), which is potentially caused by the lack of a spectral band around 665 nm in VIIRS. MSI and OLI produced the highest uncertainties for a_Gau(λ) at wavelength longer than 550 nm, but, surprisingly, a ≤50% mean UAPD was achieved for the four spectral bands of the Landsat-8 OLI sensor with the data used here, which demonstrates the potential of applying MuPI to these high spatial resolution satellite remote sensing data.

By removing and adding specific bands such as 645 and 748 nm for MODIS, and 754 nm in MERIS, differences in the a_Gau(λ) accuracy were observed. A 10% increase in mean UAPD was observed when the 748 nm band was removed from MODIS, and a 3% decrease in the mean UAPD when adding the 745 nm band to MERIS. However, with the existence of the 748 nm band, the adding and removing of the 645 nm band did not show a large influence (<3% mean UAPD variation). In addition, due to the existence of the 746 nm band in VIIRS, reasonable results were obtained in a_Gau(λ) retrievals.

Such a result suggests, at least for this dataset of bloom waters, with MuPI there is little impact on the retrieval of a_Gau(λ) when hyperspectral R_rs(λ) is degraded to multiple-spectral measurements such as MERIS and MODIS-like sensors. This is consistent with conclusions from previous studies that R_rs(λ) data do not need to be as fine as 1 nm in spectral resolution to obtain reliable retrievals of inherent optical properties [41,42,43].

When taking a closer look at the results from the LE 2014 data (Table 6, Figure 6), a mean UAPD of ≤35% was achieved for most of a_Gau(λ) from R_rs(λ) with MERIS, OLCI, and MODIS bands. The relatively low contributions of phytoplankton absorption coefficients at longer wavelengths was the main reason for the relatively larger uncertainties (e.g., a_Gau(λ) at 653 nm: ~36% for MODIS and MERIS). In addition, the lack of enough bands between 570 and 660 nm (only 620 nm for MERIS and 645 nm for MODIS) made it difficult to obtain more accurate retrievals.

3.2. Chl-a and PC from a_Gau(λ) for Lake Erie

As one of the most widely used indices for phytoplankton, concentration and light absorption coefficients of chlorophyll a (Chl-a) have been the focus of many studies [33,44,45,46,47]. Previous studies [9,10,17] have indicated that some Gaussian peak heights (a_Gau(λ)) obtained from a_ph(λ), such as peaks around 390, 413, 435 and 675 nm, represent the absorption properties of chlorophyll a. Phycocyanin (PC), the bio-marker of the blue-green cyanobacteria, is an important indicator for cyanobacteria biomass [6,7,48]. Based on the regression analysis between a_Gau(λ), Chl-a and PC from the LE dataset, the Gaussian peaks at 386.6, 414, 435 and 677 nm showed high correlation coefficients with Chl-a concentration with R² of 0.95, 0.96, 0.97, and 0.98, respectively; and the Gaussian peaks at 617.6 nm showed a 0.93 correlation coefficient with PC concentration. The empirical relationships for Chl-a and PC estimation from a_Gau(677) and a_Gau(617.6) were obtained, as shown in Figure 7. The above results further highlight the values of using a_Gau(λ) as a proxy to obtain pigment absorption or concentration, which may be used to map cyanobacteria bloom waters from ocean color imagery.

3.3. a_Gau(λ) versus Group Composition

In analysis of the a_Gau(λ) and cell counts from LE 2013 measurements, we found the spectral shapes formed by a_Gau(435), a_Gau(584.4), and a_Gau(617.6) vary with the composition of cyanobacteria species at different locations in western basin of Lake Erie, as shown in Figure 8. As shown in Table 2, the a_Gau(435), a_Gau(584.4), and a_Gau(617.6) contain absorption properties of pigments Chl-a, PE, and PC, respectively. The variation of pigment compositions, especially the different intracellular Chl-a:PE:PC ratios of the cyanobacteria at the species level is the main reason of the spectral shape variation. This can potentially be used in separating different cyanobacteria species in the bloom waters.

3.4. Pigment and a_Gau(λ) Retrieval from Satellite Imagery

3.4.1. Validation of Satellite Remote Sensing Data

Before applying MuPI to satellite imagery for a_Gau(λ) retrieval, the R_rs(λ) data from HICO and MODIS imagery for Lake Erie were first assessed with in situ measurements.

Three criteria were employed to find matchups: (1) Within ±2 days; (2) median of a 3 × 3 box, with no masks of land or clouds; and (3) coefficients of variance smaller than 0.15. The spectral response function was applied to the in situ R_rs(λ) spectrum before comparing it with satellite data. Most of the bands have a mean difference within 35% for R_rs(λ) for both sensors. At the shorter wavelengths (≤500 nm for HICO, and 412 and 443 nm for MODIS), the mean UAPD was as high as 65% and some satellite R_rs(λ) were even negative values because of the poor atmospheric correction for optically complex inland waters [49,50]. For MODIS, the mean UAPD of R_rs(λ) at 748 nm was 52%, possibly due to low R_rs(λ) values at this band and the influence from residual and likely uncorrected oxygen and water vapor absorption in the atmosphere.

Since the main focus of this work is the retrieval of a_Gau(λ) from R_rs(λ), we took advantage of the a_Gau(λ) information from the R_rs(λ) spectrum in the longer wavelengths (>500 nm for HICO and >480 nm for MODIS) in the model application, to minimize the influence of ineffective atmospheric correction on the shorter bands of satellite R_rs(λ). To validate this adjustment, MuPI was first applied to in situ measured R_rs(λ) with HICO and MODIS bands but without data from the shorter wavelengths (≤500 nm for HICO and 412 and 443 nm for MODIS). The retrieved a_Gau(λ) from such R_rs(λ) agreed quite well with a_Gau(λ) from in situ a_ph(λ) decomposition, with the mean UAPD ≤38% for HICO, and ≤48% for MODIS bands (

Table 7. The mean (Mea.) and median (Med.) of unbiased absolute percentage differences of a_Gau(λ) estimated from R_rs(λ) of HICO and MODIS satellite imagery.

Peaks	aGau(λ)	HICO		MODIS
		Mea.	Med.	Mea.	Med.
1	386.6	27	26	28	27
2	414	21	20	23	19
3	435	19	18	20	16
4	451.7	17	15	18	15
5	484	18	16	17	13
6	515.6	19	14	19	14
7	548.8	24	23	33	32
8	584.4	29	26	44	45
9	617.6	34	28	45	37
10	636	37	39	48	53
11	653	32	25	34	29
12	677	24	25	26	27
13	693.5	54	48	78	75

). The estimated Chl-a concentration has a mean UAPD of 28% which is better than the results from the standard products [45] and those shown in Pan et al. [8], and the mean UAPD for the estimated PC concentration is 32% which is consistent or even better than that reported in the literature [6,7,48]. Figure 9 presents one match-up of satellite and in situ R_rs(λ) spectra and the estimated R_rs(λ) spectra from MuPI for both HICO and MODIS, as well as, the 13 derived a_Gau(λ) from the corresponding spectrum. The a_Gau(λ) from the corresponding satellite and in situ match-up data showed the same trend.

3.4.2. Chl-a and PC from HICO and MODIS Imagery

We further explored the a_Gau(λ) distribution obtained from HICO and MODIS of the western basin of Lake Erie. The MuPI scheme was applied without using R_rs(λ) data at the shorter spectral bands (those <500 nm for HICO and <480 nm for MODIS). The power-law relationships obtained in Section 3.4 between a_Gau(677) and Chl-a, and a_Gau(617.6) and PC concentration were applied to the obtained a_Gau(λ) images to map the spatial distribution of Chl-a and PC concentration (Figure 10).

The estimated spatial distributions of Chl-a and PC from MODIS showed a similar pattern with those from HICO. A two-day difference exists between the HICO (15 August 2014) and MODIS (13 August 2014) observations due to the availability of satellite image, which explains the slightly different locations of the high biomass patches as shown in Figure 9. The non-value patches with in HICO image was a result of the poor R_rs(λ) quality due to the failing of the atmospheric correction.

3.4.3. MODIS and MERIS Imagery over Lake Erie

MERIS imagery did not coincide with the in situ observations used in this study. Thus, to provide an example of MERIS imagery, a_Gau(λ) were retrieved from MERIS imagery on 3 September 2011. MODIS imagery of this day was also considered.

As shown in Figure 11, similar patterns for high and low absorption patches were noticed for the two sensors, but the two images showed different magnitudes of a_Gau(λ). In further analysis of the R_rs(λ) from these two sensors at the same locations, good agreement was observed, as shown in Figure 12. The existence of 709 nm band in MERIS could be the main reason for the differences in the retrieval results. The patch of high values in the MERIS (marked with solid star: ★) imagery was not included in MODIS because the standard L1B to L2 processing in SeaDAS masked those pixels as clouds. Further validation and evaluation of MERIS and MODIS R_rs(λ) with in situ data are necessary to have a better understanding of the differences in the results, which is beyond the scope of the current work.

3.4.4. The Seasonal Variation of a_Gau(λ) from MERIS Imagery

We further explored the seasonal variation of a_Gau(435) and a_Gau(617.6) in Lake Erie retrieved from MERIS seasonal composed 4 km R_rs(λ) imagery (Figure 13). The differences in the spatial distribution of high pigment absorption patches in Spring, Summer, Autumn and Winter were captured. An obvious summer bloom in western basin of Lake Erie is shown in the figure.

4. Discussion

4.1. Spectral Requirements for a_Gau(λ) Retrieval

In the retrieval of a_Gau(λ) from in situ R_rs(λ) data, we found that a band around 695–715 nm is important for accurate a_Gau(λ) estimation in bloom waters. This is consistent with previous studies for the estimation of Chl-a in turbid productive waters [51,52,53,54]. Fundamentally, for phytoplankton bloom waters, the reflectance at wavelengths 695–715 nm can be augmented in the same way as occurs with terrestrial plants [51,52,53]. Lacking a proper spectral band within this spectral region for MODIS, the near-infrared band at 748 nm was included when inverting MODIS R_rs(λ). Although higher uncertainties at 748 nm band were noticed (Section 3.4.1), the inclusion of the band at 748 nm results in much more reasonable a_Gau(λ) retrieval from MODIS-Aqua measured R_rs(λ) (Figure 9, Table 7) versus without R_rs(λ) at this band. The same results were noticed in Section 3.1 for the a_Gau(λ) retrieval from R_rs(λ) at VIIRS and MSI spectral bands.

4.2. The seasonal a_Gau(λ) Variation in Lake Erie

The MERIS derived seasonal a_Gau(435) and a_Gau(617.6) variation in Lake Erie follows the pattern recorded in the literature [55,56]. In the central basin of Lake Erie, the different patterns, as shown in Figure 13, are due to the greater nutrients in the spring and autumn, and lower availability in the summer as a result of water stratification. The nutrient inputs due to agricultural activities as well as an expanding non-native mussel population, along with the light and temperature changes in different seasons form the main drivers for the different algal bloom patterns in Lake Erie during the four seasons [55,56,57,58,59]. As discovered in Moon [55], the biomass and taxonomic composition and the dominant taxa of surface assemblages varied in different seasons, which can be explained by the light, temperature and nutrient combinations in different seasons, and the strong spatial variability associated with mesoscale physical processes such as upwelling and basin-scale circulation [59].

4.3. Pigment Retrieval and HABs Detection

In this study, MuPI as a semi-analytical inversion scheme was applied to retrieve multiple Gaussian curves from satellite remote sensing data. As demonstrated in previous study [10], these Gaussian curves are related to the different phytoplankton pigments.

This phytoplankton pigment information has been used in the quantification of phytoplankton community composition, at least to a functional group level ([5] and the references therein), because many pigments are particular to specific taxonomic groups or even species [1]. However, only limited work has been conducted to obtain these phytoplankton pigments from satellite remote sensing data. Pan et al. [8] and Moisan et al. [60], as two of them, attempted to obtain 12 and 18 different phytoplankton pigments from satellite remote sensing data, respectively, but both of them are empirical approach based. Pan et al. [8] proposed using empirical relationships of 12 different phytoplankton pigments with the same band ratios of satellite R_rs(λ) around 490 nm to 550 nm. However, using 490 and 550 nm alone is not a good strategy for multiple pigments, as different pigments have different absorption peaks and troughs at different wavelengths [1]. Compared with Pan et al. [8], Moisan et al. [60] directly used the Chl-a product of satellite remote sensing to estimate a_ph(λ), then decompose this a_ph(λ) to 18 pigments based on their specific absorption coefficients. However, this satellite Chl-a product used similar empirical algorithm and spectral bands as in Pan et al. [8]. As both works were focused on coastal waters, another large uncertainty comes from non-algal particles and gelbstoff in these waters, which have a big influence on wavelengths shorter than 550 nm, and this influence cannot be eliminated by the band-ratio based algorithm. Compared with these works, MuPI not only considers the contribution of non-algal particles and gelbstoff in the coastal or inland waters, but also the different absorption properties of pigments, as shown in Table 1 and Table 2.

Another application of the pigment information is in HAB detection. The algorithms for HAB detection are usually based on Chl-a and its anomalies [47], or marker pigments, such as phycocyanin (PC) for cyanobacteria [6,7,18,48,61]. With MuPI, a successful discrimination of phycocyanin (PC) and Chl-a was achieved in this study. Limited by the dataset, the estimation of other pigment concentrations (phycoerythrin, chlorophyll b and c, carotenoids) from a_Gau(λ), as demonstrated in Hoepffner and Sathyendranath [10], could not be conducted here.

Cyanobacteria dominated HABs are increasing globally and presenting a major environmental and human health issue. Extensive Microsystis blooms with toxin production occur during summer and fall in different regions around the world and Microcystis contamination has been documented at many regions including Pinto Lake (California), Lake Erie (U.S.A.) and Lake Taihu (China) [6,7,21,61]. Other common bloom-forming pelagic genera include Aphanizomenon, Anabaena, Rhodomonas and Planktothrix. However, since toxicity is primarily associated with Microcystis, these other cyanobacteria blooms are generally considered nuisance blooms which will not cause acutely dangers to humans and wildlife [62], but they are frequently present in impacted water bodies (Figure 8). However, the ability to discriminate the different bloom-causing species is one of the challenges that existing algorithms are facing. In our analysis with a_Gau(λ) in Section 3.3, the possibility of discriminating different cyanobacteria species was shown. The variation of the spectral shape defined by a_Gau(435), a_Gau(584.4), and a_Gau(617.6) was found vary with different species and their composition in the water body as a result of the variation in pigment ratios [6]. This result showed the potential of MuPI in the application of separating species in HAB waters, which will be useful in detecting and monitoring potential toxin producers. However, because of data limitation, this potential is not fully addressed in this study, but it will be further explored with a larger dataset in the future.

5. Conclusions

In this study, the MuPI model was validated for obtaining the peak heights of Gaussian curves (a_Gau(λ)) from multi-spectral satellite remote sensing data. The model performance was validated in the retrieval accuracy of a_Gau(λ) with R_rs(λ) of six multi-spectral band configurations, and the spectral requirements were discussed. Less than 35% of mean unbiased absolute percentage differences were achieved for a_Gau(λ) from R_rs(λ) spectra with OLCI, MERIS and MODIS bands, and less than 45% for VIIRS, MSI and OLI bands. Using data from the western basin of Lake Erie as an example, the R_rs(λ) obtained from HICO and MODIS satellites were validated with in situ data over cyanobacterial bloom waters in Lake Erie, and the spatial distributions of a_Gau(λ) and the concentrations of chlorophyll a and phycocyanin were obtained, where the patches of cyanobacteria bloom were clearly presented. A seasonal distribution of pigment absorption coefficients was obtained from MERIS seasonal composed imagery of 2011 for Lake Erie. These results demonstrate that, with MuPI, it is possible to analytically retrieve information of not only Chl-a, but also PC and potentially other pigments, which will significantly enhance our capability to characterize and evaluate the status of phytoplankton blooms and discriminate phytoplankton groups using satellite ocean color remote sensing.