Spatiotemporal Variations in Biophysical Water Quality Parameters: An Integrated In-Situ and Remote Sensing Analysis for an Urban Lake in Chile

1. Introduction

Lakes in urban areas contribute to the quality of life by providing ecological and social ecosystem services, defined as benefits that people receive from the environment. These lakes represent intricate systems that ensure the continued provision of a range of ecosystem services, especially in urban lake environments, demands consistent integrated actions within the surrounding watershed or directly within the lake itself [1,2]. Ensuring the quality of surface and groundwater is an essential duty carried out by organizations responsible for water supply and management [3]. It is crucial to implement appropriate conservation and management measures to protect and improve lake water quality, thus safeguarding their long-term ecological and socioeconomic value [4] by analyzing the trophic state of a lake. The trophic state of a lake indicates its level of biological productivity and water quality, which is determined by the availability of nutrients. Lakes can be classified into different trophic states, such as oligotrophic (low productivity), mesotrophic (moderate productivity) or eutrophic (high productivity) [5]. This state significantly influences the diversity and abundance of aquatic species, as well as the presence of algae and other organisms.

Usually, to monitor water quality a traditional point-based providing high quality data is used. However, these data usually represent isolated sampling points and do not fully reflect the comprehensive spatial and temporal variations within the aquatic ecosystem [6]. To address challenges related to spatial and temporal coverage, remote sensing has been recognized as an optimal solution, leveraging advancements in sensors and methodologies. This has led to extensive utilization of satellite technology in monitoring changes in freshwater bodies [7]. The usage of satellite remote sensing systems enables the observation and study of water quality on a regional and global level [8]. Thus, remote sensing offers a comprehensive and immediate perspective of aquatic, atmospheric, and terrestrial systems, and significant advantages over traditional monitoring methods, including large-scale coverage, real-time monitoring capabilities, and cost-effectiveness. It serves as an optimal tool for continuous monitoring of the water quality in surface bodies [9,10,11]. When combined with in-situ measurements, it can deliver unparalleled, cost-efficient, and quantitative assessments of change. It provides spatial resolutions reaching the scale of meters and temporal frequency surpassing daily intervals [12]. One of the most important sources of remote sensing data is the Landsat archive. Landsat-8 (L-8) Operational Land Imager (OLI) for this study is based on its accessibility, suitable temporal frequency, spatial resolution, and the methods for utilizing the instrument to determine water quality data [13,14]. Additionally, this sensor has a high radiometric resolution, crucial for recovering constituents, since water bodies act as dark objects and absorb a larger fraction of the downward irradiance, resulting in a low signal-to-noise ratio (SNR) [15]. Satellite images acquired by Landsat satellites facilitate monthly monitoring of specific regions within inland water bodies.

When examining lake water quality through satellite imagery, chlorophyll-a (Chl-a - in mg·m^-3) and turbidity (Nephelometric Turbidity Units - NTU) stand out as two essential Apparent Optical Parameters (AOPs) to address this challenge [16,17]. Chl-a is commonly used to assess the trophic status of coastal and inland waters due to its sensitivity and rapid response to various environmental factors. Chl-a is the green pigment found in aquatic plants, such as algae and cyanobacteria [18]. An increase in Chl-a concentration can indicate the proliferation of these organisms, often associated with eutrophication, a phenomenon where excess nutrients in the water (such as nitrogen and phosphorus) promote the excessive growth of algae [19], which can have negative effects on water quality and the health of the aquatic ecosystem. Additionally, elevated turbidity levels (NTU) can result in an opaque, hazy, or cloudy appearance of the water, which is closely related to suspended sediment concentrations [10,20].

The Laguna Grande is in the commune of San Pedro de la Paz in the Central-Southern Chile. Its watershed has undergone significant modifications due to the introduction of forestry plantations (Pinus radiata and Eucalyptus globulus) and urban expansion [21]. Over time, a multitude of external factors and natural causes have contributed to the deterioration of the lagoon's water quality. Primarily, this decline can be attributed to the occurrence of eutrophication processes. Eutrophication is a natural phenomenon that typically unfolds over thousands of years. However, human activities can expedite this process, leading to the emergence of what is commonly referred to as artificial or cultural eutrophication [22]. To comprehensively examine and analyze the water quality, a few decades ago, the "Dirección General de Aguas – DGA" (Water Agency of Chile) established three sampling stations within the lagoon (Sector Desagüe, Sector Centro, and Sector Sur). These stations facilitated the collection of traditional measurements, e.g., Chlorophyll-a, turbidity, pH, CDOM, to assess the lagoon's water quality and its trophic condition. These three sampling stations operated from 1988 to 2013, but unfortunately these stations are no longer operational. Since 2014, sampling is performed only at the Sector Centro 2 station, but the available information only includes three sampling missions by year between 2014 and 2016. Assuming four measurements per year, the percentage of missing data at stations Centro, Desagüe, Sur and Centro 2 between 1988 and 2016 is 28% [23]. These factors serve as a focal point to address the main problem that this study aims to solve.

Most of the studies carried out in the "Laguna Grande San Pedro de La Paz" by remote sensing are usually done on a fixed time scale. For example, Lilo-Saavedra et al. [24] have shown the accuracy of estimating representative water quality features in the lagoon by utilizing fused satellite images. Quintana-Sotomayor et al. [25] published a study implementing an object-oriented classification methodology, for multispectral images, to quantify turbidity levels and generate a turbidity map of the lagoon. Rojas [26] focused its study on assessing the effects of land use changes in the watershed of Laguna Grande on sedimentation rates. This was achieved through a multitemporal analysis using satellite imagery. Recently, Pedreros-Guarda et al. [23], implemented six automated methods for retrieving the surface water temperature of the lagoon using satellite thermal imagery.

Thus, our main objective is to develop spectral inversion models based on in-situ data versus reflectance signal to estimate two of the most important Water Quality Parameters (WQPs) for eutrophication dynamics monitoring, i.e., Chl-a concentration and turbidity. Furthermore, these models will play a key role in filling the spatial and temporal gaps in the historical data records stored in the DGA database. The empirical retrieval models would facilitate the reconstruction of variable behavior during specific time periods when field data were not collected, thus increasing the robustness of the dataset. Our study will encompass the following key aspects: i) Implementation of atmospheric correction methods: different atmospheric correction methods will be employed to correct atmospheric distortions and the sun-glint effect in satellite data acquired from Landsat-8 OLI sensor; ii) Development of an empirical retrieval model: we will construct a spectral inversion model tailored to accurately quantify water quality parameters, specifically Chl-a and turbidity, including an assessment of statistical effectiveness; and iii) Assessment of the spatial-temporal variability: this analysis will allow to understand the eutrophication processes during different seasonal conditions and multiple years.

2. Materials and Methods

2.1. Study Area

The Laguna Grande de San Pedro de la Paz (Figure 1) is a shallow urban lake placed in the northwestern part of the Nahuelbuta mountain range (36° 51' S, 73° 06' W) at 220 m. a.s.l., south of the Biobio River, approximately 5 km from the city of Concepción. It has a perimeter of 9800 m, a maximum length of 2675 m, and a maximum width of 1375 m. It covers an area of over 155 ha and reaches a maximum depth of 13.5 meters, with an average depth of 8.1 meters [21]. On its eastern side, Laguna Grande is bordered by a low-altitude mountain range that does not exceed 360 m.a.s.l. This area is characterized by forest plantations, primarily consisting of pine and eucalyptus trees (Pinus radiata and Eucalyptus globulus, respectively). Moreover, in the hills to the north, there has been a major urban and real estate development project in recent years. Laguna Grande drains through the Los Batros estuary, receiving water from natural runoff as well as drainage from the surrounding residential areas [26]. The current trophic condition (using the Trophic State Index – TSI) of Laguna Grande de San Pedro was categorized as mesotrophic (30 < TSI < 60), based on the Carlson [5] TSI classification and the data obtained by Urrutia et al. [27].

2.2. Field Data of Water Parameters

For this project, three field surveys were conducted to collect surface water samples: i) the first one in the spring season (October 24^th, 2022); ii) the second one in the summer season (January 6^th, 2023); and iii) the last one near to the fall season (March 1^st, 2023). The selected dates were previously chosen on the current schedule of L-8 satellite overflights over the lagoon. This selection aimed to align our field trips with the satellite's image acquisition dates, minimizing the time gap between image capture and sample collection.

To determine the sampling station locations, a grid was employed. Forty (40) sampling points were distributed on the lake (Figure 1). Due to logistical constraints, we collected 40 samples in the summer field trip, 25 in the spring trip, and 20 in the fall trip, resulting in a total of 85 samples. These samples were collected using 1-liter glass bottles, typically positioned approximately 20-30 cm below the water's surface. The samples were analyzed in the EULA-Chile Center laboratory and the results are summarized in the Table 1. For Chl-a analysis, the samples were subjected to a filtration, extraction, and centrifugation, following the protocol of EPA method 445.0 [28]. Subsequently, fluorometry was used to measure Chl-a content, and unit transformations and calculations were performed using the corresponding calibration curve for the Turner Designs model 10-AU-005-CE equipment. With this process, Chl-a concentrations were obtained in milligrams per cubic meter (mg·m^-3). For the turbidity determination, a specialized turbidimeter (2100Q Portable Turbidimeter Hach) was used to obtain the results in nephelometric turbidity units (NTU).

In the last fieldtrip other important water quality parameters were measured, such as: Temperature (°C), dissolved oxygen (mg·l^-1), oxygen saturation (%), conductivity (mhos·cm^-1) and pH, which allowed us to know the trophic state of the lake of the measurements at-satellite overpass (Table 2).

2.3. Field Measurement of Water Spectral Signatures (R_rs)

We performed reflectance (R_rs) measurements at 15 sampling stations using an ASD FieldSpec-4 spectroradiometer (ASD Inc., Malvern Panalytical). This instrument is specifically designed to operate within a wavelength ranging from 350 nm to 2500 nm.

The above-water remote sensing reflectance was calculated by the equation (1):

$${\mathrm{R}}_{\mathrm{rs}}= \frac{L_w \left(\lambda \right)}{E_d \left(\lambda \right)} \left(s{r}^{-1}\right)$$

(1)

Where, R_sr = above-water remote sensing reflectance (dimensionless); L_w(λ) = spectral radiance of the target (W/cm²·sr·µm); E_d(λ)= spectral radiance of the reference plate (W/cm²·sr·µm). A Spectralon reference plate was used in the study to simulate a Lambertian surface. Measurements were made between 10:00 and 13:00, a time of optimal conditions with low solar tilt, high radiant energy, and favorable weather. Radiometric measurements followed the methods established by Milton [29] and Mobley [30], performed under conditions of low wind (0 - 4 m·s^-1), clear skies and solar angles between 0 and 30°. To minimize the sun-glint effect, the spectroradiometer was positioned at an inclination of 40° with respect to the vertical and 135° from the sun. The sensor was placed 1 meter above the water surface. Ten repeated measurements were established for each station, verifying the results directly in the field. Some anomalous spectra caused by random solar flares may be present in the measured result. However, these were suppressed, and the remaining spectra were averaged (Figure 2).

ENVI-IDL software was used to construct a spectral library of surface water samples at different Chl-a concentration conditions. This spectral library of seven spectral profiles were used to evaluate the performance of the different atmospheric correction methods in the OLI L-8 sensor scenes.

2.4. Satellite Image Data and Atmospheric Correction Methods

2.4.1. Landsat-8 OLI Satellite Imagery

The Laguna Grande occupies a geographic position at the intersection of two orbits, namely Path/Row combinations 001/085 and 001/086, resulting in a region of overlap covering approximately 31 km. It should be noted that the satellite images used in this study were obtained from the official website of the United States Geological Survey (USGS), specifically from the Collection 2, Levels 1 and 2. These images can be accessed at the following link: https://earthexplorer.usgs.gov (accessed on June 15, 2023).

Three satellite images from the L-8 OLI sensor (Table 3) were used to develop the empirical models for retrieving bio-optical parameters, i.e., Chl-a and turbidity. The empirical models were developed by a combination of the signal reflectance and in-situ data. The time difference between in-situ data collection and satellite image acquisition does not exceed more than 3 days. No rainfall events occurred during the whole data collection. Seven additional L-8 OLI images were used to study the spatiotemporal variability of bio-optical parameters (Chl-a and turbidity) in the lake during different seasonal conditions between hydrological cycle for 2014 and 2022.

2.4.2. Atmospheric Correction Methods

In this research, we conducted a comparison of four different Atmospheric Correction (AC) methods designed to obtain surface water reflectance from L-8 OLI radiance images. Several lake stations were used to compare the performance of these four atmospheric correction models:

(1)

The first tested model was ACOLITE, which is a processor developed by RBINS (Royal Belgian Institute of Natural Sciences) used to apply the AC and satellite image preprocessing approach in applications related to inland and coastal waters. ACOLITE uses the approach known as dark spectrum adjustment to perform the atmospheric correction [31,32,33,34,35]. The main objective is to remove atmospheric influences, such as scattering and absorption, to improve the accuracy of remote sensing data for different applications. Sun-glint is the specular reflection of sunlight off the surface of water. It can introduce significant errors in remote sensing data, especially in coastal and oceanic regions. ACOLITE includes algorithms to detect and quantify sun-glint in satellite imagery. Once sun-glint is detected, ACOLITE applies correction algorithms to adjust pixel values to compensate the overestimation of reflectance caused by the effect of the solar reflection. In this study the v20210114.0 version of the processor was used, and its default configuration file was modified to apply the reflection correction to the image, and also to obtain the results in GeoTIFF format.

(2)

The second model was iCOR, which is a tool developed by De Keukelaere et al. [36] that can process satellite data acquired over coastal, inland, or transitional waters and land. iCOR employs the Moderate-Resolution Atmospheric Radiance and Transmittance Model-5, known as MODTRAN5 [37], to perform radiative transfer calculations. In addition, it uses Look-Up Tables (LUTs) to speed up retrieval processes. An important aspect of iCOR is its ability to identify whether a pixel belongs to a water or land area, allowing it to apply a specific atmospheric correction accordingly [38]. The version iCOR used was v3.0 into the SNAP software.

(3)

The third model was Land Surface Reflectance Code (LaSRC) whose algorithm was developed by E. Vermote [39], National Aeronautics and Space Administration (NASA) Goddard Space Flight Center (GSFC) and was modified by the USGS Earth Resources Observation and Science (EROS) center. LaSRC generates top-of-atmosphere reflectance (TOA) and top-of-atmosphere brightness temperature (BT) using the calibration parameters provided in the metadata. Then, atmospheric correction routines are applied to the L-8 TOA reflectance data using additional information such as water vapor, ozone, and aerosol optical thickness (AOT) obtained from the Moderate Resolution Imaging Spectroradiometer (MODIS). In addition, the digital elevation model derived from the Earth Topography Five Minute Grid (ETOPO5) is used to generate surface reflectance [39]. This product can be downloaded from the following the USGS website (https://www.usgs.gov/landsat-missions/landsat-surface-reflectance (Last accessed on june 15, 2023). A scaling equation (2) was applied to this final product to normalize the reflectance values between 0 and 1, allowing comparison with the other methods.

Rrs(λ) = 0.0000275 x Pixel Value - 0.2

(2)

(1)

The last method was C2RCC (Case 2 Regional CoastColour) atmospheric correction processor, which has a deep learning approach using set of neural networks trained and linked to simulated reflectance data in water bodies, and radiances in the upper atmosphere. Its main outputs are associated with the inherent optical properties (IOPs) of water, i.e., those that depend exclusively on the absorption and scattering of its constituents [40]. This method considers three sets of neural networks for the calculation of reflectance depending on the research objective: C2RCC-Nets (Standard neural network suggested for use in eutrophic or mesotrophic water bodies), C2X-NETS (Specialized neural networks the water bodies with high concentrations of suspended matter and chlorophyll concentration) and C2X-COMPLEX-Nets (suggested mainly for use in inland waters) [41]. C2RCC can be used as a complement in the SNAP software and allows the calculation of reflectance in Sentinel 3 OLCI, Sentinel 2 MSI, Landsat-8/9, MODIS and MERIS satellite images (e.g., [42,43]).

2.5. Methodology for Water Quality Modeling

To determine the water quality parameters (Chl-a and turbidity) from in-situ and satellite data (OLI L-8 images), we implemented a complete workflow (Figure 3). This methodology includes the following main steps: i) in-situ acquisition of spectral signatures and water sample collections, accompanied by the selection of corresponding OLI L-8 scenes; ii) determination of Chl-a (mg·m^-3) and turbidity (NTU) concentrations, and evaluation of the atmospheric correction models; iii) identification of the pixels with atmospheric distortion (high aerosol levels) to avoid to introduce external error to the spectral retrieval model; iv) determination of the spectral bands (with indices and reflectance band ratios) showing the best fit regressions with the in-situ values, in which ACOLITE and C2RCC were selected as the best atmospheric correction models to estimate Chl-a and turbidity values, respectively; v) statistical assessment by employing the bootstrapping methodology to enhance the spectral retrieval model pertaining to each water quality parameter. Furthermore, we utilized a set of five distinct statistical metrics to determine the precision and accuracy of the modeling process; and vi) generation of the spatiotemporal images for the water quality parameters at different dates and seasonal conditions.

2.6. Selection of Spectral Indices

Spectral indices are used to measure the distinguishing characteristics of vegetation, water bodies and various land cover types [44]. Low to medium resolution satellites have been a widely used to apply eutrophic inland water classification algorithms and to retrieve Chl-a concentration [45]. In terms of optics, Chl-a exhibits two points of highest reflectance, located in the green and near-infrared wavelength ranges at 550 nm and 700 nm, respectively [3]. A total of three spectral indices (Table 5) have been previously used to detect the presence of Chl-a in water by remote sensing.

The Normalized Difference Vegetation Index (NDVI) is calculated by the formula NDVI = (B5 - B4) / (B5 + B4), where “B5 - NIR” represents the near-infrared spectral reflectance, predominantly influenced by plant canopy reflection, and “B4 - Red” means the reflectance within the red part of the electromagnetic spectrum, where chlorophyll absorption reaches its maximum [46]. The Green Normalized Difference Vegetation Index (GNDVI) [47] is based on the concept that an index designed to estimate chlorophyll content should remain constant regardless of pigments other than chlorophyll and should not be influenced by external factors, such as background or atmospheric conditions [48]. It is calculated by the formula GNDVI = (B5 - B3) / (B5 + B3), according to the band configuration of the OLI L-8 sensor. The Green Chlorophyll Index (GCI) [49,50] is used to assess the chlorophyll content of leaves. By including the near-infrared (NIR) and green wavelengths, this provides an improved assessment of Chl-a content [44]. It is calculated by the formula GCI = (B5 / B3) -1.

We based our evaluation of the turbidity parameter on the results of previous work [11,44,51,52]. These studies have demonstrated that reflectance values in the red and NIR bands maintain a strong correlation with suspended sediment concentration and turbidity levels. As a result, these spectral bands have become crucial indicators for assessing turbidity in water bodies [52]. With this in mind, we use empirical regression analysis techniques to establish a relationship between field measurements and satellite data. For turbidity (NTU) estimation we use different classical combinations for this parameter (Table 4), including: NIR/Red; NIR; Blue/Green; Red+NIR; together with the calculation of the Normalized Difference Turbidity Index (NDTI) [51].

2.7. Statistical Assessment

To assess the effectiveness of the selected bio-optical models, we implemented five statistical and error metrics: mean absolute error (MAE; Eq. 3), mean square error (MSE; Eq. 4), root mean square error (RMSE; Eq. 5), mean absolute percent error (MAPE; Eq. 6), and the commonly employed coefficient of determination (R²; Eq. 7). Additionally, to estimate the uncertainty in our correlation analysis, we employed bootstrapping methodology in the assessment of Chl-a and turbidity obtained by spectral indices and in-situ data. Bootstrapping involves creating multiple data sets by sampling with replacement from our original dataset, and then calculating the correlation in each of these sets. We performed 1000 iterations of the Bootstrap analysis on our data, using 80% of the dataset in each iteration to compute correlation and the remaining 20% for validation.

$$MAE=\frac{{{\displaystyle \sum }}_{i=1}^N\left|X_i^{estimated}-X_i^{measured}\right|}{N}$$

(3)

$$MSE=\frac{1}{n}\times {\displaystyle \sum }_{i=1}^n{\left(X_i^{estimated}-X_i^{measured}\right)}^2$$

(4)

$$RMSE=\sqrt{\frac{1}{n}\times {\displaystyle \sum }_{i=1}^n{\left(X_i^{estimated}-X_i^{measured}\right)}^2}$$

(5)

$$MAPE=\frac{100}{n}\times {\displaystyle \sum }_{i=1}^n\left(\frac{X_i^{estimated}-X_i^{measured}}{X_i^{estimated}}\right)$$

(6)

$$R^2=1-\frac{{{\displaystyle \sum }}_{i=1}^n{\left(X_i^{measured}-X_i^{estimated}\right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left(X_i^{measured}-{\overline{X}}^{estimated}\right)}^2}$$

(7)

3. Results

3.1. Field Measurements of Water Parameters

In total, 85 water samples were collected under different seasonal conditions. For the first phase of the study, 40 measurements of Chl-a and turbidity were carried out, allowing the creation of concentration maps by interpolating these parameters. The purpose was to identify the areas of the lake with low, intermediate, and high concentration levels, thus highlighting those areas with greater variability in those values. This information allowed us to optimize subsequent field surveys, reducing the number of samples required to a range of 20-25 stations. As nearby stations exhibited similar values, stations with the highest variations of Chl-a and turbidity were utilized to build the spectral retrieval model. Thus, in the last two field trips, we focused on the stations and areas of the lake with the highest variability in Chl-a and turbidity, to improve the efficiency of our measurements in the lake. During these three field surveys in 2022-2023, a maximum Chl-a value of 8.7 mg·m^-3 was recorded in the near-fall season, while the lowest Chl-a values, 2.2 mg·m^-3, were observed during the summer season, during 2023. The mean Chl-a values recorded are as follows: 4.7 mg·m^-3 in the spring season; 3.4 mg·m^-3 in the summer season; and a value of 6.4 mg·m^-3 during the season close to fall. According to Carlson's [5] classification, the Laguna Grande maintains a mesotrophic state, i.e., moderate productivity in all the seasons that was evaluated. The maximum value of turbidity was recorded also in the near-fall season with a value of 20.2 NTU, while the lowest value was observed during the summer season with 0.6 NTU. Mean values were 2.4 NTU (spring season), 6.4 NTU (summer season) and 7.3 NTU (near to fall season). The average values do not exceed 5 NTU in the spring and summer seasons, so it can be considered a low turbidity value, typically associated to well-preserved clean waters. However, there is an increase in turbidity values near the fall season with values above 5 NTU, representing a moderate turbidity, which could be related to an increase of suspended sediment concentrations or biological activity.

3.2. Atmospheric Correction

3.2.1. Evaluation of Aerosol Levels in the Sampling Station Points

In-situ spectral signatures of water surface were used to compare the four atmospheric correction methods. To this purpose, we selected the L-8 OLI image in a date (acquired on 9^th January 2023) closest to the collection of water spectral library (6^th January 2023), with just three days apart. Then, the product LaSRC SR_QA_AEROSOL was implemented to determine which sampling stations are suitable to compare in-situ data with reflectance band combinations (Figure 4).

Data processing started with the location of the station points on the map. To extract the reflectance values from the scene, buffer zones were calculated around each station. These zones have different sizes, including 60-meters, 40-meters, and 20-meters in diameter, instead of a punctual signal, which improves the representativity of the reflectance. In the case of stations located near the shore, we apply a coherent buffer size of 40- and 20-meters to avoid contact with the lake shore. This adjustment was intended to mitigate spectral mixing between reflectance signals from water and soil materials, to avoid wrong interpretations in our model. The LaSRC SR_QA_AEROSOL product was used to select those pixels that are not affected by high aerosol values, ensuring better correlations between the in-situ data and the results of the band combinations (Table 5).

Table 5. Adaptation of Landsat 8-9 SR_QA_AEROSOL Values table.

Attribute	Pixel Value (DN)
Low-level aerosol	66, 68, 96, 100
Medium-level aerosol	130, 132, 160, 164
High-level aerosol	192, 194, 196, 224, 228

Table 5. Adaptation of Landsat 8-9 SR_QA_AEROSOL Values table.

Attribute	Pixel Value (DN)
Low-level aerosol	66, 68, 96, 100
Medium-level aerosol	130, 132, 160, 164
High-level aerosol	192, 194, 196, 224, 228

3.2.2. Evaluation of Atmospheric Correction Methods

Figure 5 shows the evaluation of different atmospheric correction methods, using station 15 as a “clean” reference with no aerosol contamination. We compared four atmospheric correction methods: ACOLITE, iCOR, LaSRC, and C2RCC. In our analysis, we excluded SWIR bands from the OLI sensor and focused on the bands coastal, blue, green, red, and NIR, which are the most sensitive wavelengths for Chl-a and turbidity estimation. The black line depicts the measured water surface reflectance at station 15, which was collected by using the ASD Field-Spec 4 spectroradiometer.

Although seven stations of data collection were used for this comparison, for illustrative purposes only, station 15 is shown in the results. The red segmented line, which corresponds to the LaSRC model's spectral signature, appears closer to the real reflectance (black line). However, in this model an overestimation is highlighted in the B3 (green) and B5 (NIR) data, unlike to other atmospheric correction models. This indicates that the LaSRC model does not take into account the sun-glint effect on the water surface.

Comparison of satellite and field results showed that ACOLITE (cyan segmented line) yielded results more aligned with glint-free field measurements. Although the values in all bands are slightly lower than the real reflectance signature, these maintain a similar shape to the black line, as well as the LaSRC method, but without saturation in B3 (Green) and B5 (NIR). It is important to note that these Green and Near-Infrared bands are essential for creating Chl-a retrieval models. It is due to the Chl-a absorbs light in the blue and red regions of the electromagnetic spectrum, while reflects light in the green and near-infrared regions [54]. This property allows us to estimate the concentration of Chl-a in water by examining the reflectance ratio between the green and near-infrared bands. iCOR (orange segmented line) and C2RCC (purple segmented line) showed lower values compared to the real reflectance, and even to the LaSRC and ACOLITE models in all bands. C2RCC showed a low value at B1 coastal with respect to iCOR, which could indicate an overcorrection at this wavelength using this model. In general, both models (iCOR and C2RCC) present an underestimation in reflectance values. Therefore, we have chosen ACOLITE as the optimal choice to advance our analysis in Chl-a and Turbidity recovery modeling.

3.3. Empirical Retrieval Models to Chl-a and TURBIDITY

3.3.1. Chl-a Estimation Analysis

Prior to analyzing the OLI satellite image's reflectance signal, we compared the in-situ Chl-a data with the 15 measurements of the water surface signal, acquired by using the ASD FieldSpec-4 spectroradiometer. Correlation coefficients R² were calculated to examine the spectral relationship between Chl-a concentration and the reflectance at the lagoon (Figure 6). The coefficient of determination (R²) varied over the wavelength range from 350 to 1000 nm. The obtained R² values were greater than 0.6 in the red region from 650 nm, which means the relationship between Chl-a concentration and the reflectance of these wavelengths are the best. The R² values shown a remarkable decrease around 925-950 nm (a region with wavelengths longer than the OLI NIR band), but those values quickly increase up to 0.66, even at wavelengths near 1000 nm.

The Chl-a has a unique absorption band that varies with its concentration, enabling the detection of specific spectral signatures through remote sensing. The data reveals that the absolute values at wavelengths of 655-880 nm and 990 nm exhibit the highest values of R². This suggests a stronger correlation between Chl-a concentration and the reflectance at these specific wavelengths. In Figure 7, the R² values obtained for the relationship between in-situ Chl-a concentration and the reflectance values at 655 nm and 880 nm are 0.6 and 0.656, respectively. These coefficients indicate an important correlation in these spectral bands, aligning closely with the red and near-infrared (NIR) regions of the electromagnetic spectrum.

The Pearson's coefficient values were used to assess the relationship between three vegetation indices and chlorophyll, i.e., NDVI, GNDVI, and GCI, and the Chl-a concentration values, using data from the green, red, and NIR regions. Figure 8 shows a heat map with the Pearson's coefficient values of these three spectral indices. Notably, these indices (NDVI, GNDVI, and GCI) shown a statistically significant correlation (p ≤ 0.001), with Pearson's coefficient values of -0.87 for NDVI and GNDVI, and -0.88 for GCI. Hence, GCI was selected as the best option to build the retrieval model. The sensitivity of Chl-a concentration to multiband reflectance combinations is clearly demonstrated in Figure 8, highlighting its higher responsivity compared to single-band reflectance.

3.3.2. Turbidity Estimation Analysis

Turbidity estimation using satellite imagery was based on spectral indices and combinations of bands related to suspended material in the water. Several spectral indices and band combinations were evaluated: NIR, Red, NIR/Red, Red+NIR, Blue/Green, and NDTI to determinate the best fit r with the turbidity values (Figure 9). The obtained values from the lower to the higher were: NIR (Pearson's coefficient = 0.098); NIR/Red (Pearson's coefficient = -0.34); NDTI (Pearson’s coefficient = 0.62); Blue/Green (Pearson's coefficient = -0.87); Red (Pearson's coefficient = 0.87); and Red+NIR (Pearson's coefficient = 0.88). The spectral indices with the best performance were Blue/Green, Red, and notably, Red+NIR, which showed statistically significant correlation (p ≤ 0.001), with a Pearson's coefficient value of 0.88.

3.4. Statistical Evaluation and Model Robustness

3.4.1. Statistical Evaluation and Robustness of the Chl-a Estimation Model

To assess the reliability of our correlation analysis between Chl-a concentration and the GCI index (best fit R²), we applied a bootstrapping technique. Through the 1000 iterations, the distribution of the parameters of the function was obtained to optimize the selection of them between the mean and the median (Figure 10). Consequently, equation 8 represents the final Chl-a recovery model.

$$Chla=1.97\times e^{\left(1.46\times GCI\right)}$$

(8)

Statistical and correlation metrics yielded robust results with a R² = 0.88 and MAE = 0.35%; MSE = 0.20; RMSE = 0.45% and MAPE = 7.56%, ensuring the stability of the observed relationship (Figure 11.). This bootstrapping approach helped to provide a robust and reliable estimate of the relationship between the variables.

To evaluate trends in estimation errors, we examined the model residuals using scatter plots against the in-situ collection data. No discernible pattern emerged, indicating the absence of systematic trends in residual variations associated with changes in Chl-a concentration; in other words, the errors appear to be randomly distributed (Figure 12). Residuals are uniformly scattered around zero, with all residual errors plotting approximately within the range from -1 to 1.

3.4.2. Statistical Evaluation and Robustness of the Turbidity Estimation Model

Similar to Chl-a estimation model, bootstrapping method was performed using 1000 iterations, to optimize the selection of function parameters between mean and median (Figure 13). As a result, Equation 9 now serves as the final turbidity retrieval model.

$$Turbidity=0.51\times e^{\left(130.75\times \left(Red+NIR\right)\right)}$$

(9)

Statistical and correlation metrics also yielded intermediate results with R² = 0.79, MAE = 1.15%; MSE = 2.11; RMSE = 1.45% and MAPE = 54.89% (Figure 14). The latest MAPE results are unsatisfactory as the error exceeds 10-20%. This suggests that the correlation between the Red+NIR index and in-situ turbidity data is not robust enough to be considered as a reliable estimation model.

Figure 15 shows a data clustering, with one cluster with very low values (< 4 NTU) and another more dispersed cluster with slightly higher values (> 7 NTU) from 7 to 10 NTU, suggesting a seasonal behavior in the data.

It is important to note that the turbidity levels in this lake are particularly lower than those in highly eutrophic lakes [55]. This low turbidity values may result in a relatively weak spectral relationship with reflectance bands.

3.5. Spatial and Temporal Variability

Statistical results indicate a strong relationship (R² = 0.88) for the Chl-a retrieval model using the CGI index, but the turbidity retrieval model did not accurately reproduce in-situ data. As a result, our focus was on Chl-a distribution to analyze data variability across the entire lake retrospectively. Figure 16 displays Chl-a concentration distribution maps obtained from our trained model (equation 8) for different seasons in two periods: 2022 (Figure 16a-d) and 2014 (Figure 16e-h).

The data for the Laguna Grande system consistently falls within the mesotrophic classification, with Chl-a concentration levels ranging from 3 to 10 mg·m^-3 (Figure 17). These values imply a moderate level of biological productivity, which has remained stable across various seasons and periods.

A detailed analysis of the 2014 and 2022 data reveals fluctuating Chl-a mean values. They start at 6.2 mg·m^-3 in winter, drop to 5.1 mg·m^-3 in spring, increase to 5.4 mg·m^-3 in summer, and peak at 7.1 mg·m^-3 in fall (Figure 17). This pattern is consistent in both the northern (Figure 17a) and southern zones (Figure 17b) of the lake. However, in the southern part, mean values are consistently lower for each season: 5.5 mg·m^-3, 4.6 mg·m^-3, 5.2 mg·m^-3, and 6.2 mg·m^-3, respectively (Figure 17b).

An examination for the variation in Chl-a levels for the year 2014 revealed distinct patterns. In the northern zone, Chl-a values ranged from 7.4 mg·m^-3 in winter, decreased to 6.0 mg·m^-3 in spring, reached a peak of 6.7 mg·m^-3 in summer, and then decreased to 6.2 mg·m^-3 in the fall season (Figure 17a). Meanwhile, the southern zone exhibited a similar trend, albeit with consistently lower mean values. Starting at 6.8 mg·m^-3, Chl-a concentrations decreased to 5.0 mg·m^-3, increased to 6.2 mg·m^-3, and then decreased again to 5.6 mg·m^-3 (Figure 17). This comparative analysis provides valuable insights into the seasonal Chl-a dynamics in both the northern and southern regions, shedding light on how these patterns evolved over time.

Chl-a values consistently exhibit higher levels in the northern zone compared to the southern zone, regardless of the year being considered (2022 versus 2014; Figure 17). This disparity can be attributed to a more pronounced anthropic impact in the northern zone, characterized by rapid urban sprawl over recent years, while the southern zone bordered by dense vegetation has remained relatively unchanged. Recent regulations implemented by local authorities to curtail activities contributing to lake pollution have further reinforced this trend. Notably, the Chl-a values for 2022 are lower than those observed in 2014, aligning with government restrictions on lake access, such as the use of fuel-powered boats and new construction around the lake shore.

Upon closer examination, when comparing the summer and fall 2014 data, an anticipated trend is evident: a decrease in Chl-a values from 6.2 mg·m^-3 to 5.6 mg·m^-3 (Figure 17). This aligns with the expected decline in temperature between the two seasons. However, a contrasting pattern emerges when analyzing the 2022 data (Figure 17), with Chl-a values in fall surpassing those in summer. This unexpected behavior can be attributed to the fact that recent autumnal temperatures resemble historical summer temperatures, possibly driven by regional climate variability. These findings highlight the OLI sensor potential in forecasting algal bloom distribution and trophic state in the lake.

4. Discussion and Conclusions

Our main findings in the development of spectral retrieval models for water quality parameters in Laguna Grande can be categorized into three key issues:

I.

Implementation of atmospheric correction methods

We performed the first water quality study in Laguna Grande by applying an innovative methodology, which involves the collection of in-situ data (water quality parameters) combined with spectral signatures at the water surface using ground-based radiometers. These measurements were then integrated to calibrate and validate the reflectance signal of the multispectral satellite images (OLI sensor) using different atmospheric correction methods. ACOLITE provided results that closely matched glint-free field measurements. Although the values in all bands were slightly lower than the actual reflectance signatures (surface water spectral profile), they retained a similar shape. When comparing the in-situ Chl-a and turbidity data with the reflectance signals corrected by the different atmospheric methods evaluated, ACOLITE showed the highest coefficients of determination, reaching R² values of 0.88 and 0.79 for Chl-a and turbidity, respectively. These results are in line with a recent study by Rodríguez-López et al. [9], which used ACOLITE as an atmospheric distortion mitigation method on L-8 OLI scenes to estimate Chl-a concentrations in the Llanquihue Lake (southern Chile). Other water quality research has found that ACOLITE effectively corrects for atmospheric distortions in water bodies, especially in the development of Chl-a estimation models from L-8 satellite OLI scenes (e.g., [56,57,58]). Additionally, our study confirms the essential role of a rigorous atmospheric correction to mitigate errors from aerosols and the sun-glint effect during the construction of spectral inversion models from satellite imagery. As discussed in Yépez et al. [11], the utilization of the USGS LaSRC SR_QA_AEROSOL product aids in filtering out pixels affected by atmospheric distortions, improving the relationship between in-situ data and image reflectance, even in cases of temporal disparities between data collection and scene acquisition.

II.

Development of an empirical retrieval model

The Green Chlorophyll Index (GCI) was selected as the optimal choice for constructing the Chl-a retrieval model in the Laguna Grande, despite being a spectral index used in terrestrial applications [49]. Our results show a remarkable sensitivity of Chl-a concentration to various multiband reflectance combinations, underlining high responsiveness of this index. Similar results have been observed in other water quality studies, where spectral indices such as NDVI and GNDVI have been distinguished, as reported in previous works (e.g., [59,60,61]). The CGI showed the highest relationship with R² value of 0.88, making it the best choice for the Chl-a recovery model in our study lake. The robustness of the model was assessed by bootstrapping method, confirming its statistical significance. When comparing our results with previous work, a notable study that shows a relationship between the green and near-infrared bands and in-situ Chl-a values is the research conducted by Yang et al. [53]. In their research, they revealed a substantial correlation, with the B5 (near-infrared) / B3 (green) spectral ratio index derived from the OLI sensor showing the strongest association with in-situ chlorophyll-a levels (R² = 0.73).

The analysis of turbidity parameters showed that the Red + NIR spectral index achieved the highest R² value, reaching 0.79. However, a high MAPE of 54.89 % revealed model weaknesses, indicating a lack of robustness. The data also showed seasonal variations, which highlights the importance of collecting winter data (water samples) to overcome the limitation of this retrieval model. Thus, to account for the seasonality factor comprehensively, it is advisable to incorporate data from several seasons for model training.

III.

Assessment of the spatial-temporal variability

Our results clearly show a spatial pattern varying from higher Chl-a values in the northern zone compared to the southern zone. These results are in agreement with the research of Pedrero-Guarda et al., [23], which identified spatial variations in temperature in Laguna Grande, possibly attributed to a heat island effect. This effect causes temperature differences, with higher temperatures near urban areas (northern part) and lower temperatures in vegetated regions. Since Chl-a levels are affected by factors such as nutrients, phosphorus, pH, water flow and temperature, it is plausible that this temperature gradient contributes to higher chlorophyll-a levels in the northern part of the lake. A positive aspect of our analysis arises when evaluating Chl-a values estimated from OLI images between 2014 and 2022. In general, a decrease in Chl-a concentrations is observed throughout the lake, which can be attributed to recent local regulations aimed at reducing activities that cause lake pollution. However, it will be necessary to evaluate a time series on a monthly scale to determine whether this decline has been constant and is being maintained.

This study marks the first step towards the development of a web platform that will provide access to water quality data derived from satellite signals, such as Landsat-8 OLI, for the main lakes of Chile. In our ongoing analysis at the Laguna Grande, we expect to improve our spectral estimation models by incorporating advanced airborne hyperspectral technology using NEO's HySpex VNIR-1800 cameras and collecting in-situ data with other water quality parameters. The results of this study will encourage the utilization of earth observation data for monitoring inland water quality, facilitating the sustainable management of water resources.