Satellite remote sensing: trophic state, spatio-temporal distribution of chlorophyll-a, and early warning of phytoplankton blooms 

Prepared by Javier Maldonado Alfaro (IIGEO/UMSA)

& Xavier Lazzaro (BOREA/IRD)

1 - Validation of Atmospheric Correction (AC) 

This validation is crucial to correctly estimate the algorithms of chlorophyll-a and phycocyanin concentrations, based on the in situ validation of satellite images of the qualitative spatial distribution of these pigments.  We did not created a specific algorithm, but compared the histograms of the reflectance values of 11 existing atmospheric corrections for each band. We elaborated a box plot from 542 recording points on transects on 4 different dates: on April 22, May 22, June 11, and August 10, 2019. We noted that the original reflectances, without atmospheric correction (1C), presented higher values than most of the bands with atmospheric correction, with the exception of C2RCC_rtoa and ACOLITE_rhot. There are 5 histograms with a similar pattern of intermediate reflectance values: 2A, iCOR 60 m, ACOLITE_rhos, ACOLITE_rhow, and HH2. There were two histograms with very low reflectance values: C2RCC_rrs and ACOLITE_rrs. Yet, ACOLITE_rhos and ACOLITE_rhow were very similar.

We use several statistical metrics such as: 'MAPD = Mean Absolute Percentage Difference'; 'RMSD =  Root Mean Square Deviation'; 'MBE =  mean bias error (it captures the average deviations between two datasets; it has the units of the variable; values near 0 are the best, negative values indicate underestimation and positive values indicate overestimation); 'R2 = Coefficient of Determination, that provides information about the goodness of fit of the regression line; 'RMSE = Root-Mean-Square Error (measures the average difference between values predicted by a model and the actual values; and 'NRMSE = Normalized Root Mean Square Error (it calculates the absolute value between predicted and observed values using different type of normalization methods). The bands presenting the best metrics were:

•⁠  ⁠For 'MAPD', the candidates were '2a', 'iCOR_60 m', 'C2RCC_rhown', and 'C2RCC_rms'. 

•⁠  ⁠For 'RMSD', the only candidate was 'iCOR_60 m'.

•⁠  ⁠For 'Mean Bias', the only candidate was '2a'.

•⁠  ⁠For 'R2', the two candidates were 'iCOR_60 m', and 'C2RCC_rss'.

•⁠  ⁠For 'RMSE', it was only 'iCOR_60 m'.

•⁠  ⁠For 'NRMSE', it was' iCOR_60 m', by far in all bands.

Therefore, 'iCOR_60 m' was the atmospheric correction that most closely resembles HH2 (i.e. the in situ HandHeld spectroradiometer). It was thus selected to perform the following atmospheric corrections in the Titicaca Lago Menor.

2 - Concentration algorithm for chlorophyll-a

We obtained the values of the statistical metrics from transects performed in Lago Menor on April 26, 2022 (Fig. 1) and May 1st, 2022 (Fig. 2), which correspond to the start of the dry season. Along those transects, we measured the chlorophyll-a concentrations ([Chl-a], Table 1) using the YSI EXO2 multiparametric probe paired with an EXO Handheld Display that registered the GPS coordinates every second. 

The reflectances of bands B4 to B8A had a significant positive correlation with the chlorophyll-a concentration [Chl-a] only on the April 26, 2022 transect. In contrast, on the May 1st, 2022 transect there were no significant correlation, thus chlorophyll-a concentration [Chl-a] did not influence reflectance, but other parameters may have, such as turbidity for example. 

The reflectances of the B8A/B4, B8/B4, and B2/B3 band ratios present negative significant correlations with the chlorophyll-a concentrations [Chl-a] , whereas the B3/B2 band ratio presents a significant positive correlation. 

Algorithm implementation - Using Rstudio, we performed a multiple regression as 'Steptwise regression' between the predicted chlorophyll-a concentration [Chl-a] and the reflectances of all bands and all band ratios. As such, we got 12 predictive models where [Chl-a] is the dependent variable (Table 4). 

The three 'best' models are those using the LN [Chl-a] (neperian logarithm), with a p-value < 0.1. The R2 values of models R1 and R3, R2 and R4, R5 and R7, R6 and R8, R9 and R11, R10 and R12 are equal to each other; that is, models with p-value < 0.1 and with p < 0.05 are equivalent. Therefore, it is better and sufficient to use only the models with p-value < 0.05. That is, only models R3, R4, R7, R8, R11 and R12. 

Choosing the 'best' model - From the remaining 6 models with p-value < 0.05, the best model is chosen according to the statistical metrics criteria: i.e., highest R2 value, highest adjusted R2 , the lowest AIC (i.e. Akaike information criterion, an estimator of prediction error or relative quality of statistical models for a given set of data),  and the lowest RMSE (i.e. the root-mean-square error, a measure of the differences between true or predicted values). Thus, for RAW [Chl-a], the candidate is R7 (Table 5), and for LN [Chl-a] it is R8 (Table 6). 

Validation of the predictive models into a satellite [Chl-a] algorithm - The satellite image with the most data points on the transects (Fig. 5) was chosen, which corresponds to the R4 model of April 26, 2022 (Table 6). Meanwhile, the R12 model was not far behind.  

The image in Fig. 5 was created with the "R4_Exp_JAMA_XL_2023" algorithm whose equation has this form: 

[Chl-a, µg/L] = exp (- a + b *(B3/B4) + c *(B4/B5) 

- d *((B2-B4)/B3) - e *B12 + f *B2 - g *B3 + h *B5)

where a, b, c, d, e, f, g and h are the coefficients of the bands and band ratios.

Note: The values of the coefficients are not given here, because we have observed that they change with meteorological and environmental conditions, which modify the reflectance. In particular, with the current intense drought phenomenon*, it is necessary to re-validate the coefficients of this equation with new in situ measurements of [Chl-a] at the HydroMet Buoy site, as well as at the representative limnological stations of the "Minimum Network" defined in the UNDP 05-B-05 OLT Pilot Project.   

*It results in increased suspended solids, organic matter, phytoplankton biomass and [Chl-a], reduced transparency and solar radiation penetration (PAR and UV). Together, they alter the water optical characteristics, and ultimately the reflectance. 

3 - OWT classification with Sentinel-3

To visualize the spatio-temporal evolution of the trophic states (not only [Chl-a]) exhibited over the entire extension of Lago Menor, we used the Optical Water Type classification (OWT) with the Sentinel-3 mission (resolution 300 m, daily pass over the lake), as part of the Copernicus program, funded by the European Space Agency (ESA). This classification is based on the characteristics of the spectral response of the water surface remote sensing reflectance.  We specifically used the GLaSS_6C classification of ESA that comprises 6 OWT Classes based on  [Chl-a] in µg/L; CDOM (Coloured Dissolved Organic Matter) in RFU (Relative Fluorescence Unit), and TSM (Total Suspended matter) in g/m3. We used the C2RCC atmospheric correction (AC). Sentinel-3 images are located at the same site as Sentinel-2 https://scihub.copernicus.eu/dhus/#/home , and the procedure is the same. As examples, we present here the OWT maps on four occasions:

On April 12, 2022 (beginning of the dry season), Class 1 predominates (90.81%) throughout the pelagic zone of Lago Menor. This class is representative of clear waters slightly affected by [Chl-a] between 0.1 and 12.3 µ.g/L (oligo-mesotrophic). Five sectors are highlighted, where classes 2, 3 and 4 are present. These classes are representative of zones with high [Chl-a], ranging from 0.8 to 69.6 µg/L (class 2, meso-eutrophic), 1.3 to 33 µg/L (class 3, eutrophic) and 0.9 to 705 µg/L (class 4, hyper-eutrophic). Bahía Cohana has one of the highest [Chl-a], where near the mouth of the Katari River, a medium class 4 zone is observed, then a smaller class 3 zone, and a even smaller class 2 zone. This shows a strong gradient of eutrophication along the littoral zone of the North-East region, most populated and contamined by the residual waters transported through the Katari watershed from the El Alto urban areas. 

On April 27, 2022 (beginning of the dry season), the proportion area of Class 1 decreased compared to April 12, from 90.8% to 81.7%. As a consequence, the other classes expanded in the five zones, especially Class 3 from 1.3 to 33.0 µg/L. Cojata Island (Sector 1) stands out, where Class 3 extended to Cumana Bay, passing through Puerto Pérez. Classes distributions changed as well in Huarina and Chilaya, where Class 3 replaced the clear waters of Class 1.

On May 1st, 2022, the contribution of Class 1 slightly increased compared to the previous date, from 81.72% to 86.02%. In Yunguyo (Sector 5), Class 3 was fragmented, while Classes 2 and 4 developed. In Bahía Cohana, Class 3 extended towards the Taraco Peninsula. In the Cojata Island Sector 1, especially  the Class 3 extension decreased, whereas it expanded as far as Cumana Bay.

These four images over a 1.5 month period mid-April to end of May) demosntrate how fast the different OWT classes expanded in relation of each others. Class 1, representative of clearer waters poorly affected by [Chl-a] , dominated in the pelagic zone over the 4 dates with proportions varying from 90, 82, 86 and 94 %, respectively. 

More specifically, on April 27, 2022, Class 1 decreased, whereas Classes 2, 3 and 4 increased. These increases occur in Cojata island, Cohana bay, Tiahuanacu river mouth, and offshore of Zepita and Yunguyo communities. These Classes represent eutrophic to hyper-eutrophic conditions with ranges between 0.8 and 69.6 µg/L (Class 2), 1.3 and 33.0 µg/L (Class 3), and 0.9 to 705.0 µg/L (Class 4). The decreasing pattern of Class 1 and increasing patterns of Classes 3 and 4 are mainly observed near areas occupied by human populations (Yunguyo, Zepita) and in river discharge areas (Cohana bay, Cojata island, Tiahuanacu river mouth). 

This satellite remote sensing tool reveals the dynamics of eutrophication (increase in [Chl-a]), and identify areas affected by human activities, such as agriculture, aggregate movement, sewage discharge, industrial pollution. Two great advantages of the OWT GLaSS_6C classification are that: (1) it can be used even in absence of in situ water quality measurements, such as [Chl-a], and (2) using Sentinel-3, it can generate images on a daily basis. Therefore, this OWT classification provides a real-time monitoring of the trophic state evolution of the lake, and as such observing the dynamics of the different classes, especially the more eutrophic Class 4, can provide an Early Warning of phytoplankton bloom forming, identifying the more vulnerbale areas to be monitored and restored. 

To test the power of the GLaSS 6C classification as a surveillance and monitoring tool, we applied this procedure at the scale of the entire Lake Titicaca. As such, the predominance of Class 4 (dark blue, with higher [Chl-a]) is marked in areas of important human activities (Fig. 10). Thus OWT maps are particularly useful not only for the research activities on ecological processes, but also for the management of water resources and the monitoring of areas undergoing pollution and eutrophication, e.g. the anticipation of phytoplankton blooming phenomena that are very harmful.