Comparative Analysis of Empirical and Machine Learning Models for Chla Extraction Using Sentinel-2 and Landsat OLI Data: Opportunities, Limitations, and Challenges

Abstract

Remote retrieval of near-surface chlorophyll-a (Chla) concentration in small inland waters is challenging due to substantial optical interferences of various water constituents and uncertainties in the atmospheric correction (AC) process. Although various algorithms have been developed to estimate Chla from moderate-resolution terrestrial missions (∼10–60 m), the production of both accurate distribution maps and time series of Chla has proven challenging, limiting the use of remote analyses for lake monitoring. Here, we develop a support vector regression (SVR) model, which uses satellite-derived remote-sensing reflectance spectra ($R_{rs}^{\delta }$) from Sentinel-2 and Landsat-8 images as input for Chla retrieval in a representative eutrophic prairie lake, Buffalo Pound Lake (BPL), Saskatchewan, Canada. Validated against in situ Chla from seven ice-free seasons (N ∼ 200; 2014–2020), the SVR model outperformed both locally tuned, $R_{rs}^{\delta }$-fed empirical models (Normalized Difference Chlorophyll Index, 2- and 3-band, and OC3) and Mixture Density Networks (MDNs) by 15–65%, while exhibiting comparable performance to a locally trained MDN, with an error of ∼35%. Comparison of Chla retrieval models, AC processors (iCOR, ACOLITE), and radiometric products (Rayleigh-corrected, surface, and top-of-atmosphere reflectance) showed that the best Chla maps and optimal time series (up to 100 mg m⁻³) were produced using a coupled SVR-iCOR system.

RÉSUMÉ

L’extraction à distance de la concentration de chlorophylle-a (Chla) près de la surface dans les petites eaux intérieures est difficile en raison des interférences optiques importantes de divers constituants de l’eau et des incertitudes dans le processus de correction atmosphérique (CA). Bien que divers algorithmes aient été développés pour estimer Chla à partir de missions terrestres à résolution modérée (∼10–60 m), la production de cartes de répartition précises et de séries chronologiques de Chla s’est avérée difficile, limitant l’utilisation d’analyses à distance pour la surveillance des lacs. Ici, nous développons un modèle de régression vectorielle de support (RVS), qui utilise des spectres de réflectance dérivés de satellites ($R_{rs}^{\delta }$) utilisant des images Sentinel-2 et Landsat-8 comme entrée pour la récupération de la Chla d’un lac eutrophique des prairies représentatif, Buffalo Pound Lake (BPL), Saskatchewan, Canada. Validé à partir des données Chla in situ de sept saisons sans glace (N ∼ 200; 2014–2020), le modèle SVR a surpassé à la fois les modèles empiriques à réglage local, alimentés en ($R_{rs}^{\delta }$) (indice de chlorophylle par différence normalisée, bandes 2 et 3 et OC3) et les réseaux de densité de mélange (MDN) de 15% à 65%, tout en présentant des performances comparables à celles d’un MDN formé localement, avec une erreur de ∼35%. La comparaison des modèles de récupération de Chla, des processeurs AC (iCOR, ACOLITE) et des produits radiométriques (correction de Rayleigh, réflectance de surface et de la haute atmosphère) a montré que les meilleures cartes Chla et les séries chronologiques optimales (jusqu’à 100 mg m⁻³) ont été produites à l’aide d’un système SVR-iCOR couplé.

Introduction

Small inland waters (SIWs) are the predominant form of lakes globally, with 64% of basins <100 km² (Downing et al. 2006), yet they are highly subject to water quality degradation due to changes in climate and land use (Carpenter et al. 1998, Delpla et al. 2009). Despite recognition of the problem for decades, the water quality of SIWs continues to degrade, resulting in harmful algal blooms (HABs) composed of cyanobacteria (Walker 2019). The frequency, magnitude, and persistence of HABs have also increased globally due to atmospheric warming (Ho et al. 2019; Hayes et al. 2020). A change in the near-surface concentration of chlorophyll-a (Chla) is one of the most reliable proxies of algal bloom intensification retrievable from satellite analyses, as Chla is present in all phytoplankton, including cyanobacteria (Roesler et al. 2017), and has unique absorption features (peak at ∼430 and ∼670 nm in live organisms) that can be detected through optical imaging (Kutser 2009).

Accurate Chla retrieval from optical radiometry is affected by the interplay between the inherent optical properties (absorption, scattering) of pure water, its dissolved or suspended constituents, and solar photons in water-leaving radiance ($L_w$). In particular, reflectance is affected strongly by phytoplankton density, colored dissolved organic matter (CDOM), and non-algal particles (NAP) (Babin et al. 2003). Further, $L_w$ is modulated by atmospheric properties of the transmission path to satellite sensors, consequently, atmospheric correction (AC) processors are used to convert top-of-atmosphere reflectance (${\rho }_{\mathit{TOA}}$) to satellite-derived remote sensing reflectance ($R_{rs}^{\delta }$) to retrieve Chla. Estimates of ${ R}_{rs}^{\delta }$ include uncertainties in the AC and sensor radiometric measurements, as well as the effect of surface-reflected radiance (sun glint; Bulgarelli and Zibordi 2018), but are useful estimates of remote sensing reflectance ($R_{rs}),$ defined as the ratio of water-leaving radiance to the total downwelling irradiance just above water. Once $R_{rs}^{\delta }$ is approximated, a wide range of algorithms, including semi-analytical, empirical, and machine-learning (ML) models, can be applied to retrieve Chla from reflectance measurements (Carder et al. 1999; Morel 1980; Odermatt et al. 2012).

Semi-analytical models retrieve water absorption and scattering properties from $R_{rs}$ measurements and can be used to estimate Chla (Gons 1999; Lee et al. 2002; Schroeder et al. 2007). While accurate in some circumstances (Santini et al. 2010; Van Der Woerd and Pasterkamp 2008), these models are sensitive to the form of AC and require accurate estimates of optical water parameters (IOCCG 2006; Odermatt et al. 2012). In contrast, empirical models (differential/ratio-based indices) based on blue-green wavelengths (e.g., NASA’s OCx models) tend to perform well in phytoplankton-dominated aquatic ecosystems (O'Reilly et al. 1998; O'Reilly and Werdell 2019). Further, various red- and near infrared (NIR)-indices have been developed and validated for ocean color sensors, including the 2band, 3band, and Normalized Difference Chlorophyll Index (NDCI) (Dall’Olmo and Gitelson 2005; Mishra and Mishra 2012; Moses et al. 2009) for use with data from the Medium Resolution Imaging Spectrometer (MERIS). Models based on red or NIR bands may be less sensitive to uncertainties in AC, especially when closely spaced (Moses et al. 2009); nonetheless, model performance depends on the range of Chla variation, the amount of interference from other constituents (e.g., backscattering NAP), and the band configuration of sensors (Gitelson 1992; Gitelson et al. 2007). Instead, empirical ML algorithms, especially neural networks (NN), are widely used to retrieve Chla over geographically-extensive regions using large synthetic or in situ radiometric measurements from diverse optical water types (OWTs) (Doerffer and Schiller 2007; Hu et al. 2021; Pahlevan et al. 2020; Schroeder et al. 2007).

To date, remotely-sensed Chla estimates have been applied successfully to large waterbodies, including the open ocean (Bryan et al. 2005; O'Reilly and Werdell 2019), coastal waters (Moses et al. 2012; Werdell et al. 2009), and large lakes (Binding et al. 2021; Binding et al. 2011; Gons et al. 2008; Schaeffer et al. 2018), using ocean-color sensors, such as MERIS and the Sea-viewing Wide Field-of-view Sensor (SeaWiFS). In contrast, Chla retrieval for SIWs has been challenging because the optical regimes of inland waters are influenced by particulate organic and inorganic particles, as well as CDOM (Mobley 1994). Generally, ocean-color sensors lack sufficiently high spatial resolution (<100 m) to sample SIWs (Ansper and Alikas 2018; Philipson et al. 2014). Likewise, very high-resolution sensors are not widely used to monitor water quality mostly because they do not offer much improvement in terms of spectral resolution and signal to noise ratio (SNR), despite their commercial nature. Instead, Multi-Spectral Instrument (MSI) and Operational Land Imager (OLI) sensors onboard Sentinel-2 (S2) and Landsat-8 (L8) satellites show potential for sensing Chla from SIW, as they provide excellent global coverage with spatial resolutions from 10 to 60 m. Although designed for land observations, these sensors may be applicable to small aquatic ecosystems (Cao et al. 2019; Pahlevan et al. 2014; Xu et al. 2020), primarily because of their radiometric performance and stability (Claverie et al. 2018; Helder et al. 2018; Pahlevan et al. 2019; Wulder et al. 2015) compared to heritage Landsat-class missions (Allan et al. 2011; Tebbs et al. 2013; Yacobi et al. 1995). Once combined, S2 and L8 images are available at sub-weekly revisit rate in high-latitude regions (Li and Roy 2017), favoring their use for water quality and HAB monitoring.

A broad range of Chla models has been used for MSI and OLI images, including 2band, 3band, and NDCI (Ansper and Alikas 2018). While MSI has been utilized for detecting cyanobacterial blooms and retrieval of Chla in subalpine lakes (Bresciani et al. 2018), studies suggest that current approaches have limitations at the extremes of the observed Chla range, e.g., Chla < 10 mg m⁻³ or Chla > 100 mg m⁻³ (Toming et al. 2016; Dörnhöfer et al. 2016; Kutser et al. 2016). Instead, the application of Mixture Density Networks (MDN) to a large dataset of in situ radiometry and Chla measurements has allowed the development of models which outperformed other state-of-the-art algorithms for a wide range of Chla concentration (0.1–100 mg m⁻³) using MSI (Pahlevan et al. 2020), as well as OLI data (Smith et al. 2021). Additionally, Cao et al. (2020) developed BST, a model based on the Gradient Boosting Tree algorithm (XGBoost) (Chen and Guestrin 2016), and successfully tested it on OLI data taken from lakes in eastern China. As another example of ML models employed to retrieve Chla, Support Vector Machines/Regressions (SVM/SVR) (Vapnik 2013) have been applied to oceanic (Camps-Valls et al. 2006; Hu et al. 2021; Kwiatkowska and Fargion 2003; Martinez et al. 2020) and inland waters (Tian et al. 2022).

Despite recent developments, reliable estimates of Chla in SIWs remain challenging when based on moderate-resolution (∼10–60 m) satellite data. Empirical models leverage only a limited range of the spectrum and may not optimally solve ill-posed conditions (O'Sullivan 1986) that are common to inverse problems, such as Chla retrieval (Defoin‐Platel and Chami 2007; Pahlevan et al. 2020; Sydor et al. 2004; Werdell et al. 2018). Similarly, while globally trained ML (GML) models (e.g., MDN) can leverage the full visible and near-infrared spectrum (VNIR) and may handle non-linear and ill-posed problems (Pahlevan et al. 2020), they can be susceptible to uncertainties in AC that could reduce their suitability under sub-optimal atmospheric or aquatic conditions (Pahlevan et al. 2020; Smith et al. 2021). These observations suggest that the development of locally trained ML (LML) models using $R_{rs}^{\delta }$ measurements might be suitable solution for optimal monitoring of Chla at local scales.

Here, we employed an ML approach based on SVR to retrieve robust and reliable Chla time series and maps for Buffalo Pound Lake, Saskatchewan, Canada, using MSI and OLI imagery. Our SVR model was trained and validated with ∼200 co-located in situ Chla measurements with corresponding $R_{rs}^{\delta }$ observations. We compared model performance against several state-of-the-art algorithms, including OC3, MDN, 2band, BST, and LMDN—a locally trained MDN—in terms of its quantitative (general and stratified) performance, as well as its spatial and temporal consistency. Then, we assessed the robustness of the model for uncertainties from two AC processors (i.e., iCOR and ACOLITE) and different broadly-defined OWTs, to assess its potential utility for other small eutrophic lakes. Our overall objective was to develop a reliable baseline model for BPL that might also be suitable for other regional lakes exhibiting similar HABs and water conditions.

Materials and methods

Study site

Buffalo Pound Lake (BPL) is a long (∼30 km), narrow (<1 km), and shallow (<6 m) lake located in the Qu’Appelle River watershed, Saskatchewan, Canada (Figure 1, Table 1). Currently, the basin is eutrophic, with summer blooms occurring during June-September and peak surface populations of phytoplankton during July–August (Kehoe et al. 2019). Continuous monitoring for over 25 years shows that cyanobacteria are the predominant phytoplankton during July-September (Swarbrick et al. 2019; Vogt et al. 2018). The lake landscape orientation parallel to the direction of prevailing winds means that the water column is polymictic, experiencing frequent mixing periods with only intermittent vertical stratification (Dröscher et al. 2008).

Figure 1. Map and location of Buffalo Pound Lake (BPL), Saskatchewan, Canada. (a) Location of the Qu’Appelle River watershed within Canada. (b) Location of BPL within Qu’Appelle River watershed. (c) A Landsat-8 RGB image of BPL is overlaid on a bathymetric map on which sampling stations are also shown (solid black triangles numbered 1–11).

Table 1. Buffalo Pound Lake characteristics and water quality parameters at station 1 (averaged from late May to early September  2014–2020).

Physical characteristic	Value	Water quality parameter	Mean value
Altitude (m)	501	Water temperature (°C)	19
Mean depth (m)	3.8	Total P (mg m^−3)	88
Maximum depth (m)	5.8	Chla (mg m^−3)	31
Length (km)	30	TSS (g m^−3)	11.8
Average width (km)	0.9	$a_{\mathit{CDOM}}$ (440) (m^−1)	0.28
Surface (km^2)	30	SDD (m)	0.98
Volume (m^3)	92 × 10^6	Turbidity (NTU)	9.2

P stands for phosphorus, TSS for total suspended solids, and SDD for Secchi Disk Depth. $a_{\mathit{CDOM}}$ (440) is CDOM absorption measured at 440 nm.

Several attributes make BPL suitable for the development of remote sensing models of Chla. First, the lake is an important freshwater resource as it supplies drinking water to one-quarter of the provincial population, including the nearby cities of Regina and Moose Jaw (Hosseini et al. 2018). Second, multi-decadal records demonstrate elevated Chla content during late summer, with abundant surface blooms of toxic cyanobacteria (Kehoe et al. 2015; Hayes et al. 2020). Third, the lake size and elongated shape produce large gradients and patches of differing Chla concentration (10–100 mg Chla m⁻³), as well as regions of contrasting optical properties (NAP turbidity, HABs) that are suitable for analysis with spatially resolved MSI and OLI platforms. Finally, BPL is representative of many other prairie lakes in terms of physical, biological, and chemical properties (Finlay et al. 2015; Hayes et al. 2020), suggesting that models developed in this site may have regional suitability for water quality monitoring.

BPL exhibits two distinct OWTs (Appendix A). OWT1 characterizes the southern basin (stations 1–8), where Chla are elevated and optical characteristics are similar to those recorded in phytoplankton-rich systems elsewhere (OWT4 in Pahlevan et al. 2021; OWT8 in Spyrakos et al. 2018). In contrast, the northern basin (stations 9–11) exhibits of suspended sediments and lower Chla values (Table A1 in Appendix A), similar to OWT5 in Pahlevan et al. (2021) or OWT4 in Spyrakos et al. (2018).

Data

Although there is a long history of recorded in situ data in BPL (Swarbrick et al. 2019), we selected the period of 2014–2020 to match Landsat-8 and Sentinel-2 missions.

In situ Chla data

In situ Chla data originated from multiple datasets (Table 2). At station 1, autonomous, on-site fluorescence probes were available through deployment on a buoy. These fluorometric measurements were calibrated following Chegoonian et al. (2022). In addition, discrete water samples were collected from the lake surface and 0.8-m depth, with Chla collected on Whatman GF/F frozen and later extracted following Wintermans and De Mots (1965) and analyzed using a UV-visible spectrophotometer (Shimadzu UV-1601-PC). Samples from station 2 were obtained from the water treatment plant intake at ∼3 m depth in this polymictic lake. Samples from the intake were filtered onto a 0.45-µm pore filter, extracted in 90% acetone, and analyzed via spectrophotometry following standard methods (Eaton et al. 2017).

Table 2. Details of in situ Chla measurements employed in this study.

Station	Measurement method	Available years	Sampling interval	Sampling time	Depth
1	Field fluorometry	2014–2020	10 min	All-day	0.8
1	Lab spectrophotometry	2017–2020	∼Once/week	Daytime	0, 0.8
2	Lab spectrophotometry	2014–2020	Once/week	∼5 a.m.	2.8
3	Lab HPLC	2014–2019	Twice/month	∼11 a.m.	0–1*
4–11	Lab spectrophotometry	2015–2020	Once/month	9 a.m.–1 p.m.	1

The unit for depth values is meter.

*Depth-integrated sampling.

Phytoplankton from station 3 was collected on GF/C glass-fiber filters (nominal pore size 1.2 µm) following Swarbrick et al. (2019). Briefly, surface water (∼0.5-m depth) and depth-integrated samples were filtered through GF/C filters and frozen (−10 °C) until analysis for Chla (mg m⁻³) through standard trichromatic assays (Jeffrey and Humphrey 1975) and biomarker pigment (nmoles pigment L⁻¹) analysis by HPLC following Leavitt and Hodgson (2001).

Samples from stations 4 to 11 were collected during monthly field visits at a 1-m depth using a Niskin bottle. Sub-samples for Chla analysis were transferred into laboratory bottles, stored in dark cool, containers, and analyzed using Eaton et al. (2017) method 10200H. Briefly, samples were filtered at low vacuum through 0.45 µm nitrocellulose filters, and pigments were extracted using a 90% acetone solution by mixing. The resulting samples were steeped for <24 hours before Chla values were calculated following Jeffrey and Humphrey (1975).

Satellite images

Cloud-free Level-1C MSI images acquired by the Sentinel-2A/B satellites with a 2–3 days revisit time during the open water season were identified manually and downloaded for the period 2017–2020. The MSI sensor collects data in 13 spectral bands from 443 to 2190 nm at spatial resolutions of 10, 20, and 60 m, and with a 12-bit radiometric resolution (Li et al. 2017). In addition, cloud-free OLI Level-1 images from Landsat-8 satellite (launched 2013) were downloaded for the period 2014–2020. The spatial resolution of the optical channels of OLI is 30 m, and the satellite overpasses the study site every ∼8 days. Appendix B compares the MSI and OLI’s spectral configuration with Chla spectral reflectance (Figure B1), including reflectance spectra for samples with different Chla measured in BPL using an ASD spectrometer (Analytical Spectral Devices, ASD Inc., Boulder, CO, USA).

Methodology

A similar data analysis workflow (Figure B2) was used for all analyses in this study, although algorithms (e.g., AC processors, Chla retrieval models) and train-test split approaches differed between experiments (see Table D1).

Data preprocessing

All images were corrected for atmospheric effects to produce two different reflectance quantities, namely satellite-derived remote sensing reflectance ($R_{rs}^{\delta }$) and Rayleigh-corrected reflectance (${\mathrm{\rho }}_{rc}$). We selected ACOLITE (v20210114.0) (Vanhellemont 2019; Vanhellemont and Ruddick 2014) and iCOR (version 3) (De Keukelaere et al. 2018) as AC processors since they outperform other processors in inland waters with OWTs similar to BPL (Pahlevan et al. 2021), especially when red-NIR wavelengths are used (Ilori et al. 2019). Visual inspection of images showed no significant sunglint effect in BPL; besides, a sunglint correction in the presence of adjacency effect (AE) may result in overcorrection (Vanhellemont 2019). The use of iCOR applies the SIMilarity Environment Correction (SIMEC) algorithm (Sterckx et al. 2015) to reduce AE which may be an issue for BPL due to its narrow width. Although ACOLITE lacks an inherent AE correction in the current version, a low threshold for the SWIR band (top of atmosphere reflectance at 1609 nm = 0.0215) was set to remove pixels highly impacted by AE and sunglint, as well as land pixels (Vanhellemont 2019). Furthermore, thanks to the dynamic band selection, the dark spectrum fitting (DSF) algorithm used in ACOLITE selects other bands (typically blue or red which might be unaffected by AE from nearby dark vegetation) if the NIR/SWIR adjacency effects are severe (Vanhellemont 2019). Regardless of AC processors, all MSI spectral bands were then resampled to a 60-m grid to be consistent for further steps (Ansper and Alikas 2018). Analysis of model performance using images resampled at 10, 20, and 60 m resolution (Figure B3) demonstrates that resampling at 60 m did not affect retrieval accuracy, yet improves time efficiency in model development.

Optically deep waters are the focus of this study; hence, Chla samples for which Secchi Disk Depth (SDD) measurements equal to bottom depth were excluded. This was to ensure that bottom reflection is avoided in our assessments. In situ samples (1394 station-day samples) were then collated with the closest matching satellite-derived $R_{rs}$ products to create co-located $R_{rs}^{\delta }$−Chla matchups. The maximum time span between field sampling and image acquisition was 3 days (median = 0 day). While longer than the ±3 hours interval recommended for oceanic waters (Werdell and Bailey 2005), this value is much shorter than the interval needed (up to ±7 days) for reliable retrieval from optically-stable inland waters (Ansper and Alikas 2018; Dörnhöfer et al. 2018; Lunetta et al. 2015; Tang et al. 2003). However, to minimize potential mismatches in the sampling date, we used continuous Chla from the buoy to exclude matchups for which Chla at the time of satellite overpass differed from in situ values by >20%. Representative $R_{rs}^{\delta }$ spectra for matchups were chosen to be the median of 3 × 3-element windows centered around the matchup locations.

Both AC processors mask land and clouds automatically; however, we manually deleted matchups that were contaminated by thin clouds/haze and cloud shadow through a visual assessment of images. Both processors occasionally overcorrect for atmospheric effects, mostly due to aerosol contribution, resulting in negative reflectance, especially in the 443 and 490 nm bands. However, in this study, there were few instances of negative reflectance values (∼5%) and these were excluded after inspection. Finally, we implemented an outlier detection algorithm to remove samples whose $R_{rs}^{\delta }$ deviated from the mean values of $R_{rs}^{\delta }$ by more than ±3σ. Approximately 200 matchups (varies by sensor type and AC processor) were selected for algorithm development and evaluation (Table D1). The distribution of $R_{rs}^{\delta }$ derived from ACOLITE and iCOR ($R_{rs}^{\delta ,\mathit{ACL}}$ and $R_{rs}^{\delta ,\mathit{iCOR}},$ respectively) is shown in Figure B4.

Model development

Input and output Chla values were log₁₀-transformed in the SVR model (see Appendix C). We allowed some outliers using a C = 2.5 parameter (regularization term) to decrease the chance of overfitting. We also employed a Radial Basis Function kernel (RBF) with γ = 0.14 and 0.25 for MSI and OLI data, respectively, to handle non-linearity in the feature space. These hyperparameters (C, γ, and kernel type) were tuned using a grid-search cross-validation process that minimizes model errors (mean absolute error) on a validation set. Here, the validation set was one-fifth of the training data (see Table D1) that was repeatedly set apart for hyperparameter tuning. After identifying optimized hyperparameters, the validation set was merged with the whole training data and fed into the model for a final training process.

Using $R_{rs}^{\delta }$−Chla matchups, we calibrated several state-of-the-art empirical Chla retrieval algorithms for use with MSI and OLI spectral bands, namely OC3 (O'Reilly et al. 1998; O'Reilly and Werdell 2019), 2band (O'Reilly et al. 1998), 3band (Dall’Olmo and Gitelson 2005), and NDCI (Mishra and Mishra 2012) for MSI, and OC3, as well as FLH-blue (Beck et al. 2016) for OLI. Although OC3 was originally developed for clear oceanic waters, this model is commonly used as a benchmark for Chla retrieval in inland waters (e.g., Pahlevan et al. 2020). After a log₁₀ transformation, these differential/ratio-based indices implied a linear relationship with log₁₀-transformed Chla. The exceptions were 2band and 3band for which we added a power-of-two term to better fit the data. The tuned formula and coefficients for empirical models are presented in Appendix D (Table D2).

We also applied MDN and BST models as representatives of state-of-the-art ML models developed for MSI and OLI. MDN was implemented using the code available via https://github.com/STREAM-RS/STREAM-RS (Pahlevan et al. 2020; Smith et al. 2021). In addition, we implemented a locally trained MDN (LMDN) using local $R_{rs}^{\delta }$−Chla matchups. A similar process was conducted for the BST model (Cao et al. 2020) using the BST-OLI package (https://github.com/zgcao/bst_oli) and a locally trained XGBoost model, LBST. The reflectance spectra imported into these LML models (LMDN and LBST) were identical to our SVR model; i.e., $R_{rs}^{\delta }$ derived from the first seven and four spectral bands (400–800 nm) for MSI and OLI, respectively.

Model assessment

Chla retrievals were assessed from three different aspects; quantitative performance, spatial integrity, and temporal validity. We also examined the robustness of our proposed model under various scenarios, including changes in water type, AC processors and radiometric products, and remote sensing data types.

For the quantitative assessment, the MSI datasets are selected as the main data source. Matchups were split into training and test datasets for the following experiments to estimate general and stratified performance, spatial and temporal integrity, model sensitivity to sensors and AC processors, and model transferability; however, the method used to do so differed among the experiments to assure a complete assessment of our model. Table D1 summarizes the evaluation approaches (training-test splitting), as well as the number of training/test matchups, available for each experiment.

Assessment of general performance (section General performance) and model transferability (section Model transferability over water types) was based on a cross-validation approach in which the matchups were categorized either annually (Table D3) or geographically (southern/northern basins). In each run, $R_{rs}^{\delta }$−Chla matchups related to a single year (or basin) were put aside as test data before the model was trained with the remaining data and used to assess model performance.

To gain insight into the model performance in two eutrophic conditions (OWTs; stratified performance hereafter; section Stratified performance), model sensitivity to the two AC processors (section Model sensitivity to AC and radiometric products), and its robustness for each sensor (section Model sensitivity to sensor type), we used a 5-fold cross-validation approach to randomly select among $R_{rs}^{\delta }$−Chla matchups. This approach ensures sufficient, equal training/test data for each run.

Assessment of model capability in generating Chla maps (section Spatial integrity) using both MSI and OLI images was based on images from a single date (July 16, 2020) when we had both cloud-free images from both sensors (∼10 minutes apart) and the maximum number of coincident (within 2 hours) in situ Chla samples (nine total), spanning a broad range of Chla (∼10–100 mg m⁻³). The corresponding matchups were considered equivalent to unseen test data, and the models were trained with the remaining matchups (184 matchups for MSI and 169 for OLI) (Table D1). In addition, to assess the stability of Chla retrieval over time (section Temporal validity), MSI-derived $R_{rs}^{\delta }$−Chla matchups corresponding to the continuous measurements of the buoy in 2020 were considered as unseen test data, and the remaining matchups were used to train the models.

Accuracy metrics

Both linear and log₁₀-transformed metrics were examined to assess model accuracy. In general, metrics calculated in log₁₀-transformed space (i.e., RMSLE, SSPB, and MdSA) are believed to provide a better assessment due to the log-normal distribution of Chla (O'Reilly and Werdell 2019; Seegers et al. 2018). The performance metrics for accuracy assessment were estimated as follows: (1) $$\mathit{RMSE}={\left[{\sum }_{i=1}^N{(P_i-M_i)}^2/n\right]}^{1/2}\left(\mathrm{mg}\mathrm{ }{\mathrm{m}}^{-3}\right)$$(1) (2) $$\mathit{RMSLE}={\left[{\sum }_{i=1}^N{\left({{\log }_{10}(P}_i\right)-{{\log }_{10}(M}_i\left)\right)}^2/n\right]}^{1/2}$$(2) (3) $$\mathit{MAPE}=100\times \mathit{median}\left(\right[|P_i-M_i|/M_i\left]\right)$$(3) (4) $$\mathit{SSPB}=\mathrm{\hspace{0.17em}100}\text{ × sign}\left(\mathrm{z}\right)\left({10}^{\left|z\right|}-1\right), z=\mathit{\hspace{0.17em}median}\left(\log 10\left(P_i/M_i\right)\right)\left(\%\right)$$(4) (5) $$\mathit{MdSA}=100\text{ × }\left({10}^y-1\right),\mathrm{ }\mathrm{ }y=\mathit{median}\left|\log 10\left(P_i/M_i\right)\right|\left(\%\right)$$(5) where $P_i$ and $M_i$ stand for predicted and measured Chla, respectively. RMSE is the root mean squared error, RMSLE is the root mean squared log-error, MAPE is the median absolute percentage error, SSPB represents symmetric signed percentage bias, and MdSA is the median symmetric accuracy, computed in log-space (Morley et al. 2018).

SSPB and MdSA were expressed as percent (%), expected to be resistant to outliers, zero-centered, and easily interpretable (Pahlevan et al. 2020). While SSPB measures the bias of a model, MdSA is believed to be an indicator of its precision. Because SSPB and MdSA are relatively new indices, we also estimated RMSE, RMSLE, and MAPE to facilitate the comparison with earlier studies. Finally, models were evaluated using Slope and Model Win Rate (MWR) criteria, wherein Slope was used to compare the results with earlier studies, while MWR, expressed in %, was used to determine which model performed better in pair-wise comparison of the residuals (Seegers et al. 2018).

Results

Quantitative assessment of the model on MSI data

Quantitative assessments were conducted using both general and stratified performance. Here, general performance analysis employed all matchups, whereas stratified analysis was conducted separately on two OWTs and provides insights into the use of SVR models in eutrophic conditions.

General performance

The overall accuracy of models for retrieving Chla from MSI-derived $R_{rs}^{\delta ,\mathit{ACL}}$ values was computed over all stations, and the whole Chla range (∼1–125 mg m⁻³; Table 3, Figure 2). Results show that LML models (SVR, LMDN) significantly outperformed (>15% improvement in MdSA) all other empirical and GML models. In particular, SVR outperformed all empirical models as reported via MWR, representing >60% of retrievals. Compared to LMDN, SVR performed marginally better (∼3% improvement in MdSA) but returned equal estimates of bias (as SSPB). The slope for SVR (0.78) demonstrates reasonable performance through the whole range of Chla in BPL. Among other models, the performance of OC3 was poor, as expected because of its dependency on blue-green band ratios, while other empirical models for eutrophic waters (2band, 3band, and NDCI) performed better and similarly in BPL, with the 2band algorithm generally outperforming other empirical models.

Figure 2. Matchup analysis of Chla derived from different algorithms applied on MSI-derived $R_{rs}^{\delta ,\mathit{ACL}}$ data and near-coincident, co-located in situ Chla samples in BPL. Year of data acquisition indicated by colored solid circles.

Table 3. Evaluation metrics (general performance) for Chla retrieval models on MSI and in situ Chla matchups (N = 193).

	MdSA (%)	SSPB (%)	RMSLE	MAPE (%)	RMSE (mg m^−3)	Slope	R^2	MWR (%)
OC3	99.1	−6.2	0.43	55.71	29.53	0.11	0.17	73.6
2band	50.3	−12.9	0.31	36.39	18.95	0.7	0.55	60.1
3band	54.2	−10.2	0.36	38.33	355.5	1.06	0.45	61.1
NDCI	52.3	−14.3	0.32	37.63	27.13	0.66	0.52	61.1
MDN	55.9	28.2	0.35	48.32	36.81	1.23	0.57	62.2
LMDN	38.9	3	0.29	33.22	17.51	0.83	0.63	52.9
SVR	35.6	3.4	0.27	31.5	13.96	0.78	0.66	N/A

The Model Win Rate (MWR) is computed relative to SVR as the reference model. Highlighted cells indicate the best score for the corresponding metrics.

The MDN model, trained on global $R_{rs}$ data, exhibited comparable precision to empirical models (∼56% error), albeit with a high bias (SSPB = 28%) and a tendency to overestimate Chla (Slope = 1.23) reflecting its sensitivity to $R_{rs}^{\delta }.$ The LMDN showed good performance, implying MDN’s strong performance even with a relatively small training sample size (∼10% of matchups used by Pahlevan et al. 2020).

Visual inspection of scatter plots revealed that SVR and LMDN predictions were less biased than other models when based on annual sampling. Scatter plots of Chla retrievals also illustrated a reasonable overall performance of all models (except OC3) for Chla > 10 mg m⁻³ ($\overline{\mathit{MdSA}}$ = 39.9 ± 10.32%). However, retrieving Chla < 10 mg m⁻³ was less accurate ($\overline{\mathit{MdSA}}$ = 57.8 ± 6.95%), with most models overestimating Chla in this lower range ($\overline{\mathit{SSPB}}$ = 46.8 ± 12.1%). Nonetheless, SVR and LMDN models exhibited a substantially better performance ($\overline{\mathit{Slope}}$ = 0.47 ± 0.04) when compared to the other models ($\overline{\mathit{Slope}}$ = −0.19 ± 0.09). While all models failed to retrieve Chla < 2 mg m⁻³, the absence of data in this range (2 matchups) prevented a detailed evaluation of performance. Instead, we infer from Figure 2 that empirical models underestimated Chla > 30 mg m⁻³ ($\overline{\mathit{SSPB}}$ = −43.4 ± 7.1%) especially values >100 mg m⁻³, while ML models (SVR and LMDN) did not ($\overline{\mathit{SSPB}}$ = −7.2 ± 5.2%), possibly because the latter uses at least four additional MSI spectral bands.

Stratified performance

Analysis of stratified performance (OWT1 vs. OWT2) suggests that SVR significantly outperformed all other algorithms in the southern basin, which is almost 80% of the lake area (Table 4). SVR also excelled relative to other algorithms in the northern basin, when considering most performance metrics, including MdSA and MWR. However, LMDN performance was comparable to that of SVR in the northern basin and even surpassed it in terms of SSPB and Slope. The reasonable performance of LMDN with few data (e.g., in the northern basin with only 45 samples for training) was unexpected; however, results should be treated with caution due to the low availability of test data (11 samples in each run; Table D1). Scatter plots in Figure D1 (Appendix D) further demonstrate that empirical models failed to estimate Chla in the northern basin (Slope < 0.1), while LML models provide better estimates of Chla in turbid water (Slope = 0.3).

Table 4. Evaluation metrics for Chla retrieval models on MSI and in situ Chla matchups based on water type.

		MdSA (%)	SSPB (%)	RMSLE	MAPE (%)	RMSE (mg m^−3)	Slope	MWR (%)
OWT1 N = 137 $\overline{\mathit{Chla}}$ = 26 $\overline{\mathit{TSS}}$ = 7	OC3	108.7	−14.4	0.46	64.55	34.99	0.04	80.3
	2band	36.84	−4.27	0.25	29.55	16.49	0.72	54.7
	3band	46.93	−5.84	0.31	33.28	30.19	0.64	56.9
	NDCI	37.7	−5.6	0.26	30.69	22.19	0.76	56.2
	MDN	50.4	29.9	0.28	39.6	29.04	1.19	59.1
	LMDN	38.8	6.4	0.25	29.91	15.34	0.81	52.1
	SVR	32.7	4.2	0.24	29.57	15.43	0.82	N/A
OWT2 N = 56 $\overline{\mathit{Chla}}$ = 15 $\overline{\mathit{TSS}}$ = 15	OC3	59.3	−10.5	0.39	39.86	10.75	0.01	60.7
	2band	84.3	−22.4	0.81	48.02	12.33	−0.16	66.1
	3band	69.2	−20.3	1.95	46.16	12.27	−0.15	66.1
	NDCI	74.6	−16.6	0.4	47.45	10.51	0.04	60.7
	MDN	87.8	10.9	0.48	55.02	51.04	2.6	66.1
	LMDN	53.1	3.2	0.34	39.9	8.94	0.34	51.6
	SVR	51.4	−6.9	0.33	38.04	8.85	0.29	N/A

See Figure A1 and section Study site for more details of each water type and the classification approach. The Model Win Rate (MWR) is computed relative to SVR as the reference model. $\overline{\mathit{Chla} }$ and $\overline{\mathit{TSS}}$ are the median of Chla and TSS in associated stations. Units are mg m⁻³ and g m⁻³ for Chla and TSS, respectively. Highlighted cells mark the best score for the corresponding metrics in each OWT.

Except for OC3, Chla retrieval was more accurate (15–50% improvement in MdSA) in the southern basin compared to the northern site. The higher concentration of suspended sediments and NAP in the northern basin, which leads to a higher Chla interference by NAP backscattering particularly at longer wavelengths (red-NIR), likely explains the lower accuracy of Chla retrieval at that location. This pattern may also explain the higher accuracy of OC3 in the northern basin; given that it was the only model that did not use red-NIR bands.

Model sensitivity to AC and radiometric products

Model performance was assessed over two different AC processors (ACOLITE, iCOR) and three radiometric products ($R_{rs}^{\delta },$ ${\mathrm{\rho }}_{rc},$ and ${\mathrm{\rho }}_{\mathit{TOA}}$) applied to MSI data (Figure 3). While ACOLITE provided all three products, iCOR only returns $R_{rs}^{\delta }.$ Overall, SVR and LMDN manifested robust outputs for both AC processors and all the radiometric products ($\overline{\mathit{MdSA}}$ = 43.7 ± 3.7%). In contrast, the mean of variability for empirical models was almost 2-fold greater (±7.8%) than these models, with a maximum for OC3 (±14.8%) and a minimum for 2band (±3.5%). SVR-$R_{rs}^{\delta ,\mathit{ACL}}$ exhibited the best performance of all combinations of retrieval models and AC processors. SVR’s superiority was also evident when employing $R_{rs}^{\delta ,\mathit{iCOR}}$ or ${\mathrm{\rho }}_{rc},$ with only ${\rho }_{\mathit{TOA}}$ showing comparable results to those obtained with LMDN (<2% difference).

Figure 3. Median Symmetric Accuracy (MdSA) for Chla retrieval algorithms when applied to MSI-A/B data processed to produce different radiometric products ($R_{rs}^{\delta },$ ${\mathrm{\rho }}_{rc},$ and ${\mathrm{\rho }}_{\mathit{TOA}}$) with different AC processors (ACOLITE and iCOR). Note that ${\mathrm{\rho }}_{rc}$ is generated with ACOLITE and theoretically is not different when using iCOR. N is the total number of matchups. See Table D1 for the detailed training/test split process.

No single AC processor or radiometric product performed best in all Chla retrieval models. For example, OC3 and 3band worked better with iCOR as the AC processor, while the others (2band, NDCI, LMDN, SVR) all presented better results with ACOLITE. For these latter models, $R_{rs}^{\delta }$ displays the highest accuracy compared to the other products (${\mathrm{\rho }}_{rc},$ ${\mathrm{\rho }}_{\mathit{TOA}}$), suggesting that ACOLITE outperformed iCOR whenever it successfully carried out aerosol correction (${\mathrm{\rho }}_{rc}$→$R_{rs}^{\delta }$). Our results also show that Rayleigh correction (${\mathrm{\rho }}_{\mathit{TOA}}$→${\mathrm{\rho }}_{rc}$) as implemented in ACOLITE reduced Chla retrieval accuracy except for OC3, confirming that this procedure over-corrects reflectance in red-NIR wavelengths while remaining suitable for use with blue-green bands. On the other hand, declining accuracy after aerosol correction in OC3 applications indicates that the AC processors failed to accurately remove aerosol effects in blue-green bands, a task that has proven to be challenging elsewhere (Pahlevan et al. 2021).

Model transferability over water types

Model transferability over two OWTs in BPL was assessed using $R_{rs}^{\delta ,\mathit{ACL}}$- Chla matchups derived from MSI images (see section Model assessment) (Figure 4). All empirical algorithms (OC3, NDCI, 2band, and 3band) failed to retrieve Chla when they were trained by matchups from a different, but similar, OWT (MdSA > 100%, Slope < 0.2). Additionally, LMDN showed poor transferability over both water types (MdSA > 200%, Slope < 0.2). In contrast, SVR maintained a reasonable transferability over two OWTs (MdSA = 61%, Slope = 0.35) compared to alternate models. Although the error and bias increased ∼2- to 4-fold compared to instances where both OWTs were used to train the SVR model (MdSA = 61 vs. 36% and SSPB = 15.8 vs. 3.4%) (see section General performance), they remained within an acceptable range for many applications. SVR’s high transferability might be related to its proven resistance to overfitting, thanks to the regularization parameter C.

Figure 4. Scatter plot of in situ Chla versus predicted Chla from MSI-A/B images. Chla values in the northern basin (OWT2, red solid circles) are predicted using a model trained with southern basin matchups (OWT1, blue solid circles) and vice versa.

Model sensitivity to sensor type

Matchups of $R_{rs}^{\delta ,\mathit{ACL}}$–Chla derived from OLI images were employed to retrieve Chla in BPL. LMDN outperformed SVR in most metrics when using OLI data, by ∼5% in MdSA and with a 2-fold greater Slope (Table 5). MDN displayed an overall error of 95% and a bias of ∼50% reflecting the training of this global model with in situ $R_{rs}$ rather than $R_{rs}^{\delta }.$ Additionally, even though OLI lacks spectral bands at red-edge wavelengths, a red-NIR empirical model (FLH-blue) outperformed the blue-green-based index of OC3 by ∼10%. A global BST model (Cao et al. 2020) failed to estimate Chla in BPL (results not shown here), similar to what has been observed elsewhere (Smith et al. 2021), likely due to much lower CDOM absorption in BPL compared to the waterbodies that were used to train BST ($\overline{a_{\mathit{CDOM}}^{(440)}}$ = 0.28 m⁻¹ vs. $\overline{a_{\mathit{CDOM}}^{(440)}}$ = 0.8–1 m⁻¹). Finally, the LBST model exhibited poor performance (MdSA = 77%), possibly because the boosting algorithms degrade in the presence of outliers and errors in training data (Li and Bradic 2018).

Table 5. Evaluation metrics for Chla retrieval models on OLI and in situ Chla matchups (N = 178).

	MdSA (%)	SSPB (%)	RMSLE	MAPE (%)	RMSE (mg m^−3)	Slope	MWR (%)
OC3	75.7	3	0.35	56.41	31.65	0.08	0.57
FLH-b	64.7	−6.6	0.33	45.94	29.05	0.2	0.6
LBST	77.2	2.2	0.35	53.08	30.06	0.19	0.59
MDN	94.98	48.9	0.42	64.7	38.4	0.14	0.65
LMDN	55.02	8.1	0.32	40.23	27.43	0.45	0.47
SVR	60.27	2.8	0.32	43.78	27.83	0.23	N/A

Each model was trained and tested using a 5-fold cross-validation approach. The MWR was computed relative to SVR as the reference model. Highlighted cells mark the highest score for the corresponding metrics.

Overall, Chla retrieval using OLI data (Table 5; $\overline{\mathit{MdSA}}$ = 71.3 ± 13.2%) appeared less accurate than that based on MSI summarized in Table 5 ($\overline{\mathit{MdSA}}$ = 55.2 ± 19.3%). OLI’s poor performance was also inferred from low Slope (< 0.5), likely due to the absence of a red-edge band. Similar to MSI, LML models exhibited better performance than empirical and GML models when applied to OLI data. The analysis of scatter plots (Figure D2) also revealed that all models failed to estimate Chla values <10 mg m⁻³ and concentrations >100 mg m⁻³. Although the former limitation was also observed when using MSI data (see section General performance), the latter might be intensified because OLI does not possess a spectral band in the domain of Chla fluorescence (680–710 nm).

Spatial integrity

Chla maps for BPL were generated from an MSI image taken on July 16, 2020 (Figure 5). All model-processor combinations suggested Chla as low as ∼10 mg m⁻³ in the north basin, whereas some models/processors (e.g., SVR-iCOR) predicted Chla values up to ∼100 mg m⁻³ in the south basin. Regardless of the AC processor used, ML models (SVR and LMDN) seem to deliver overall smoother maps (less noise) compared to the 2band output, probably due to leveraging all spectral bands.

Figure 5. Chla maps for BPL derived from different retrieval algorithms/AC processors couples applied on MSI-A image acquired on July 16, 2020. The markers in the insets represent examples of the location of in situ data, collected on the same date, and employed as unseen test data. The color bars and associated numbers beside the markers show estimated Chla concentration in mg m⁻³. In situ Chla concentration in points A, B, and C are 41.2, 66.2, and 102.8 mg m⁻³, respectively. 2band was used as the best representative of empirical models.

Visual comparison of Chla maps based on near-coincident in situ measurements revealed that the SVR model, coupled with iCOR processor, had the highest consistency with in situ measurements (Figure 5). Although all models/processors showed a reasonable and similar performance in mapping moderate Chla concentrations (Figure 5, upper insets), they differed more substantially in estimating high Chla values at the south of the lake. SVR tended to estimate higher Chla concentrations than did LMDN and 2band models, regardless of AC processors (lower insets in Figure 5). SVR-iCOR also seemed to be more capable of detecting high spatial gradients in Chla, as it is the only combination to capture large gradients of Chla at two nearby stations (Chla = 66.2 to Chla = 102.8; lower insets Figure 5). Such high-frequency changes in Chla may be related to the surface patchiness of cyanobacteria.

SVR results appeared prone to mixed pixels compared to LMDN and 2band models. Although this effect was limited to 1–2 pixels close to the shoreline, this issue should be treated with caution when producing maps of nearshore Chla. Similarly, despite being very eutrophic ($\overline{\mathit{SDD}}<1 \mathrm{m}$), SDD measurements across the lake show that a small portion of the lake area in the north basin can be considered as optically shallow waters, mostly in very early or late summer (i.e., May or October). Consequently, the elevated Chla estimates produced for the northern basin by models are probably influenced by very shallow depths (<2 m) or a high density of rooted aquatic macrophytes. While maps were produced using MDN and empirical models, none outperformed the above-mentioned models. For instance, MDN returned some unrealistically high Chla values, and OC3 routinely and significantly underestimated Chla.

As it is sometimes more important to reconstruct spatial patterns of Chla than accurately estimate absolute concentrations, we normalized the predicted Chla vector of unseen matchups for stations 4–11 (a longitudinal transect along the lake), by dividing by the vector norm to better evaluate which algorithms recorded spatial patterns of Chla in BPL (Figure 6). Overall, normalization did not reveal a single superior model/processor in terms of retrieving spatial gradients of Chla. While SVR-iCOR provided the most similar pattern to measured Chla gradients in the northern basin (#station > 8), SVR-ACOLITE demonstrated good performance in retrieving Chla changes in the southern stations 5–8. In contrast, the 2band model performed well at stations 4–5 whereas LMDN performed poorly at stations 4–6 and 10–11. Together, these patterns suggest that SVR showed the highest overall capability in retrieving the Chla gradient along the lake.

Figure 6. Spatial profile of normalized Chla along the lake (south to north) for July 16, 2020, derived from in situ measurement Chla (solid line) as well as predicted Chla from algorithms applied on MSI image (dashed lines). X-axis denotes station number (see Figure 1).

We also mapped Chla over the lake using OLI data for the same date (July 16, 2020) using FLH-Blue, LBST, LMDN, and SVR models (Figure 7). Maps from LBST and LMDN were markedly noisy, whereas LMDN showed reasonable quantitative performance for OLI data (Table 5), and FLH-Blue and SVR generated smooth maps. The SVR model exhibited more consistency with in situ data (marked points in Figure 7), while LMDN retrieved Chla values higher (120 mg m⁻³) than observed in situ, and the other algorithms underestimated Chla. In terms of reconstructing the spatial pattern of Chla, LMDN seems to provide the best performance, consistent with its higher Slope (Slope = 0.45) (Table 5). Appendix E also shows more examples of the produced Chla maps for BPL as well as a pixel-by-pixel comparison of MSI- and OLI-derived Chla in same-date images over BPL.

Figure 7. Chla map for BPL derived from different algorithms applied on OLI image acquired on July 16, 2020. The markers in the insets represent examples of in situ data, collected on the same date, and employed as unseen test data. The color bars and associated numbers beside the markers show estimated Chla concentration in mg m⁻³. In situ Chla concentration in points A, B, and C are 41.2, 66.2, and 102.8 mg m⁻³, respectively.

Temporal validity

Robust retrieval of Chla over time is a daunting task in a eutrophic waterbody due to high variations in surface bloom densities, resultant water optics, and atmospheric conditions. Comparison of SVR-iCOR, SVR-ACOLITE, LMDN-ACOLITE, and 2band-ACOLITE processing couples at station 1 in BPL revealed that SVR-iCOR tracked in situ Chla measurements better than the other model/processor combinations (Figure 8), with particularly good capture of intense summer blooms in July to September. Although none of the models accurately captured the peak of Chla (>100 mg m⁻³) over the investigated period, SVR-iCOR followed the shape and magnitude of the measured time series with a ∼15% underestimation of peak Chla values. In contrast, couples based on ACOLITE failed to deliver consistent Chla values on August 7, 2020 when cloud shadow contaminated images. For more moderate Chla concentrations (20–60 mg m⁻³), SVR-ACOLITE displayed better performance than SVR-iCOR. Overall, a correlation analysis between the time series of measured and predicted Chla showed that SVR-iCOR (ρ = 0.798) outperformed other models (ρ = 0.684–0.728) in retrieving Chla time series.

Figure 8. Time series of Chla in station 1 in BPL for summer 2020, derived from in situ measurement Chla (solid line) as well as predicted Chla from algorithms applied to MSI images.

Discussion

Analysis of Landsat 8 and Sentinel 2 images using locally trained machine-learning models, particularly those based on SVR, provided robust retrieval of Chla for a small eutrophic lake using MSI and OLI images (sections General performance and Stratified performance). These models also generated realistic annual time series and spatial gradients of Chla of scales appropriate to the prairie lake (sections Spatial integrity and Temporal validity). Overall, these models were robust to variations in AC processors (ACOLITE vs. iCOR) and sensor types (MSI vs. OLI). Together, our analysis suggests that pre-trained SVR models may provide important information on spatial and temporal patterns of water quality and HABs in regional lakes, provided that optical water types and atmospheric conditions are similar. However, the fact that the results here are based on a single lake study may necessitate further investigations of the presented model in other regional lakes.

Uncertainties in Chla and radiometric data

Although we attempted to reduce the uncertainties associated with in situ Chla data, any comparison of remotely sensed images and discrete lake measurements can be complicated due to the high variability of in situ data (Clay et al. 2019; Qiu et al. 2021). Here, we tried to reduce random noise in Chla measurements by conducting each measurement several times and averaging values. However, our in situ data originated from different laboratories using contrasting measurement techniques (field fluorometry, laboratory spectrophotometry, HPLC), instrumentation, calibration, and field sampling (surface 1 m vs. depth-integrated). While these factors may affect model performance, they also suggest that our algorithms exhibit minimal overfitting and systematic errors in performance assessment, and may be generalizable to other regional lakes.

Several lines of evidence suggest that potential outliers and other uncertainties in Chla measurements did not alter the results of comparative assessment of retrieval models. First, we used the median symmetric accuracy (MdSA) as the main metric to compare the models, as it is highly robust to potential outliers in in situ Chla measurements. Second, we conducted various experiments with different numbers and combinations of matchups, and in all cases, SVR showed robust and similar results, meaning that uncertainties in lake production do not substantially alter results. Moreover, given that BPL is well mixed vertically (Dröscher et al. 2008), we expect that differences in sampling protocols may not greatly affect our findings. Finally, earlier studies suggest that SVR can handle diverse, highly uncertain datasets because they use only a part of the data (support vectors) for learning (Chegoonian et al. 2017; Chegoonian et al. 2021; Foody and Mathur 2006; Hu et al. 2021; Nikparvar and Thill 2021). Handling uncertainties of in situ data becomes crucial when input data to observatory systems originate from diverse field and laboratory sources.

Comparable accuracy of Chla retrieval accuracy obtained from 10, 20, and 60 m MSI data (Figure B3) supported the use of 60 m data, which substantially (>5 times) reduced the time needed to develop a reliable model. This finding is consistent with studies that utilized 60 m resampled data (Ansper and Alikas 2018) and averaged pixel values over an equal-sized window (Pahlevan et al. 2020; Werther et al. 2022). Although the use of 60 m data may increase the likelihood of having mixed pixels in our model, lower spatial resolution also reduces random noise arising from very fine resolution (10–20 m) imaging of aquatic environments. Moreover, our sampling stations (Figure 1) were on the central axis of the lake, at least 500 m from shore, which eliminates the possibility of mixed aquatic-terrestrial pixels even with 60 m data. Our field observations of bloom formation in the lake during the summers of 2014–2021 also indicate that very small scale patchiness of phytoplankton blooms are rare, consistent with the similar performance of models based on 10, 20, and 60 m MSI data.

Merits of locally trained ML models

When compared with traditional empirical models (e.g., OC3, 2band), LML models exhibit several clear advantages, particularly with regard to SVR models. First, their ability to leverage all spectral bands and the capability to learn and model diverse uncertainties (in situ data, non-linearity, non-Chla constituents) is an advantage over traditional empirical/physical models and led to 15–65% error reduction. Such performance might be improved further when using models, such as LMDN that can deal with ill-posed problems (Pahlevan et al. 2020).

Currently, the uncertainties in AC processors are the major hurdle for employing GML models in inland waterbodies (Pahlevan et al. 2020). These models are often trained with in situ radiometric measurements and can be degraded when fed by satellite-derived measurements. LML models that can learn AC uncertainties ($R_{rs}^{\delta }$) specific to a lake of interest may be an important solution for application to local and regional resource management issues, such as blooms of toxic cyanobacteria near recreational areas or drinking water inlets. Meanwhile, the development of global models based on satellite-derived reflectance or including ancillary data may provide an opportunity to expand the geographic range of applications of ML models (Smith et al. 2021).

Presently, the need for substantial training data is a major obstacle to the development of local ML models. Fortunately, here we demonstrate that LML models (SVR and LMDN) were trainable with ∼200 matchups (section General performance), while the stability of the results even with only ∼50 matchups was encouraging (section Stratified performance), as many regional agencies in Europe and North America conduct routine monitoring (e.g., Soranno et al. 2017). Ideally, such locally trained models should possess reasonable generalization to retrieve reliable Chla in nearby lakes where optical conditions, water type, and atmospheric conditions differ only slightly. Our results suggest that SVR models exhibit adequate transferability when trained and tested with two different (but similar) water types in BPL (section Model transferability over water types). Although this capability is in agreement with SVR resistance to overfitting (Kwiatkowska and Fargion 2003; Mountrakis et al. 2011; Zhan et al. 2003), it is still essential to further validate our results using a more consistent and systematically collected/calibrated in situ Chla dataset.

This study was also the first independent assessment of the global MDN model in a small eutrophic lake. Although MDN is not expected to outperform locally trained models, it showed errors within ∼60% of in situ measurements. Nonetheless, MDN tended to significantly overestimate Chla (high bias) relative to locally trained, $R_{rs}^{\delta }$-fed models. Substantial uncertainties in AC process, which can be seen in drastically different $R_{rs}^{\delta }$ distributions from ACOLITE and iCOR (Figure B4), or low performance of the model with respect to spectral ambiguities, may explain MDN overestimation.

Atmospheric correction

Algorithms developed to retrieve downstream products, such as Chla, always should exhibit consistent performance with different intermediate processors, specifically AC processors. Here, we demonstrate the robustness of the SVR model when data is processed using ACOLITE and iCOR, and three different radiometric products ($R_{rs}^{\delta },$ ${\mathrm{\rho }}_{rc},$ and ${\mathrm{\rho }}_{\mathit{TOA}}$). The fact that ${\mathrm{\rho }}_{rc},$ and ${\mathrm{\rho }}_{\mathit{TOA}}$ exhibit reasonable results—especially when using red-NIR bands—is in agreement with findings from previous studies (Matthews et al. 2012; Matthews and Odermatt 2015; Wynne et al. 2010) and can support using atmospherically-uncorrected data commonly-available on global scale (e.g., Google Earth Engine). Furthermore, our results show that the accuracy of Chla estimates was generally greatest when using $R_{rs}^{\delta }$ for retrieval models, other than those based on OC3 and 3band for which ${\mathrm{\rho }}_{rc}$ generated more accurate products. We also note that empirical algorithms using blue-green bands (e.g., OC3) significantly benefited from Rayleigh correction for blue-light scattering. While Rayleigh correction did not appear to increase the accuracy of the models that were based on red-NIR bands, further evaluations are needed to evaluate this finding.

Modeling results were consistent with those of Pahlevan et al. (2021), who recently conducted a comprehensive comparison between AC processors in retrieving $R_{rs}$ using an extensive global dataset. For example, we observed that 2band and NDCI—two algorithms that use only 665 and 704 nm bands—performed better when they are coupled with ACOLITE than with iCOR. Similarly, OC3 and 3band models that use blue (443 or 492 nm) and 740 nm bands showed better performance with iCOR when compared to ACOLITE (Figure 3). We interpret the high consistency between the assessments of downstream products (Chla concentration) and satellite-derived reflectance as an indicator of the effectiveness of AC process on the accuracy of downstream products. However, we also recognize that further examination of the effectiveness of AC will require a separate estimate of retrieval uncertainty from the AC process; an assessment that needs field radiometric measurements which were not sufficiently available in our study.

Comparisons among experiments in this study, as well as drastically different $R_{rs}^{\delta }$ distributions derived from ACOLITE and iCOR (Figure B4), suggest that different AC processors may lead to significant differences in retrieval performance. Thus, the algorithms for retrieval of downstream products should be examined as retrieval models/AC processors. For instance, the SVR model shows greater accuracy when used in conjunction with ACOLITE (Figure 3), and more temporal stability when using iCOR as the AC processor (Figure 8). However, the comparison between ACOLITE and iCOR is not entirely equivalent due to differences in the number of matchups (15 more for iCOR when masked by ACOLITE); thus, other studies (Ilori et al. 2019; Pahlevan et al. 2021; Warren et al. 2019; Xu et al. 2020) are needed for a more comprehensive comparison of AC processors.

Conclusion

This paper presents a machine-learning model based on support vector regression (SVR) to retrieve Chla concentration from satellite-derived reflectance measurements ($R_{rs}^{\delta }$) of Sentinel-2 (MSI) and Landsat-8 (OLI). The proposed model was trained and evaluated using a dataset of near-coincident, co-located in situ Chla and $R_{rs}^{\delta }$ observations (N ∼ 200), collected in a mid-latitude eutrophic lake from 2014 to 2020. Comparison of the SVR model against state-of-the-art, commonly used alternates revealed that SVR outperformed all other algorithms when using MSI data. This superiority is seen in both general (entire samples, Chla = 1–125 mg m⁻³) and stratified levels (two distinct optical water types).

The proposed model also showed superiority in retrieving time series of Chla and producing Chla maps, two important applications of remote sensing in monitoring and mapping of harmful algal blooms. The superiority of SVR was also demonstrated by the return of robust and similar results following the alteration of AC processors (ACOLITE vs. iCOR). The model was also stable when fed with different radiometric products ($R_{rs}^{\delta },$ ${\rho }_{rc},$ and ${\rho }_{\mathit{TOA}}$). Quantitative evaluation of SVR also showed a promising transferability among two optical water types common to this study region, particularly in comparison to standard models.

Together, these findings reveal the high potential of SVR models to retrieve Chla in small waterbodies, even using data from multi-spectral terrestrial missions, such as MSI and OLI. Although results are presented only for BPL, the fact that the lake is broadly representative of over 100 regional lakes within a 240,000 km² area (Finlay et al. 2015; Hayes et al. 2020) suggests that our findings may be generalized to other eutrophic mid-latitude waterbodies of similar optical water types. Development of such models for consistent retrievals from long-term observational records of satellite missions, such as Landsat and Sentinel increases the potential for monitoring and mapping the extent and intensity of harmful algal blooms in an era of global warming.

Author contributions

Conceptualization: AMC; limnological in situ data: PRL, HMB, and JMD; method development: AMC, NP, and KZ; data analysis: AMC; manuscript preparation: AMC; editing and approval: All authors; funding: CRD, PRL, HMB, and NP.

Acknowledgments

We thank members of University of Saskatchewan (US) Global Institute for Water Security and the University of Regina (UR) Limnology Laboratory for field data collection, as well as David Vandergucht and seasonal staff with the Saskatchewan Water Security Agency, and the Buffalo Pound Water Treatment Plant (BPWTP). Operations of the monitoring buoy were supported by Jay Bauer, Katy Nugent, Cameron Hoggarth, and staff of BPWTP. We gratefully acknowledge that the field research and UR analyses took place on Treaty 4 territory, homelands of the Cree, Saulteaux, Lakota, Dakota, and Nakota peoples, as well as the Metis/Michief nation. US is located on Treaty 6 territory, while University of Waterloo is located on the traditional territory of the Neutral, Anishinaabeg, and Haudenosaunee peoples.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

Buoy operations and associated research supported by grants to HMB from the Natural Sciences and Engineering Research Council of Canada (NSERC), Canada Foundation for Innovation (CFI), Global Water Futures (GWF)-Canada First Research Excellence Fund (FORMBLOOM project), Global Institute for Water Security, and Buffalo Pound Water Treatment Plant. Qu’Appelle Long-term Ecological Research program (QU-LTER) was supported by grants to PRL from NSERC, CFI, Canada Research Chairs, the Province of Saskatchewan, and University of Regina. NP was supported under the NASA ROSES contract #80HQTR19C0015, Remote Sensing of Water Quality element, and the USGS Landsat Science Team Award #140G0118C0011. AMC, CRD and KZ were supported by GWF TTSW (Transformative sensor Technologies and Smart Watersheds for Canadian Water Futures) project.

Appendix A

 

Chla in BPL is higher than the average for freshwaters (Filazzola et al. 2020), as opposed to Dissolved Organic Carbon (DOC) which is in a low/medium amount for freshwaters (Toming et al. 2020). It can be claimed that particles, especially algal particles, mostly control the optical characteristics of water in BPL. This hypothesis can be confirmed by Figure A1 (upper diagonal) where the optical characteristics of water samples e.g., turbidity and Secchi Disk Depth (SDD) are highly correlated with Total Suspended Solids (TSS) and Chla. However, Figure A1 (lower diagonal) also reveals that the relationship depends on the station location as southern stations (1–8) show a stronger relationship between Chla and optical characteristics. This can be justified based on the distribution of water constituents, plotted in Figure A1 (diagonal), where northern stations (9–11) contain more sediments, as opposed to the southern stations which are more dominated by Chla.

Figure A1. Pair plots of some optically-derived/driven parameters in BPL (averaged on stations 4–11 from late May to early September  of 2014–2020). Diagonal elements are the distribution of each parameter, color-coded according to station numbers. Upper-diagonal elements are the scatter plot of paired parameters. Lower-diagonal charts are the contour plots showing the relationship between the parameters in northern and southern stations. N and ρ are the number of samples and correlation coefficients, respectively. Units are mg m⁻³, g m⁻³, NTU, and m for Chla, TSS, turbidity, and SDD, respectively.

Appendix B

 

Figure B1. Comparison of MSI bands (red boxes) and OLI bands (blue boxes) in wavelengths <800 nm. The spectra are from three different samples measured at BPL using an ASD spectrometer and display how water spectra change with changes in Chla content.

Figure B2. Overview of a workflow developed in this study.

Figure B3. Scatter plot of in situ Chla versus predicted Chla from MSI-A/B images, resampled to 10, 20, and 60 m. Data from Buffalo Pound Lake 2017–2021. Model performance does not vary consistently with image resolution.

Figure B4. Normalized frequency distributions of MSI-derived $R_{rs}^{\delta }$ spectra for the matchups processed via ACOLITE (N = 193) and iCOR (N = 208) processors.

Appendix C

 

SVR uses an ε-insensitive cost function (with ε as a threshold) in which errors ($e_i)$ up to ε are not penalized, whereas larger errors are penalized using a linear function, i.e., $L\left(e_i\right)=\max (|e_i|-\epsilon , 0).$ Thus, compared to traditional ML models (e.g., multilayer-perceptron neural networks), SVR is more robust to small errors and inherent uncertainties in training data (Zhan et al. 2003). Weights (ω) are estimated in the linear regression problem (Equation A1(A1) $${\mathit{Chla}}_i=\sum _{j=1}^M{\omega }^T\varphi \left(R_{{rs}_i}^{\delta }\right)$$(A1) ), where $i$ is the number of training data, $j$ is the number of predictors (spectral bands), and $\varphi$ is a kernel (a non-linear mapping function). SVR minimizes Equation (A2)(A2) $$\mathit{Cost}= \frac{1}{2}{\Vert \omega \Vert }^2+C\sum _{i}L\left(e_i\right)$$(A2) , where ${\xi }_i$ are $\left|e_i\right|>\epsilon$ and C is the regularization parameter, balancing the minimization of errors and generalization capabilities (Camps-Valls et al. 2006; Smola and Schölkopf 2004). Figure C1 depicts a schematic view of the regression between Chla and reflectance measurements using SVR. (A1) $${\mathit{Chla}}_i=\sum _{j=1}^M{\omega }^T\varphi \left(R_{{rs}_i}^{\delta }\right)$$(A1) (A2) $$\mathit{Cost}= \frac{1}{2}{\Vert \omega \Vert }^2+C\sum _{i}L\left(e_i\right)$$(A2)

Figure C1. Graphical depiction of principles of support vector regression (SVR). (a) Schematic view of regression between Chla and $R_{rs}^{\delta }$ using SVR. (b) Loss function defined for SVR; while errors less than $\epsilon$ are not penalized, larger errors are penalized by a linear function.

Appendix D

 

Figure D1. Matchup analysis of measured and predicted Chla from in situ Chla and MSI-A/B images for two different regions in BPL, categorized based on optical water type. For each optical water type, a model is trained and tested using a 5-fold cross-validation approach.

Figure D2. Matchup analysis of Chla derived from different algorithms applied on OLI data and near-coincident, co-located in situ Chla samples in BPL. The results are from a 5-fold cross-validation approach.

Table D1. Assessment approaches (training-test split) as well as the number of $R_{rs}^{\delta }$–Chla training/test matchups available for each experiment in this study.

Experiment	Evaluation approach	$N_{\mathit{Total}}$	$N_{\mathit{Training}}$	$N_{\mathit{Test}}$
General performance	Cross-validation on yearly-grouped matchups	193	∼150	∼50
Stratified performance	5-fold cross-validation	137*	110*	27*
		56**	45**	11**
Model sensitivity to AC processors	5-fold cross-validation	193^†	154^†	39^†
		208^‡	166^‡	42^‡
Model transferability over water types	Cross-validation on OWT-grouped matchups	193	137	56
			56	137
Model sensitivity to sensor type	5-fold cross-validation	178	142	36
Spatial integrity	Matchups for a specific date are set apart as test data	193	184	9
Temporal validity	Station 1 matchups collected in 2020 are set apart as test data	193	173	20

*Matchups for OWT1.

**Matchups for OWT2.

^†Matchups for ACOLITE.

^‡Matchups for iCOR.

Table D2. Formulas and coefficients of empirical models employed in this study. b# and w# are the reflectance and wavelength at specified bands, respectively.

Sensor	Algorithm	Formula
MSI	OC3 (O'Reilly et al. 1998)	$y=2.198-2.404x,$ $x=lo[\max (b1,b2)/b3]$
	2band (Moses et al. 2009)	$y=-1.165+4.146x-{1.337x}^2,$ $x=\log [\mathrm{b}5/b4]$
	3band (Dall’Olmo and Gitelson 2005)	$y=1.094+1.631x-{0.65x}^2,$ $x=\log [({b4}^{-1}-{b5}^{-1})\mathrm{\hspace{0.17em}*\hspace{0.17em}}\mathrm{b}6]$
	NDCI (Mishra and Mishra 2012)	$y=1.082+2.283x,$ $x=\log [(\mathrm{b}5\text{-}\mathrm{b}4)/(b5+b4)]$
OLI	OC3 (O'Reilly et al. 1998)	$y=0.807-1.886x,$ $x=\log [\max (b1,b2)/b3]$
	FLH-blue (Beck et al. 2016)	y $=0.799+98.256x,$ $x=b3-(b2+(b4-b2\left)\mathrm{\hspace{0.17em}*\hspace{0.17em}}\right(\frac{w3-w2}{w4-w2}))$

y (output of the equations) is Chla on a logarithmic scale.

Table D3. Annual frequency and statistics of $R_{rs}^{\delta }$–Chla matchups derived from MSI sensor.

	2017	2018	2019	2020
N	33	51	48	61
Mean Chla	17.3	19.6	18.4	20
Median Chla	15.4	26	18.8	19.3
Standard deviation Chla	2.8	2.6	4.3	2.3

Chla concentration unit is mg m⁻³.

Table A1. Statistics for water constituents associated with two distinct OWTs in BPL (averaged from late May to early September  2017–2020).

	OWT1 (stations 1–8)		OWT2 (stations 9–11)
	Chla	TSS	Chla	TSS
Median	17.1	6.5	11.4	11
Mean	25.9	7.2	14.9	14.8
Standard deviation	24.7	4.1	12.2	13.1
N	159	104	201	125

Comparing the values in this table with values presented elsewhere (Pahlevan et al. 2021, Spyrakos et al. 2018) revealed two distinct OWTs in BPL. Units are mg m⁻³ and g m⁻³ for Chla and TSS, respectively.

Appendix E

 

Figure E1. Chla maps for BPL derived from SVR algorithm applied on MSI (left) and OLI (right) images acquired on the same dates from 2017 to 2020. The color bars show the estimated Chla concentration in mg m⁻³.

Figure E2. Matchup analysis of Chla derived from SVR algorithm applied on the three same-date MSI and OLI images captured over the BPL from 2017 to 2020.