Remote Sensing of Chlorophyll-A in Clear vs. Turbid Waters in Lakes

1. Introduction

Recent increases in the frequency, prevalence, and intensity of harmful phytoplankton blooms in lakes worldwide have raised significant environmental concerns [1,2,3]. These blooms contribute to issues such as toxin production, water odor, taste degradation, and aesthetic decline [4,5]. They also reduce biodiversity through habitat degradation and alter energy and nutrient fluxes within food webs [5,6,7]. Addressing these challenges requires a thorough understanding of past, current, and potential future occurrences of phytoplankton blooms, necessitating regular, lake-wide monitoring of phytoplankton biomass [8]. However, traditional in situ monitoring methods are costly, time-consuming, and limited in spatial and temporal coverage [9].

Satellite sensors offer a cost-effective solution for monitoring chlorophyll-a (Chl-a), a proxy for phytoplankton biomass, across lakes by mapping its concentration [10]. Estimating Chl-a from satellite spectral data relies on the absorption and scattering of light by phytoplankton at specific wavelengths [11]. In waters with straightforward optical properties, phytoplankton and related substances predominantly determine the color of the water and thus the signals received by satellite sensors [12]. However, in optically complex waters like many inland and coastal areas, non-phytoplankton optically active constituents (OACs) such as detritus and colored dissolved organic matter (CDOM) also influence the color of water, complicating Chl-a estimation from satellite data [12]. These OACs vary widely in type, quantity, and band ratios across different water bodies [13], posing challenges for developing universal Chl-a prediction models [14,15].

Successful Chl-a estimation from satellite data requires the development of region-specific models, known as regional models, which account for the optical characteristics of the target water body [11,16,17]. Efforts to create such models often rely on lake-wide approaches that overlook within-basin variations in water optical properties [18].

This study aims to propose a method for developing simple yet accurate Chl-a prediction models that address the optical heterogeneity of lakes. Using Lake Winnipeg (LW), Canada, as a case study, we hypothesize that improving Chl-a predictions can be achieved by separating the lake into its optically distinct basins. We compare "basin-specific" models against conventional "lake-wide" models to evaluate their efficacy. Enhanced Chl-a predictions using basin-specific models could facilitate Chl-a modeling on platforms like Google Earth Engine (GEE), which is crucial for processing and analyzing long-term satellite data at broad spatial scales using medium- to high-resolution imagery, often computationally intensive outside cloud-based environments.

2. Materials and Methods

2.1. Case Study Area

Lake Winnipeg, with a surface area of 23,750 km², is recognized as the world’s 10th largest, North America’s 7th largest, and Canada’s 3rd largest freshwater lake [19,20]. It is a shallow, wind-mixed, polymictic lake, divided into two physically and biogeochemically different basins: the north basin (NB) and the south basin (SB) (Figure 1) [19]. Compared to the SB, the NB is larger, encompassing 74% of the total lake area, and deeper, with an average depth of 13.3 m and a maximum depth of 18.0 m. The water in the NB is relatively transparent. In contrast, the SB, encompasses 11% of the total lake area and is shallower, with an average depth of 9.0 m and a maximum depth of 12.0 m. The Red River transports a large amount of sediments to the SB, combined with sediment resuspension, making the SB highly turbid. A deep, narrow channel, known as the Narrows, separates the NB and SB. The Narrows is characterized by its relatively rapid water flow, with a maximum depth of approximately 60.0 m and a width of around 2.6 km [19].

2.2. LW Turbidity

We used in situ turbidity data (Nephelometric Turbidity Units; NTU) measured throughout LW, provided by Manitoba Agriculture and Resource Development during the growing seasons (May to October) from 2002 to 2018 (Figure 1). A non-parametric, one-way analysis of variance by ranks (Kruskal–Wallis H test) was employed to determine whether the two LW basins (NB and SB) and the Narrows differed in their turbidity levels [21]. The choice of a non-parametric test was based on the non-normal distribution of the data (Shapiro-Wilcoxon test; p > 0.05). A significant result from the Kruskal-Wallis test (p < 0.05) indicated that at least one basin differed from the others. We then applied a multi-step a posteriori pairwise testing procedure based on studentized range statistics (non-parametric multiple comparison Dunn’s procedure) to identify which basins were significantly different from each other [22].

2.3. Modeling Chl-a

2.3.1. In Situ Data Acquisition

In situ Chl-a measurements (μgL⁻¹) were provided by the Manitoba Sustainable Development, Water Quality Management Section in 2019. Chl-a data were collected from 32 stations distributed throughout the lake (14 stations in the NB, 15 in the SB, and three in the Narrows; Figure 1) during the peak phytoplankton biomass season (July to October) between 2002 and 2018 [19]. All Chl-a data were measured using spectrophotometric methods (American Public Health Association (APHA) 10200H [23]). Since satellite-derived Chl-a information is limited to the uppermost water layer, we only used samples taken from the surface of the lake (≤ 0.1 m depth). Although using Chl-a samples from deeper layers (≤ 0.5 m depth) could have substantially increased the number of in situ and satellite matchup points, it reduced the performance of Chl-a models (Figure S1). Our criteria for in situ Chl-a sample selection resulted in no samples from the Narrows, as all samples there were from deeper than 0.1 m. Due to the log-normal distribution of Chl-a in nature [24], all Chl-a data were logarithmically transformed (using natural log) prior to analysis.

2.3.2. Landsat Data

Landsat 4-5 Level 1 Collection 2 Thematic Mapper (TM) and Landsat 7 Enhanced Thematic Mapper Plus (ETM+) images covering the entirety of LW were used in this study (https://www.usgs.gov/landsat-missions/landsat-collection-2). Both sensors (TM and ETM+) have a spatial resolution of 30 m. We used four visible bands: blue (B, 450–520 nm), green (G, 520–600 nm), red (R, 630–690 nm), and near-infrared (NIR, 770–900 nm). Matchups between Chl-a sample times and Landsat image acquisition times were extracted from 9 TM and 20 ETM+ images (6 TM and 6 ETM+ path/rows). Images with more than 90% cloud or cloud shadow cover were excluded.

2.3.3. Pre-Processing Landsat Level 1 Data

Partial atmospheric correction was applied to Landsat Level 1 data (Figure 2) because Level 2 products may introduce higher uncertainty in aerosol calculations over inland waters, leading to erroneous results [25,26]. All pre-processing steps for transforming Landsat data to partially corrected satellite top-of-atmosphere (TOA) reflectance were performed in GEE (Figure 2, Supplementary Material). Land pixels and those covered by cloud or cloud shadow during satellite image acquisition were masked out using a standard LW boundary shapefile and Pixel Quality Assessment bands provided with Landsat Level 2 products. The individual pixels were then pre-processed in three steps.

Step 1. Radiometric correction

Landsat Level-1 data are stored as digital numbers (DNs). During the radiometric correction, we recalibrated the DNs to the satellite TOA radiance (Lλ) using the standard equation [27] (1):

$$L_{\mathsf{\lambda }}=\left({DN}_{\mathsf{\lambda }}\times G_{\mathsf{\lambda }}\right)+B_{\mathsf{\lambda }}$$

(1)

Where:

L_λ is TOA radiance for band λ,

DN_λ is DN for band λ,

G_λ is the multiplicative rescaling factor for band λ, and

B_λ is the additive rescaling factor for band λ.

Step 2. Partial atmospheric correction:

L_λ values were corrected for Rayleigh scattering using a partial atmospheric correction algorithm [28]. We used an inverse algorithm based on a simplified radiative transfer model to calculate Rayleigh scattering radiance (or Rayleigh path radiance) for each band (L_r(λ)) [28] (2):

$$L_{r\left(\lambda \right)}=\left({\displaystyle \frac{{ESUN}_{\lambda }\times cos{\theta }_s\times P_r}{4\pi \times \left(cos{\theta }_s+\mathit{cos}\theta \right)}}\right)\times \left(1-\mathit{exp}\left({-\tau }_r\left(\lambda \right)\times \left(\left({\displaystyle \frac{1}{cos{\theta }_s}}\right)+\left({\displaystyle \frac{1}{\mathit{cos}\theta }}\right)\right)\right)\right)\times t_{oz}\uparrow \left(\lambda \right)\times t_{oz}\downarrow \left(\lambda \right)$$

(2)

Where:

ESUN is the mean solar exo-atmospheric irradiance,

P_r is the Rayleigh scattering phase function,

θ_s is the solar zenith angle (degrees),

θ is the satellite viewing angle (degrees) (set to 0°),

τ_r is the Rayleigh optical thickness, and

t_oz↑ and t_oz↓ are upward and downward ozone transmittance, respectively.

The phase function P_r was calculated using (3):

$$P_r={\displaystyle \frac{3}{4}}\times {\displaystyle \frac{1-\gamma }{1+2\gamma }}\times \left(1+{cos}^2\Theta \right)+{\displaystyle \frac{3\gamma }{1+2\gamma }}$$

(3)

Where:

Θ (the view zenith angle) is the scattering angle (180° − θ_s),

γ = δ/(2- δ) [29], and

δ is the depolarization factor [30], which denotes the polarization of anisotropic particles at right angles and is dependent on the wavelength, atmospheric pressure, and air mass (with the last two variables being constant values) [31].

The Rayleigh optical thickness, τ_r, was calculated using [32,33] (4):

$${\tau }_r=0.008569{\lambda }^{-4}\times \left(1+0.0113{\lambda }^{-2}+0.00013{\lambda }^{-4}\right)$$

(4)

Where:

λ is the band specific mid-wavelength value (nm).

t_oz↑ and t_oz↓ were calculated using [34] (5,6):

$$t_{oz}\uparrow =\mathit{exp}\left(-{\tau }_{oz}\right)$$

(5)

$$t_{oz}\downarrow =\mathit{exp}\left({\displaystyle \frac{{-\tau }_{oz}}{cos{\theta }_s}}\right)$$

(6)

Where:

t_oz is the constant ozone optical thickness [35].

Finally, the Rayleigh-corrected radiance; ${\widehat{L}}_{\lambda }$was obtained by subtracting the Rayleigh scattering radiance (Lr) from the TOA radiance (L_λ) (7):

$${\widehat{L}}_{\lambda =}{L}_{\lambda }-L_r$$

(7)

Step 3. Converting radiance to reflectance:

TOA reflectance for each band ($\rho \lambda$), was calculated using (8):

$$\rho \lambda ={\displaystyle \frac{\pi \times {\widehat{L}}_{\lambda }\times d^2}{{ESUN}_{\lambda }\times {cos\theta }_s}}$$

(8)

Where:

ρλ is TOA reflectance for wavelength range or band λ, and

d is the Earth-to-sun distance in astronomical units [27].

2.3.4. Extracting Matchups

We extracted remote sensing reflectance (R_rs) for each band (B, G, R, and NIR; Section 2.3.3) and paired it with the corresponding sampling points (Section 2.3.1) in GEE (See Supplementary Material). To increase the chance of matching satellite observations with sampled Chl-a, we considered Landsat images within ±3 days of the sampling dates. Larger temporal windows were not used due to the highly dynamic nature of blooms in LW [36]. Satellite data were extracted from corresponding pixels to sampling point locations (single pixels; spatial windows of 1 × 1). To mitigate potential errors in digitization [37], we also extracted R_rs using the median of pixels corresponding to sampling locations and their eight (spatial window of 3 × 3) and 24 neighboring pixels (spatial window of 5 × 5). Based on these criteria, a total of 42 samples (18 in the NB and 24 in the SB) from 28 stations (13 in the NB and 15 in the SB) were identified (Figure 1, Table 1). We found 9 Landsat images for the NB, 11 for the SB, and 20 for the entire LW matching the time and location of the sampling points, providing 23 matchup points in the NB, 30 in the SB, and 57 for the entire LW (Table 1). We created 36 sets of matchups with different temporal windows (0 to ±3 days), different spatial windows (1 × 1, 3 × 3, and 5 × 5), and in each LW basin (NB-specific and SB-specific) as well as for the entire lake (LW-specific).

2.3.5. Model Calibration

For each set of matchups (Section 2.3.4), we tested the performance of 44 models, including all possible bands and band ratios, as well as various commonly used multiband combinations (Supplementary Material, Table S1) to predict Chl-a in LW using regression analyses. Each model was fitted using the matchup natural log-transformed in situ Chl-a values against the partially corrected R_rs (with and without natural log-transformation) using simple linear regression equations (Table S1). We then identified the Best Performing Model (BPM) for each set of matchups as the model that exhibited the highest coefficient of determination (R²) and was statistically significant (p < 0.05). Model calibration on the BPMs was performed using 5-fold cross-validation, chosen to balance bias and variance given the number of matchups available [38,39]. During model calibration, for each fold of the 5-fold cross-validation, the models were fitted to the training dataset (natural log-transformed in situ Chl-a values against partially corrected R_rs) using simple linear regression equations. We then estimated the average predictive performance of each model during the calibration stage using the average R² of the five cross-validation folds.

2.3.6. Model Validation

To assess the predictive capacity of the BPMs identified during the initial calibration phase for each set of matchups (Section 2.3.5), we validated each model by applying the slope and intercept from each calibration to its associated test matchup set in the cross-validation fold (k = 5). The predictive performance of the BPMs was measured the predictive performance of the BPMs using Root Mean Square Error (RMSE; µgL^-1) for each cross-validation fold as follows (9):

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\sum _{i=1}^n\left({\displaystyle \frac{\left({\left(y_o-y_p\right)}^2\right)}{n}}\right)}$$

(9)

Where:

$y_o$ represents the observed Chl-a values, and

$y_p$ represents the predicted Chl-a values.

We normalized the RMSE (NRMSE) for each cross-validation fold as follows (10):

$$\mathrm{N}\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}={\displaystyle \frac{\sqrt{\sum _{i=1}^n\left(\frac{\left({\left(y_o-y_p\right)}^2\right)}{n}\right)}}{\overline{yo}}}$$

(10)

Where:

$\overline{yo}$is average observed Chl-a.

We calculated the Root Mean Squared Logarithmic Error (RMSLE; µgL^-1) as follows (11):

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{L}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{n}}\sum _{i=1}^n\left(\mathit{log}\left[{1+y}_o\right]-log{[1+y}_p\right]{)}^2}$$

(11)

We also calculated the Mean Absolute Error (MAE; µgL^-1) and Mean Absolute Percentage Error (MAPE; %) for the five cross-validation folds to measure the predictive performance of the models as follows (12,13):

$$\mathrm{M}\mathrm{A}\mathrm{E}={\displaystyle \frac{\sum _{i=1}^n\left|y_o-y_p\right|}{n}}$$

(12)

$$\mathrm{M}\mathrm{A}\mathrm{P}\mathrm{E}=100\times \mathrm{m}\mathrm{e}\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{o}\mathrm{f}\left[{\displaystyle \frac{\left|y_o-y_p\right|}{y_p}}\right]for1=1,\dots ,n.$$

(13)

Bias was measured as follows (14):

$$\mathrm{B}\mathrm{i}\mathrm{a}\mathrm{s}={\displaystyle \frac{\sum _{i=1}^n\left(y_o-y_p\right)}{n}}$$

(14)

Finally, we compared the performances of the models based on the average errors for all the cross-validation folds [40].

2.4. Application of the BPMs to Landsat OLI

The Landsat TM and ETM+ were launched in 1984 and 1999 (respectively), and the Landsat Operational Land Imager (OLI) was launched in 2013. In our study, all matchups since 2013 were discarded because the in situ Chl-a data in those years were taken from depths below 0.1 m, which did not meet our matchup selection criteria (Section 2.3.1). To test the applicability of our BPMs for predicting Chl-a using OLI, we applied the models to concurrent OLI and ETM+ samples (with 1 to 8 days of difference between the two sensors overpassing times) from 19 stations distributed throughout LW between 2013 to 2017 (Figure 1).

For each OLI matchup, Level 1 Collection 2 OLI data (https://www.usgs.gov/landsat-missions/landsat-collection-2) for the bands used in the BPMs (i.e., B (450–515 nm), G (595–600 nm), and Red (630–680 nm)) were pre-processed in GEE using the same approach as for Landsat TM and ETM+ (Section 2.3.3). After applying the BPMs to ETM+ and their concurrent OLI data, we compared the modeled Chl-a derived from each sensor (a total of 26 matching points) using a simple linear regression model.

2.5. Comparing Chl-a Predictions Using Basin-Scale vs. Lake-Scale BPMs

To compare Chl-a predictions using basin-specific versus lake-specific BPMs, we applied the BPMs from each approach to all Landsat images during the peak Chl-a period (July–October) from 1984 to 2023 in GEE. We then compared the ability of the predictions from the two approaches to capture spatial and temporal changes in peak Chl-a in LW. All Landsat images were pre-processed to partially corrected R_rs in GEE (section 2.3.3) prior to model application.

To study spatial changes in peak Chl-a in LW, we calculated the median of peak Chl-a for each pixel within each basin for all years from 1984 to 2023. To study temporal changes in peak Chl-a in LW, we calculated the median peak Chl-a for all pixels within each basin for each year from 1984 to 2003.

3. Results

3.1. Spatial Heterogeneity in LW Turbidity

Our results revealed significant differences in turbidity levels among different regions of LW during the growing season (May-October). Specifically, turbidity in the NB was consistently lower than in the SB and the Narrows (p < 0.05), whereas no significant difference was observed between the SB and the Narrows (p ≥ 0.05; Table 2).

3.2. Best Performng Models

The performance of basin-specific models in predicting Chl-a was superior to that of lake-specific models (Figure 3 and Table S1). The NB-specific and SB-specific BPMs achieved R² values of 0.85 (n = 10) and 0.76 (n = 7), respectively (p < 0.05; Figure 3 and Tables S2, S3). In contrast, the lake-specific BPM showed a lower R² of 0.38 (p < 0.05; n = 17; Figure 3 and Table S4).

The NB-specific BPM (R² = 0.85) used matchups from ±1 day with a 1 × 1 spatial window (Figure 3 and Table S2). Increasing the spatial window to 3 × 3 and 5 × 5 decreased the R² to 0.70 and 0.60, respectively, while increasing the temporal window to ±2 and ±3 days decreased the R² to 0.75 and 0.57, respectively. Despite these changes, increasing the temporal window to ±2 days still provided a high R² of 0.75 with reduced errors (NRMSE decreased from 1.98 to 0.73; bias reduced from 9.39 to 2.04).

The SB-specific BPM (R² = 0.76) used matchups from ±1 day with a 1 × 1 spatial window (Figure 3 and Table S3). Increasing the spatial window to 3 x 3 and 5 x 5 decreased the R² to 0.74 and 0.63, respectively, while increasing the temporal window to ±2 and ±3 days decreased the R² to 0.61 and 0.32, respectively. The temporal window of ±2 days provided 17 matchups with an R² of 0.61 and reduced bias (bias decreased from 0.25 to 0.23).

Like the basin-specific models, the LW-specific model (R² = 0.38, based on matchups from ±1 day with a 1 × 1 spatial window) exhibited decreasing R2 values with larger spatial and temporal windows (Figure 3 and Table S4).

3.3. Best Chl-a Prediction Models

From the BPMs in section 3.2, the basin-specific models (NB and SB) using a 1 × 1 spatial window and a temporal window of ±2 days were selected for their balanced performance (high R² and low errors) and increased confidence due to a larger number of matchups used.

In the NB, a 2-band model based on the R_rs of the band ratio of G to B showed the highest performance (R² = 0.74) (Table 3 and Figure 4a) (15):

$$\ln \mathrm{C}\mathrm{h}\mathrm{l}-a=8.153\times \left[G/B\right]-5.622$$

(15)

In the SB, a 2-band model based on the R_rs of the band ratio of R to G performed best (R² = 0.62) (Table 3 and Figure 4b) (16):

$$\ln \mathrm{C}\mathrm{h}\mathrm{l}-a=4.313\times \left[R/G\right]-2.603$$

(16)

In the entire LW, a 2-band model based on the R_rs of the band ratio of G to B achieved the highest performance (R² = 0.38) (Table 3 and Figure 4c) (17):

$$\ln \mathrm{C}\mathrm{h}\mathrm{l}-a=5.836\times \left[G/B\right]-4.005$$

(17)

Validation tests using 5-fold cross-validation further supported the superiority of basin-specific models over the LW-specific model (Table 3).

3.4. Application of the Models to Landsat OLI

Cross-comparison between ETM+ and OLI sensors indicated a strong linear relationship between Chl-a predictions using OLI (dependent variable) and those using ETM+ (independent variable) (R² = 0.87, p < 0.05) (Figure 5). Notably, the correlation was stronger in the NB (R² = 0.83, p < 0.05) compared to the SB (R² = 0.75, p < 0.05) (Figure 5).

3.5. Comparing Chl-a Predictions Using the Basin-Scale vs Lake-Scale BPMs

Figure 6 illustrates the median (1984-2023) of natural log-transformed peak Chl-a (µgL^-1) in the NB and SB using basin-specific (Figure 6a) and lake-specific BPMs (Figure 6b). Basin-specific models predicted higher Chl-a concentrations in the NB compared to the SB (Figure 6a), while lake-specific models showed the opposite pattern (Figure 6b)

Figure 7 illustrates the annual time series of peak Chl-a in the NB and SB derived from basin-specific and lake-specific BPMs. In the NB, lake-specific models consistently predicted lower Chl-a concentrations compared to basin-specific models (Figure 7). Conversely, in the SB, lake-specific models predicted higher Chl-a concentrations compared to basin-specific models (Figure 7). Temporal variations in peak Chl-a from 2002 to 2019 also exhibited distinct patterns between basin-specific and lake-specific models, reflecting different spatial dynamics within LW (Figure 7).

4. Discussion

Remote sensing is a crucial tool for the rapid and consistent monitoring of spatial and temporal changes in surface water Chl-a. Limnologists often use global satellite-derived level-2 or level-3 Chl-a products designed for optically complex waters [36,41,42] without considering variations in water optical properties among lakes [40]. Chl-a prediction using remote sensing reflectance can be improved by developing models that account for potential within-lake variance in optical properties. Such models can be delineated based on optical properties or, as in our study, by morphometric characteristics that distinguish lake basins with unique optical properties. Our study demonstrates that by considering within-basin differences in optical properties and dividing the lake into optically distinct parts, Chl-a prediction can be improved.

4.1. Basin-Specific Chl-a Prediction Models

The physical and biogeochemical differences between LW’s NB and SB [19,36] create different optical properties, favoring basin-specific models over lake-specific models in both basins. BPMs in both basins used band ratios. Compared to single bands, band ratios enhance Chl-a estimations by reducing the atmospheric, irradiance, and air-water surface effects on reflectance [43]. Band ratios also minimize the impact of non-phytoplankton OACs found in inland waters, which can alter the spectral profile observed by the sensor [44,45]. This is significant because the varying concentrations of these constituents can explain the lack of correlation between observed Chl-a and single bands in lakes [43], as observed in our study (Supplementary Material).

Many models based on band ratios have been proposed for estimating Chl-a in inland, coastal, and open ocean waters [11,46]. These models leverage the reflection and absorption patterns of sunlight by Chl-a. With increasing Chl-a, the relative energy between 520 and 600 nm (Landsat G band) increases, while the energy between 450 and 520 nm (Landsat B band) decreases [47]. In our study, while BPMs in both NB and SB included the G band, the second band in the ratio differed between the two basins.

In the NB, where the water is deeper and clearer [19], Chl-a reaches higher concentrations during bloom events [19,35]. In this basin, the B to G band ratio provided the highest correlation with observed Chl-a (R² > 0.80). The B to G band ratio is widely used in estimating Chl-a in lakes with low concentrations of inorganic suspended matter or in eutrophic and hypertrophic waters where phytoplankton are the predominant constituents and non-phytoplankton constituents have a smaller effect on reflectance [48,49,50].

In the SB, where the water is shallow, well-mixed, and highly turbid [19,36], the R to G band ratio provided the highest correlation with observed Chl-a (R² > 0.70). The SB creates optically complex water conditions, where non-phytoplankton OACs interfere with spectral bands – particularly the B band. The concentrations of these constituents do not always correlate with Chl-a [51]. This interference complicates Chl-a estimation in optically complex turbid waters using the B band [44,52,53]. To account for the overlaps of the absorption of non-phytoplankton particles, the R to G band ratio is commonly used for estimating Chl-a in highly turbid waters [43,54,55,56]. Many Chl-a prediction models developed for turbid, optically complex waters are also based on the red and near-infrared regions of the electromagnetic spectrum [43,57]. Like the red region, the near-infrared region (near 700 nm wavelength) shows less sensitivity to non-phytoplankton particles' interference compared to the blue region [58]. However, near-infrared-based models are successful with satellites having a band centered at the red- near-infrared edge, like MERIS [58,59]. Multispectral sensors like Landsat have broad NIR bands (770–900 nm in TM and ETM+ and 845–885 nm in OLI), which have been shown to be unsuitable for estimating Chl-a [43]. In our study, we could not improve the relationship between observed and modeled Chl-a using NIR in either basin-specific or lake-specific models.

4.2. Application of the Models to Landsat OLI

Our matchup selection criteria excluded in situ matchups concurrent with the period of Landsat OLI operation (2013 to present). OLI differs from TM/ETM+ in image quality and the number and width of spectral bands. Despite these differences, the similar spatial and spectral resolutions allow for intercalibration between these sensors [60]. Our cross-comparison analysis showed that the models developed in this study are applicable for predicting Chl-a using OLI. Other studies have also successfully applied models developed based on TM or ETM+ to OLI images or vice versa [3,61].

The matchup selection criteria in our study left out in situ matchups concurrent with the period of Landsat OLI operation (2013 to present). OLI is different from TM/ETM+ for image quality and number and width of spectral bands. Despite the differences between TM/ETM+ and OLI, the similar spatial and spectral resolutions allow intercalibration between these sensors [60]. Our cross-comparison analysis showed the models developed in this study are applicable to predict Chl-a using OLI. Similarly, other studies successfully applied models developed based on TM or ETM+ to OLI images or vis versa [e.g., 3,61].

4.3. Future Applications

Overall, the high correlation between modeled and observed Chl-a obtained using the basin-specific models suggests these models can reliably capture relative changes in Chl-a in LW. Therefore, the models can be used for predicting Chl-a in LW with Landsat data to study spatial and temporal variabilities of Chl-a [61,62]. However, the lack of in situ Chl-a data with concentrations above 10 µgL-1 in the SB due to our matchup selection criteria may affect the accuracy of absolute Chl-a predictions in this basin [61]. More in situ Chl-a data can also enable Chl-a predictions using an optical water types (OWTs) approach, which requires large in situ Chl-a data for training. Using this approach, lakes are classified into OWTs before prediction [40,63,64], with unique Chl-a prediction models developed for each water type. Our matchup selection criteria also excluded any samples from the Narrows. Like the SB, the Narrows is optically complex, where non-phytoplankton particles are primary determinants of water transparency [19]. However, applying the SB-specific model to the Narrows requires in situ Chl-a data to validate these models' performance in this part of the lake.

Applying models to Landsat imagery (available since 1984) in GEE (see Supplementary Material for the code) can provide the longest available Chl-a time series for most water bodies. The length of this time series (40 years) allows for the study of long-term climate trends or oscillations' impacts on changes in Chl-a. At least 30-40 years of continuous and consistent data on water temperature (a critical indicator of climate change) and Chl-a is required to separate climate change signals from natural variability [64,65,66,67]. This information will also provide ecosystem managers and scientists with a better understanding of the sensitivity of different regions of their study area to environmental perturbations. This can lead to more accurate predictions of future lake ecosystem responses to increasing anthropogenic pressures in the face of a changing climate. Ecosystem managers can use this knowledge to design contextual and more efficient management strategies towards eutrophication and phytoplankton blooms in lakes.

5. Conclusion

In this study, we predicted relative concentrations of Chl-a in a large lake using simple empirical equations applied to Landsat data. Our results highlight the importance of considering optical differences within lakes when developing Chl-a prediction models. By separating a lake into basins with distinct water optical properties, we significantly improved the relationships between in situ and modeled Chl-a. These models can be easily applied in cloud-based platforms like Google Earth Engine (GEE) for computationally intensive modeling of Chl-a on broad spatial and temporal scales. The user-friendly GEE code requiring minimal adjustments of the times and locations of in situ Chl-a samples that was developed for this study is freely available, facilitating the development of models for optically distinct lakes world-wide.