Retrieving Water Turbidity in Araucanian Lakes (South-Central Chile) Based on Multispectral Landsat Imagery

: Remote sensing was used as an early alert tool for water clarity changes in ﬁve Araucanian Lakes in South-Central Chile. Turbidity records are scarce or unavailable over large and remote areas needed to fully understand the factors associated with turbidity, and their spatial-temporal representation remains a limitation. This work aimed to develop and validate empirical models to estimate values of turbidity from Landsat images and determine the spatial distribution of estimated turbidity in the selected Araucanian Lakes. Secchi disk depth measurements were linked with turbidity measurements to obtain a turbidity dataset. This in turn was used to develop and validate a set of empirical models to predict turbidity based on four single bands and 16 combination bands from 15 multispectral Landsat images. The best empirical models predicted turbidity over the range of 0.3–12.3 NTUs with RMSE values around 0.31–1.03 NTU, R 2 (Index of Agreement IA) around 0.93–0.99 (0.85–0.97) and mean bias error (MBE) around ( − 0.36–0.44 NTU). Estimation maps to analyze the temporal-spatial turbidity variation in the lakes were constructed. Finally, it was found that the meteorological conditions may affect the variation of turbidity, mainly precipitation and wind speed. The data indicate that the turbidity has slightly increased in winter–spring. These models will be used in the future to reconstruct large datasets that allow analyzing transparency trends in those lakes. To analyze climatic factors inﬂuencing variations in


Introduction
The characteristics of water regulate the metabolism of lakes; the modifications of the aquatic environment are produced as a response to climatic and geographical variations [1]. Water quality monitoring in continental aquatic ecosystems using satellite image processing is a tool that can be employed as an early alert system for changes in lakes [2]. One important characteristic in an aquatic environment is the water clarity. Generally, the clarity is quantified by the Secchi disk depth (SDD) or turbidity parameter and is mainly caused by increases/decreases in the concentration of suspended sediments, algae, and organic matter [3,4]. Reductions in water clarity can inhibit sunlight absorption in lakes, thereby slowing radiant heating of lakes, reducing light penetration, and primary productivity [5]. Halted or reduced primary productivity means a decrease in plant survival and dissolved oxygen output. Other impacts are the regrowth of pathogens [6], habitat quality [7], and water bodies' recreational use [8]. In Landsat images, high concentrations of algae are associated with decreased radiance because of the chlorophyll absorption in this region [9].
Turbidity is an optical characteristic of water clarity [1] and is a measurement of the amount of light that is scattered by particles in the water column. When higher the intensity

Sampling Measurements
Seven monitoring campaigns during the 2015 year were conducted by the EULA Center of the University of Concepcion. Therefore, all seasons were covered. The monitored parameters were surface temperature (2250 B standard methods 22 Ed. thermometry, method used as a reference for temperature analysis according to NCh 2313 compendium), surface Chl-a (fluorometric method), transparency (measured by SDD), profiles Chl-a and temperature, (Multiparametric Probe Seabird 19 Plus). Data were collected from different sampling stations for each lake in addition to stations in the respective tributaries (tributaries and effluents) according to the characteristics of these aquatic ecosystems showed in Figure 2. The locations were selected in the field considering the following criteria: 1.
Morphology of the lake.

2.
Presence of tributaries (away from their influence).

3.
Presence of industrial effluents or urban discharges.
At each of the lake stations, water samples were collected at 5 depths using a 5 L Niskin sampling bottle. The samples were stored and transported in thermally insulated boxes, duly refrigerated with ice at a temperature of approximately 5 • C for subsequent analysis. The chemical analyses were carried out in the chemistry laboratory of the EULA-Center Chile. This laboratory accredited by the National Institute of Normalization for the Chilean Norm NCh ISO 17.025 of 2005.

Depth.
At each of the lake stations, water samples were collected at 5 depths using a 5 L Niskin sampling bottle. The samples were stored and transported in thermally insulated boxes, duly refrigerated with ice at a temperature of approximately 5 °C for subsequent analysis. The chemical analyses were carried out in the chemistry laboratory of the EULA-Center Chile. This laboratory accredited by the National Institute of Normalization for the Chilean Norm NCh ISO 17.025 of 2005.

Image Acquisition and Processing
A total of 15 multispectral Landsat images, 13 by L8 OLI (Operational Land Imager) and 2 by L5 TM (Thematic Mapper) with a 30 × 30 m resolution, obtained from the official site of the United States Geological Survey (USGS) Earth-Explorer (https://earthexplorer. usgs.gov/) accessed on 20 December 2020 were used. The study area is covered by three scenes with the following path/row: 233/88, 233/87 and 232/88. The selected images were chosen according to quality, availability, a low percentage of cloud cover (0-12%) and proximity to the sampling date (within a ±14 days difference of satellite overpasses). To mask clouds, cirrus and shadows, the quality assessment (QA) bands were used and confirmed with a visual inspection. Only data that was outside of cloud cover were used. The lake's edge was established using geospatial information from the DGA. Although the satellite images covered an area greater than the lake, only the water body was considered in the analysis.
Images preprocessing was carried out using ENVI 5.3 (Environment for Visualizing Images), ArcGIS (ESRI's v. 10.8.1) software tools and began with geometric correction (reprojection to the UTM-WGS84 coordinates system, zone 19S) [26]. Radiometric calibration is a common pre-processing step that attempts to compensate for radiometric errors from sensor defects, variations in scan angle, and system noise to produce an image that represents the true spectral radiance at the sensor [27][28][29]. Therefore, radiometric calibration to L5 TM images was carried out in accord with [28,29] and to L8 OLI images, across the instructions proposed by the USGS [30] in Landsat 8 Data Users Handbook v.5, converting the original digital number (DN) of each band to radiance values and then transforming them to top-of-atmosphere (TOA) reflectance values, that include a correction for solar angle and use of the rescaling factors and parameters found in the MTL.txt file [30,31]. The atmospheric correction continued with the Dark Object Subtraction (DOS) technique [32]. This method is based on some pixels in the image being in complete shadow and their radiances received at the satellite are due entirely to atmospheric scattering. Therefore, this value is subtracted from each pixel value in the image [31]. DOS is widely used in research applied to the water quality of aquatic ecosystems with remote sensing, such as [31,33,34]. With the corrected bands the calculation of simple bands, spectral indices and band combination was selected according to a literature review [5,[33][34][35]. Table 2 present the Landsat images selected and in situ measurements date.

Standardization of Existing Water Clarity Measurements
Field measurements of water clarity via SDD were acquired by the DGA between 2011 and 2012 (Calafquén, Panguipulli, Riñihue) and 2012-2013 (Puyehue, Neltume). Measurements of water clarity via SDD and NTU were acquired on 25 January 2011 (Calafquén, Panguipulli, Riñihue) (n = 28) and 8 January 2013 (Puyehue, Neltume) (n = 28). To standardize the different measures of water clarity, we developed a linear function between SDD and NTU measures using the samples collected in January corresponding to 2011 and 2013 following the methodology proposed by [5]. The resulting function was used to convert in situ SDD measurements into correlative NTU.

Band Combinations, and Turbidity Index
Correlative NTU values were correlated to the corresponding cell value of four (n = 4) independent bands (blue (B), green (G), red (R) and Near Infrared (N)) and sixteen (n = 16) different band combinations consisting of bands B to N, including the Normalized Difference Turbidity Index (NDTI) which is calculated as (R − G)/(R + G) [5]. A simple regression to develop a relation between in situ NTU measurements and space-based observations was used [33,36]. These bands and band combinations are regularly employed to detect turbidity in water bodies [4,5,37]. The Pearson correlation coefficient (r) was used to evaluate the correlation between NTU values and 20 independent bands, band combinations and NDTI.

Empirical Models and Validation
In situ datasets and the preprocessed Landsat 8 images (11) during the monitoring campaigns in 2015 were used to develop a set of empirical models for obtaining turbidity or NTU. To reduce the possible errors in geometric correction of Landsat images and the dynamics of water bodies, we used a sampling window with the 3 × 3 pixel mean to extract the surface reflectance of the bands, band combinations, and NDTI spectral index (i.e., the preprocessed Landsat data) according to [4]. The data selection for internal validation/training and external validation/prediction was conducted according to the methodology proposed by [38], which was the most useful of the several data selection strategies that were tested in this study. The resulting data set consisted of the NTU response variable (from the SDD) and best band or band combinations predictors for each lake of the 4 selected images (2011-2013), resulting in 4 × 3 matrices. The strategy consists of randomly selecting 70% of the data as the calibration series, with the remaining 30% serving as the validation series. The models were evaluated based on the coefficient of determination (R 2 ), index of agreement (IA), mean bias error (MBE) and root mean square error (RMSE). The RMSE is frequently used to compare the forecasting errors of different NTU models [37,39,40]. In general, the RMSE quantifies the dispersion between simulated and measured data, where low RMSE values indicate better model performance. Negative (Positive) MBE indicates that the model results underestimate (super-estimate) measurements. An IA value close to 1 indicates a more efficient model [41]. All statistical analysis and figures in this study were realized in Origin Pro 2021 version 9.8.0.200 (Academic) software. Finally, estimated turbidity maps were obtained from the models, using ArcGIS 10.8.1 software.

Water Clarity and Meteorological Conditions
Seasonal meteorological factors such as wind speed (WS) and precipitation impact water turbidity [42], due to their ability to affect the movements of sediments. For example, in Lobo reservoir [43] it was observed that maximum particulate and organic matter resuspension is caused by turbulence on days of high wind speed. In reference [44] it is said that the turbidity of the surface layer is attributed to the decrease in incoming water caused by the decrease in rainfall, which, in turn, causes a decrease in the buoyancy of the hydrothermal plumes. In addition, it could cause a reduction in the supply of sediment from runoff. To analyze climatic factors influencing variations in Remote Sens. 2021, 13, 3133 8 of 16 estimated turbidity, meteorological data on precipitation and wind direction/speed were used. The meteorological data analyzed were obtained from Dirección Meteorológica de Chile (DMC, https://climatologia.meteochile.gob.cl/) accessed on 20 December 2020 and Instituto Nacional de Investigaciones Agrarias (INIA, https://agrometeorologia.cl/) accessed on 26 December 2020. The wind direction/speed and precipitation data (e.g., on the day, on the last 3 days when the satellite image was captured) were taken from stations close to study lakes, Santa Carla station (−39.67 • S, 72.61 • W, 264 m.a.s.l) representative to (Calafquén, Neltume, Panguipulli, and Riñihue lakes) and Rupanco or Rucatayo stations (−40.74 • S, −72.66 • W, 272 m.a.s.l) representative of the conditions of Puyehue lake (see Figure 1). Table 3 shows the average of SDD, Chl-a, NTU and temperature values obtained in the monitoring campaigns during 2015. The behavior of temperature fluctuated between 14.9 ± 5.6-16.1 ± 4.4 • C with a maximum value of 25.9 • C in Neltume lake. Water clarity during all campaigns was measurement via SDD and indicate that the average fluctuated between 11.9 ± 5.8 m (Calafquén) to 6.8 ± 2. The NTU values during 2015 were estimated through the relationship between SDD and NTU. The average values of NTU ranged from 1.5 ± 1.0 NTU until 9.2 ± 1.6 NTU for Puyehue and Neltume lakes. Hence, the trophic state of the studied lakes can be classified as oligotrophic or ultra-oligotrophic, with low nutrient levels and productivity.

In Situ Secchi Disk Depth (SDD) Measurements into Correlative Nephelometric Turbidity Units (NTU)
A strong negative correlation (r > −0.92) was found between SDD measurements and simultaneous NTU measurements (see Table 4). This relation allowed us to derive NTU values for all SDD data collected during other sampling periods.  Figure 3 shows the Pearson correlation between the turbidity variable and 4 band/ 16 band combinations/NDTI from the processing of the 11 satellite images. Significant band correlations included various combinations of the blue, green, red and near-infrared bands. For Calafquén lake, twelve significant turbidity relationships were ranging from −0.67 to 0.97 (p-value ≤ 0.001). Neltume and Panguipulli lakes had 16 significant correlations (p-value ≤ 0.001) and three (p-value ≤ 0.01) ranging from −0.67 to 0.88 and 0.60 to 0.85, respectively. Riñihue lake had four significant correlations with a significance level (p-value ≤ 0.001) ranging from −0.76 to −0.83 and three correlations of (p-value ≤ 0.01). Meanwhile the Puyehue lake only had two significant correlations r = −0.91 (p-value ≤ 0.001). For more detail on the significance level look at Figure S1.  The performance of the five best models was evaluated based on determination coefficient (R 2 ), RMSE, IA and MBE as shown in Table 5. The empirical models showed highest R 2 and IA values ranged between 0.93 to 0.99, and 0.85 to 0.97 respectively, which The performance of the five best models was evaluated based on determination coefficient (R 2 ), RMSE, IA and MBE as shown in Table 5. The empirical models showed highest R 2 and IA values ranged between 0.93 to 0.99, and 0.85 to 0.97 respectively, which indicates a good score for all models. An inspection of Table 5 reveals that the model underestimates the NTU measurements at Calafquén, Neltume and Puyehue lakes and super-estimate in Panguipulli and Riñihue lakes. However, the MBE values were smallest. Similarly, the RMSE values were lowest ranged between 0.31 NTU to 1.03 NTU. Therefore, the estimated NTU values are close to the observed values.  Figure 4 presents the scatter plots comparing satellite-derived and in-situ measured turbidity. The values showed a high coefficient of determination (R 2 > 0.93), suggesting that all models predict NTU values very well. Therefore, the models can be used as an early warning tool for changes in water clarity. Indeed, reached a better or similar R 2 than other regression equations for the NTU models from Landsat-8 shown by [5,37,39].   Figure 5 indicates the spatial distribution of turbidity derived from Landsat 8, with NTU estimated by developed empirical models of each lake in two different seasons (sum-  Figure 5 indicates the spatial distribution of turbidity derived from Landsat 8, with NTU estimated by developed empirical models of each lake in two different seasons (summer and spring of Southern Hemisphere). The value of NTU is representative of the upper layer (euphotic zone, 90% of the incident radiation).  In general, turbidity levels are low, with an increase in the winter and spring seasons. In all lakes the turbidity increments were different (see Figure 5). It can note that the highest turbidity is reported in Neltume lake with values in the range of 10.3-12.3 NTU). This could be because it is the shallowest lake (86 m depth), with the smallest area (9.86 km 2 ) and it has four tributaries that contribute sediment, increasing turbidity. While the lowest NTU values are reported for Puyehue lake (0.8-1.3 NTU).

Influence of Meteorological Conditions in the Study Area
For a better analysis of physical forcing and their interplay in the distribution of turbidity at the local scale of the five lakes, we compared a short time series of turbidity, precipitation, WS and direction in the period covered by Landsat 8 imagery. The meteorological variables study was taken from stations Santa Carla representative to (Calafquén, Neltume, Panguipulliand Riñihue lakes) and Rupanco station representative of the conditions of Puyehue lake (see Figure 1). Figure 6 shows the behavior of monthly values of accumulated precipitation during 2015. The maximum accumulation rate of precipitation occurs between the winter and spring months (rainy seasons), while the minimum accumulated precipitation values were obtained in the summer months (0 mm, dry season). The highest accumulated value (444.5 mm) during 2015 was reported in the month of July (see Figure S2) in the area. The highest turbidity values correspond to the rainy period (winter and spring seasons) and the lowest to the dry period (summer season). The precipitation variable is directly correlated (R 2 = 0.91) with the turbidity reported in the study lakes (with p-value ≤ 0.001). In general, turbidity levels are low, with an increase in the winter and spring seasons. In all lakes the turbidity increments were different (see Figure 5). It can note that the highest turbidity is reported in Neltume lake with values in the range of 10.3-12.3 NTU). This could be because it is the shallowest lake (86 m depth), with the smallest area (9.86 km 2 ) and it has four tributaries that contribute sediment, increasing turbidity. While the lowest NTU values are reported for Puyehue lake (0.8-1.3 NTU).

Influence of Meteorological Conditions in the Study Area
For a better analysis of physical forcing and their interplay in the distribution of turbidity at the local scale of the five lakes, we compared a short time series of turbidity, precipitation, WS and direction in the period covered by Landsat 8 imagery. The meteorological variables study was taken from stations Santa Carla representative to (Calafquén, Neltume, Panguipulliand Riñihue lakes) and Rupanco station representative of the conditions of Puyehue lake (see Figure 1). Figure 6 shows the behavior of monthly values of accumulated precipitation during 2015. The maximum accumulation rate of precipitation occurs between the winter and spring months (rainy seasons), while the minimum accumulated precipitation values were obtained in the summer months (0 mm, dry season). The highest accumulated value (444.5 mm) during 2015 was reported in the month of July (see Figure S2) in the area. The highest turbidity values correspond to the rainy period (winter and spring seasons) and the lowest to the dry period (summer season). The precipitation variable is directly correlated (R 2 = 0.91) with the turbidity reported in the study lakes (with p-value ≤ 0.001). The wind forcing values were taken in a three-day window before the Landsat 8 acquisitions, and are represented in the four diagrams in Figure S3. The wind speed at the Santa Carla station was 6-8 km/h, while at the Rupanco station was 4-6 km/h. The Santa Figure 6. Accumulation monthly rate of precipitation, average monthly wind speed, and temperature in all study lakes and estimated turbidity.
The wind forcing values were taken in a three-day window before the Landsat 8 acquisitions, and are represented in the four diagrams in Figure S3. The wind speed at the Santa Carla station was 6-8 km/h, while at the Rupanco station was 4-6 km/h. The Santa Carla wind rose ( Figure S3) shows that the typical Puelche wind that blows from the N to the NE is the predominant wind characteristic of South-Central Chile, followed by winds from the NW sector in spring ( Figure S3). Furthermore, the wind analyses for Rupanco ( Figure S3) clearly show that the N winds for this season are the strongest and occur with the highest frequency during 80% of the time in the investigated period. The wind increase shows a bimodal distribution between the NE and NW winds except for Rupanco in the summer season. The diagrams are quite different, mainly due to the different capacities of each station to detect winds. This is essentially related to geographical factors. We do not find a good correlation between wind direction and estimated turbidity, but we obtain a good correlation with the WS (R 2 = 0.75). Similar results were reported by [4] in Diefenbaker lake.

Discussion
Water clarity is an important aspect of the freshwater lake system. Increases or decreases in water clarity can negatively affect the biological component of the system which can be adapted to specific light penetration conditions [9]. Remote sensing was used as an early alert tool for water clarity changes in an inland aquatic ecosystem in South-Central Chile. Particularly in Chile, this study and methodology are relevant because there are not enough data to assess the state of water resources, mainly in lake systems https://snia.mop.gob.cl/ accessed on 20 January 2021. The monitoring lake network only covers 20 of the 375 lakes with areas larger than 3 km 2 [22]. On the other hand, due to the climatic conditions of the south, the satellite images captured in the winter period present extensive cloud cover that makes it impossible to acquire specific data from the North Patagonian aquatic ecosystems. Hence the marked importance in developing empirical models for estimating parameters of water quality, applicable for any period of the year.
Water turbidity was modeled in five lakes in South-Central Chile through seven monitoring campaigns during 2015. From 15 multispectral Landsat images (13 L8 OLI, 2 L5 TM) with a resolution of 30 × 30 m, turbidity index values were predicted based on four simple bands (green, blue, red and infrared) and a combination of 16 bands. For modelling we used 100 records of NTU as response variable and 20 spectral band combinations as predictor variable. Linear regression models were constructed using the best correlation between NTU and band combination. The empirical models were evaluated using four statistical indicators (R 2 , RSME, IA and MBE) to select the best model per lake. The models present a good agreement between NTU modeled and in situ measurements with low RSME values. Indeed, they had higher coefficients of determination and IA compared to those used in other studies [17,42,45]. The reported MBE indicates a low bias for all models.
Thanks to the unprecedented spatial and radiometric resolution of the L8 sensor, we then mapped the turbidity estimated in 2015 to use as an early warning tool for changes in water clarity. Turbidity for the Araucanian lakes studied is low (range between 0.3-12.3 NTU), with marked seasonal differences of a slight increase in winter and spring. The lake with the highest turbidity was Neltume lake (12.3 NTU, spring season). This behavior could be influenced by the hydrodynamic condition of these lake systems, where they present a summer stratification, for which the euphotic surface layer is separated by a temperature gradient from the deeper or lower layer (aphotic), which prevents the circulation of the water column during this season of year. Also, in the summer months (higher surface temperature) the rainfall regime is lower than in the other seasons of the year. After the breakdown of the stratification, a gradual mixture of the water column would be generated in the autumn-winter and spring seasons, where a hydrodynamically mixed column of water is present, with the suspension of particles such as total or dissolved solids and other contributions of the basin, which allows a greater movement of the water column by internal circulation, and in these months greater contributions are evidenced due to a higher rainfall regime and increased surface runoff [42].
The maximum accumulation rate of precipitation occurs in the winter and spring months, coinciding with the highest values of NTU. In both meteorological stations, the precipitation is greater towards the first and last days of every month (all seasons), which coincides with the dates of the processed Landsat images; therefore, this shows the influence of the precipitation on the turbidity behavior. To know how meteorological conditions contribute to the variability of turbidity were correlated with precipitation, wind speed and wind direction. The precipitation and wind speed showed high correlation (R 2 = 0.91 and 0.75) with the turbidity reported in each lake. Precipitation in conjunction with wind speed may have induced sediment resuspension during the rainy season which in turn contributed to turbidity in the upper reaches of lakes.
It is relevant to point out that the remote-sensing tool allows monitoring of representative limnological variables of the trophic state in the lakes and improves the spatial representativeness of the measurements [14,16,17]. In the planning and management of the water quality of the oligotrophic lakes studied here, it is important to consider the precautionary principle, that is, to have as an objective the monitoring of these lakes that have an oligotrophic water quality to prevent them from evolving irreversibly to a deteriorated water quality. Recently, there has been evidence that Panguipulli lake is showing the first signs of eutrophication [16], this being particularly noticeable in an increase of microalgae biomass in the lake bay, which coincides with a deterioration of water and sediment quality in this bay. Although the lakes have good water quality conditions, we recommend gradually reducing the pressure of anthropogenic activities around the basin and in the lake, to maintain what are still good water quality conditions. We suggest gradually reducing the pressure of anthropogenic activities, for example, land use change, considering scenarios of climate change and water deficit in this region [46,47]. Decreased precipitation and deforestation could mean a decrease in base flows in these lakes.
In situ turbidity measurements and monitoring campaigns are often costly and timeconsuming, limiting their spatial and temporal representation in large and remote lakes. In this way, remote sensing processing and its relation to in situ measurements allows modeling to provide early information on water quality changes in aquatic ecosystems. The use of the models and maps developed in this work can be applied to obtain a real-time assessment of the turbidity parameter representative of water quality and provide information to assess water clarity for governmental decisions and environmental protection planning.

Conclusions
To our knowledge, this is the first time that a model to estimate turbidity in Araucanian lakes has been used incorporating in situ measurement Secchi disk depth, and Landsat images, rather than using single bands and a band combination. The lakes studied have maintained their oligotrophy, related to a high quality of the water bodies and a low productivity present in the in-situ measurements. The accuracy of the empirical models was generally high, with RMSE values ranging from 0.31 to 1.03 NTU values over the range of 0.3 to 12.3 NTU in the validation data not used for model development. Therefore, the models can be used as an early warning tool for changes in water clarity. Spatio-temporal variation shown in turbidity maps was low with a slight increase during spring months. The assembled meteorological data (precipitation and wind speed) had a direct influence on water turbidity, following a seasonal cycle. The results provide information to evaluate the water quality for governmental decisions and environmental protection planning. As well as reconstructing large turbidity datasets, this allows analyzing transparency trends in those lakes.