An UAV and Satellite Multispectral Data Approach to Monitor Water Quality in Small Reservoirs

: A multi-sensor and multi-scale monitoring tool for the spatially explicit and periodic monitoring of eutrophication in a small drinking water reservoir is presented. The tool was built with freely available satellite and in situ data combined with Unmanned Aerial Vehicle (UAV)-based technology. The goal is to evaluate the performance of a multi-platform approach for the trophic state monitoring with images obtained with MultiSpectral Sensors on board satellites Sentinel 2 (S2A and S2B), Landsat 8 (L8) and UAV. We assessed the performance of three di ﬀ erent sensors (MultiSpectral Instrument (MSI), Operational Land Imager (OLI) and Rededge Micasense) for retrieving the pigment chlorophyll-a (chl-a), as a quantitative descriptor of phytoplankton biomass and trophic level. The study was conducted in a waterbody a ﬀ ected by cyanobacterial blooms, one of the most important eutrophication-derived risks for human health. Di ﬀ erent empirical models and band indices were evaluated. Spectral band combinations using red and near-infrared (NIR) bands were the most suitable for retrieving chl-a concentration (especially 2 band algorithm (2BDA), the Surface Algal Bloom Index (SABI) and 3 band algorithm (3BDA)) even though blue and green bands were useful to classify UAV images into two chl-a ranges. The results show a moderately good agreement among the three sensors at di ﬀ erent spatial resolutions (10 m., 30 m. and 8 cm.), indicating a high potential for the development of a multi-platform and multi-sensor approach for the eutrophication monitoring of small reservoirs. ﬀ ective in discriminating chl-a values with this sensor. Three of the developed algorithms were combined in a preliminary ensemble model. Further research will be done to validate these models with di ﬀ erent ranges of chl-a data to apply them to intermediate and higher chl-a classes.


Introduction
One of the main problems currently affecting water resources is the eutrophication of water bodies, both coastal and freshwater [1,2]; it was already known ten years ago that around 20% of European lakes were suffering from nutrient enrichment [3]. In the European Union, the management of reservoirs and natural waterbodies is directly linked to the Water Framework Directive [4], which requires periodic monitoring; this being a common trait all over the world [5]. Most of the reservoir monitoring programs are currently based on measurements taken in the field making it difficult to capture the spatial and temporal variability of phenomena such as algal blooms, with irregular Figure 1. Localization of the study site in Europe (A), in the NW of the Iberian Peninsula (B) and aerial picture of the reservoir with in situ sampling points both from the periodic monitoring [28] and from the two Unmanned Aerial Vehicle (UAV) flights (same sampling points for 2017 and 2018) (C).

In Situ Data
The in situ chl-a data used to calibrate and validate the satellite models were taken from the Galician Platform for Environmental Information (GAIA), supplied freely and on a regular basis, coming from the cyanobacterial monitoring program in the reservoirs of the district [28]. It maintains 2 permanent sampling points in the waterbody, one close to the dock (BP) and the other near the The main problem facing the reservoir is eutrophication. It has a very small catchment area with agricultural and livestock farms that, in many cases, generate an excess of slurry responsible for the existence of diffuse organic pollution and an excess of nutrients in the water. It is a system affected by recurrent cyanobacterial blooms with periodic high chl-a concentration, especially in late summer-early autumn (Table 1), documented by the Environmental Information Service of Galicia (SIAM) since 2012 (www.siam.xunta.gal). In addition, its proximity to a highway makes it sensitive to potential spills from it. Taking this into account, and the sensitivity due to its role as the supply of drinking water and a leisure and recreational area, this reservoir is considered a very good example of a system with the need for preventive monitoring and periodic control. Table 1. Summary statistics of the chl-a values (µgr/L) registered in the two sampling points of the reservoir during the periodic monitoring during the 3 years of the study period. (n is the total number of samples available for each year and sample point BC and BP). Source: [28].

In Situ Data
The in situ chl-a data used to calibrate and validate the satellite models were taken from the Galician Platform for Environmental Information (GAIA), supplied freely and on a regular basis, coming from the cyanobacterial monitoring program in the reservoirs of the district [28]. It maintains 2 permanent sampling points in the waterbody, one close to the dock (BP) and the other near the uptake system for drinking water (BC) (Figure 1). For the purpose of this study, we established two matchup points for the satellite reflectance data retrieval, as near as possible as these two ones, in a place still considered representative of the chl-a values, while far enough from the dock and the wall dam to be as unaffected as possible from adjacency effects. The point near the dock (BP) is located in the direction of dominant winds and tends to accumulate a high concentration of cyanobacteria during bloom events. Under such conditions, this point was not considered representative of the free water pixel used for matchup with satellite images and not used in the study.
Chl-a concentrations were retrieved from filtered surface water samples followed by acetone extraction and fluorometric measurements [29]. All the samples were taken from 10:00 am to 2:00 pm (local time). The sampling period used for calibration spans from January 2017 to December 2017, with a dataset of 45 samplings and 86 potential matchup points. The available data for validation is from the previous and following year: 2016 (from 4th January to 19th December) and 2018 (5th March to 26th November) with a dataset of 34 and 45 potential matchups, respectively.
During the field data campaigns, simultaneously to the UAV flights (19th September 2017 and 2th October 2018), water samples were also taken in 5 points of the reservoir ( Figure 1). As the 90% of the radiation received by the remote sensors comes from the first optical thickness, which corresponds to 0.6 times the Secchi disk depth (SD) ( [24,30], integrated water samples were taken with a hosepipe from the surface layer at a depth equal to 0,6*SD. These water samples were then taken to the lab, stored in the dark in a cooler and analyzed for chlorophyll a (chl-a), phycocyanin (PC), Total Suspended Solids (TSS) and Dissolved Organic Carbon (DOC). The chl-a was determined by acetone extraction (90%) followed by spectrophotometric measurements and the application of Jeffrey and Humpfrey [31] equations following the procedure described in [32]. Phycocyanin determination was made by the method described by Quesada [33]. The method for obtaining TSS was nº2540D of the APHA Standard Methods [32] and the DOC concentration was obtained through combustion and catalytic oxidation followed by NIR with a Skalar Formacs HT TOC analyzer following UNE EN 1484:1998. Remote Sens. 2020, 12, 1514 5 of 33 Five profiles were also done with a YSI6600 multi parameter sonde with the following attached sensors: temperature ( • C), electrical conductivity (µS/cm.), pH/ORP, dissolved oxygen (% and mg/L), turbidity (NTU), chlorophyll (µgr/L) and phycocyanin (BGA-PC, cells/mL). The sample taken in the central point of the reservoir was also analyzed for principal nutrients (Total Phosphorous TP, Total Nitrogen TN, Nitrate, Nitrite and Ammonium). TP was analyzed by ICP-MS following an internal procedure of the support services for research (SAI) of the University of A Coruña (nº P-SAI-UEPM-09) accredited by ENAC, with high-resolution equipment (ELEMENTXR/ELEMENT2; Thermo Finnigan). Determination of TN was done through catalytic combustion and chemiluminiscence detection following prior oxidation with a Skalar Formacs HT analyzer in samples filtered through membrane (Millipore Millex-HM, 0.45 µm) following UNE-EN 12260:2004. Colorimetric ammonium determination was performed by AquaKem 250 (Labmedics) in water samples previously filtered through Millex-HN membrane (0.45 µm) following SAI internal procedure: P-SAI-UAA-15. Nitrate and nitrite were determined by ionic chromatography with a chromatographer Metrohm 850 Professional IC equipped with a column Metrohm Metrosep A Supp 7 250/4.0 mm, following the internal SAI procedure: P-SAIUAA-20.

Remotely Sensed Data
Three different multispectral sensors on board four platforms were used as well as three satellite platforms (S2A, S2B and L8) and one UAV platform (Octocopter Atyges FV8). The characteristics (spatial resolution and bandwidths) of the sensors (namely MSI, OLI and Rededge Micasense) are described in Figure 2 and Table 2).
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 of 33    All available images of both sensors (MSI and OLI) acquired during 2017 were considered for calibration. The first selection criterion was a maximum interval of 3 days with respect to the acquisition of in situ data. Images with cloud cover were discarded, achieving a final number of 21 a priori valid images for Sentinel 2A and 2B. Once processed and analyzed, those images taken during winter (December, January and February) were considered not valid because of the shadow projected from the hills besides the reservoir. The imagery of the rest of the year, captured under higher sun elevation angles, were not affected by this artifact. The final number of Sentinel 2 MSI valid images for 2017 were 12 (24 potential matchups). Landsat 8 OLI images were obtained following the same criteria as for Sentinel 2, with a final number of 7 valid scenes for calibration (14 matchup points).
The validation dataset with the previous and following year images (2016 and 2018) followed the same rules. For Sentinel 2 MSI, 4 images of 2016 and 7 images of 2018 were considered valid. The 2016 dataset has few images because only S2A was operational and the high cloud cover in this area; moreover, some of the cloud-free images were invalid by a very strong sun-glint effect. The S2 MSI validation dataset was made of 15 match-up points. In case of Landsat 8, the validation dataset was made of 5 images of 2016 and 6 corresponding to 2018, with a total number of 14 match-up points.
The Sentinel 2 MSI product used was L1C scenes in the SAFE format, i.e., orthorectified, geolocated and radiometrically calibrated top-of-atmosphere (TOA) reflectances in Universal Transverse Mercator (UTM) projection with the WGS84 datum. S2 MSI images were subsequently corrected to bottom-of-atmosphere (BOA) reflectance using the iCOR (Image Correction for atmospheric effects) atmospheric correction algorithm (AC processor), a single atmospheric correction implementation for land and water. It derives aerosol optical thickness from the image and allows for adjacency correction with SIMilarity Environmental Correction (SIMEC) over water. In iCOR, land pixels are used to derive the AOT based on an adapted version of the SCAPE-M method and uses MODTRAN5 LUT to perform the atmospheric correction. Solar and viewing angles (Sun zenith angle (SZA), view zenith angle (VZA) and relative azimuth angle (RAA)) and a DEM are used as auxiliary data and can originate from external sources or be derived from the image itself [34]. Default settings for this processor were used and the SIMEC option was enabled. iCOR is implemented as a plugin in the SNAP (Sentinel Application Platform) toolbox version 6.0. Once atmospherically corrected, all the used spectral bands were resampled to a 10-m. pixel size by the nearest neighbor resampling method in SNAP 6.0.
The Landsat 8 OLI product used was Landsat 8 Surface Reflectance (L8SR). The L8SR product [35] is processed with the LaSRC algorithm and made available by the United States Geological Survey (USGS). The images are radiometrically calibrated and orthorectified, projected in Universal Transverse Remote Sens. 2020, 12, 1514 7 of 33 Mercator (UTM) with the WGS84 datum. The L8SR product provides the surface reflectance image rescaled by a factor of 10,000. This product has been proven effective in aquatic applications [36][37][38]. For both satellite products, surface reflectance data were transformed into Rrs dividing by π.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 33 where is the predicted values, is the measured values and n is the number of samples.

Water Quality
The studied reservoir is an eutrophic waterbody with recurrent cyanobacterial blooms since 2012, documented through the Galician Platform for Environmental Information (GAIA) (www.gaia.xunta.es) ( Table 1). This periodic monitoring showed during 2017 Microcystis sp. as the dominant cyanobacteria, even though Chorococcus sp., Aphanocapsa sp. and Woronischia sp. were also present. The chl-a reached a maximum in the BC point (which is considered more representative of the main waterbody) of 21.26 µgr/l during May being the mean value for this year in BC sampling point 11.11 µgr/l. The highest cyanobacteria cell concentration was registered in October (9750 cel/ml) being Aphanocapsa sp. the dominant taxa. In the BP sampling point, the maximum chl-a was registered on 25th September (89.70 µgr/l), in this case, Microcystis sp. Being the dominant cyanobacterial taxa (115,500 cel/ml), and the mean annual value 20.96 µgr/l [28]. During 2018 (when the second UAV

UAV Imagery Pre-Processing
In order to increase the spatial and temporal resolution obtained with freely available satellite imagery, a complementary system was implemented based on the analysis of multispectral images taken with a commercial sensor (Rededge Micasense) (5 bands; Figure 2, Table 2) on board an octocopter (Atyges FV8). The first sampling took place on 19/09/2017and a second one, with the initial objective of model validation, was done during the same phenological period in 2018 (2th October).
The flights time was 12:00 am, established regarding the sun angle (between 30 • and 60 • to minimize sun-glint effects and specular reflections). Nine ground control points were distributed evenly in the flight area and x,y,z coordinates were surveyed with a GPS with submetric precision. The flight mission was planned as a simple grid with a distance between flight lines of 45 m. and a flight height of 120 m; the speed of the octocopter was set to 4 m/s.
No atmospheric correction was applied on the UAV imagery, as there were still not commercially available atmospheric corrections designed for aquatic purposes and UAV close-range remote sensing. Besides, the purpose of this work was to deliver an operational solution with no dependencies on Remote Sens. 2020, 12, 1514 8 of 33 concurrent field radiometric sampling, and therefore the effect of the thin atmosphere layer between the sensor and the water surface was neglected. The georeferencing, radiometric processing and mosaicking of the image were performed with the Pix4D software (Pix4D, Lausanne, Sw.). Raw multispectral images were corrected to ground absolute reflectance using the irradiance measures taken along with each snapshot by a specific sensor onboard the UAV and spectrally rescaled against a calibration panel with values around 70 % of diffuse reflectance in the visible-infrared spectrum. The extraction of the image information was done at different sampling sizes corresponding to the median of all the pixel values included in circular buffers of 0.5, 1.0 and 1.5 m. around each match-up point using QGIS 2.18.15. These buffers were set in order to test the coherence and representativeness of the data given by the sensor against a different size of sampling areas around the in situ sampling points.
In order to compare the reflectances recorded by the S2 MSI and UAV sensors, we used two simultaneous images of both sensors (10/02/2018). We made a focal analysis on the Micasense bands in which the extent of each of the S2 MSI pixels covering the reservoir was used as focal unit, once filtered those prone to spectral mixtures from the shore area. From these extents (squares of 10 × 10 m. labelled with an unique ID), a synthetic value of the Micasense reflectance for each band was calculated as the mean value of the Micasense image pixels covered by the area of each individual S2 MSI image pixel extent.

Data Analysis
For developing the empirical model, linear fits were applied expressing the relationships between Rrs at different wavelengths in every match-up point and their corresponding "in situ" chl-a data as a linear function of the format: where a and b are the model parameters (intercept and slope, respectively) and SBC is the tested Spectral Band Combination (Table 3). SBCs are either an index, a ratio or an algorithm. Some of these combinations (of two and three bands) are those proposed by different authors and adapted to OLI, MSI and Rededge spectral bands ( Figure 2, Table 2). As shown in Table 3, in addition to the standard 3 bands and 2 bands algorithms (3BDA and 2BDA) [21,22], we have included one modification of the 3 bands algorithm (3BDA_MOD) adapted to UAV-Rededge Micasense data.
These both algorithms (2BDA and 3BDA) come from a common initial conceptual model, relating remotely sensed reflectance and pigment content in vascular plant leaves [44]. The model was initially conceived to isolate the absorption coefficient of the pigment of interest from the reflectance spectrum. The 3BDA relationship used between the pigment content in the leaves and the reflectance was as follows: When the 3BDA model is tuned for the retrieval of chl-a in turbid productive waters, as described by [21,22,45], the spectral regions are defined as: • λ 1 : Spectral region such that R −1 λ1 is maximally sensitive to the absorption chl-a but still affected by the absorption of other constituents and backscattering: λ 1 = 660-690 nm. • λ 2 : Spectral region such that R −1 λ2 is minimally sensitive to chl-a for which the absorption by other constituents is almost equal to that at λ 1 : 710-λ 2 -730 nm. The 2BDA algorithm is described by the authors as a particular case of the 3BDA, which only uses the λ 1 and λ 3 terms of the equation and gives good results when a chla (λ 1 ) b b and a chla (λ 1 ) a tripton (λ 1 )+ a CDOM (λ 1 ) [39], any a being the absorption coefficients and b b the backscattering coefficient.
Two formulations for Sentinel 2 MSI data were tested for 2BDA and 3BDA, using Rrs from Bands 6 and 7 as R λ3 . Band 6 Rrs data for MSI have been mostly used in 3BDA and 2BDA algorithms formulated by other authors (e.g., [46]) as it offers more reliable data for water studies inside the NIR region of the spectrum than the rest of NIR bands in MSI. Yet, this band is centered exactly in 740 nm, which is lower than the theoretical wavelength specified by [23] who established the optimal wavelength above 740 nm and obtained the best results with 740-750 nm. Therefore, just one version of these two algorithms with Rrs data from Band 7 ( Figure 2, Tables 2 and 3) was tested.
The computation (R −1 λ 1 -R −1 λ 2 ) is related to the content of chl-a, but it is still affected by the variability in the backscattering of the medium; however, in our case, we also tried to apply this part of the algorithm alone with the UAV Multispectral Imagery data.
Other algorithms using red and NIR bands were tested, like the Normalized Difference Chlorophyll Index (NDCI) [50] and the Normalized Difference Vegetation Index (NDVI). The first one was inspired in NDVI and developed to retrieve chl-a emulating MERIS channels, thus using Rrs (708) and Rrs (665). NDVI was used in case of L8 OLI images in our study in absence of 708 nm band in this sensor, and also tested for Rededge Micasense.
Algorithms and indices using known features in blue and green spectral channels were also evaluated. We can divide them into two groups.
The first group used green and blue bands but also red and/or NIR. In this group, we have tested Surface Algal Bloom Index (SABI) [48] initially defined with MODIS bands to delineate the spatial distribution of floating micro-algal species in ocean waters with a NIR response similar to that of land vegetation, using in this case MODIS Rrs (869 nm) as NIR data. KIVU, developed for a freshwater lake with chl-a values below 6 µgr/L using Landsat 5 TM images [49], is also in this group. This algorithm uses blue and green bands (TM1 and TM2) and also the red band (TM3) to correct for the additional radiance caused by the scattering of non-organic suspended matter (see Table 3 for formulation and references).
In the last group, we tested Kab1 and Kab2, algorithms developed including blue and green bands and originally formulated for coastal areas using Landsat 7 ETM+ bands. Finally, in this group, we also tested simple ratios (GB1, GB2) and normalized indices (B3B1, B3B2) using only the green and blue bands available for each sensor.
Linear fits were adjusted both to chl-a data and Ln (chl-a) data, this approach also being tested by other authors to account for the non-normal distribution of the chl-a values [47,51].
The validation of the best performing models was done, only in case of satellite data, with images of the previous and next years (2016 and 2018). For both years, the validation was done with data in the chl-a range of 1-40 µgr/L for which the models were calibrated. A range of 4.38-20.42 µgr/L was used in case of Landsat 8 OLI, and 4.38-34.64 µgr/L in case of Sentinel 2 MSI.
where y p is the predicted values, y m is the measured values and n is the number of samples.

Water Quality
The studied reservoir is an eutrophic waterbody with recurrent cyanobacterial blooms since 2012, documented through the Galician Platform for Environmental Information (GAIA) (www.gaia.xunta.es) ( Table 1). This periodic monitoring showed during 2017 Microcystis sp. as the dominant cyanobacteria, even though Chorococcus sp., Aphanocapsa sp. and Woronischia sp. were also present. The chl-a reached a maximum in the BC point (which is considered more representative of the main waterbody) of 21.26 µgr/L during May being the mean value for this year in BC sampling point 11.11 µgr/L. The highest cyanobacteria cell concentration was registered in October (9750 cel/mL) being Aphanocapsa sp. the dominant taxa. In the BP sampling point, the maximum chl-a was registered on 25th September (89.70 µgr/L), in this case, Microcystis sp. Being the dominant cyanobacterial taxa (115,500 cel/mL), and the mean annual value 20.96 µgr/L [28]. During 2018 (when the second UAV flight was done) the general condition of the reservoir was very different, with bigger cyanobacterial blooms, the mean and maximum values of chl-a being higher than in 2017 (BCmean = 64.47 µgr/L, BPmean = 130.31 µgr/L). This documented highly changing condition in space and time makes this drinking water reservoir a suitable example of a sensible area which needs to be continuously monitored ( Table 1).
The water sampling (5 sampling points) made synchronous with the first UAV flight (19th September 2017) showed a situation of low nutrient concentration and phytoplankton biomass (low both chl-a and PC), even though the Secchi Depth values together with turbidity and TSS indicate a slightly turbid state (Table 4). Even though the phenological phase was considered the same in 2018 and in 2017, the second UAV flight was done in very different water quality conditions with a cyanobacterial bloom in the reservoir and much higher chl-a and PC concentrations (Table 4). Curiously, turbidity was higher in the first flight. Even though the same number of samples (5) was taken in either campaign, we have 4 data in 2018 as one of the chl-a samples had to be discarded.

Spectral Band Combinations for the Retrieval of chl-a
The empirical approach was selected to build models for the retrieving of chl-a from the three types of multispectral images. All the indices and band ratios tested were adapted to the characteristics of the bands of each sensor (Tables 2 and 3, Figure 2). In case of satellite imagery, the time window used for the development of the algorithms (3 days) can be a source of error in case of sudden changes in the development of an algal bloom. The chosen year for the calibration of the models was the one with more potential matchups (more in situ data) and stable conditions regarding chl-a concentration (cf. Table 1), in order to reduce the uncertainty in the calibration of the models regarding this time difference. To evaluate this effect, we plotted the scatter plot graphics for the calibrated models for retrieving chl-a with a color code reflecting this time difference (0-3 days) (e.g., Figure 4F).

Landsat 8 Imagery
In the calibration to find a suitable linear model for the determination of chl-a (Equation (1)) with Landsat 8 OLI multispectral data, most of the algorithms and band indices tested were positive and significatively correlated (Pearson correlation) with chl-a concentration ( Table 5). The linear fit was made with both chl-a data inside an interval of 5.65-18.09 µgr/L, and Ln transformed chl-a data, giving the latter the best results.
Given the inherent complexity of inland waters, the combinations using red and NIR bands, specifically 2BDA, SABI and NDVI, were the best performing Spectral Band Combinations. The Surface Algal Bloom Index (SABI) which uses green and blue bands in addition to red and NIR [48] also gave good results. These three had the strongest (Pearson r > 0.8) and the more significative correlations with Ln(chl-a) data (Table 5, Figure 4). The results obtained with un-transformed Ln data did not show high differences, the Pearson r and R 2 being 0.858 and 0.753 for 2BDA, 0.844 and 0.712 for SABI and 0.834 and 0.6953 for NDVI models, respectively. Kab_1 was also tested as its results were considered good, but the scatterplot showed a relationship based in two defined and separated group of points.
Chl-a models fitted for the 2017 data were subsequently applied to Landsat 8 OLI datasets from 2016 and 2018. The validation, with chl-a from 4.38 to 20.42 µgr/L, showed very similar results for the three best performing indices ( Figure 5, Table 6). However, all showed a tendency to overestimate chl-a values (bias ranging from 0.94 to 1.13). The lowest root mean square error (RMSE) was obtained with the NDVI index, which also gave the lowest NRMSE and MAPE, very similar to 2BDA, but the highest bias. Despite the good results in the calibration, the validation did not show a good agreement between measured and predicted values. Even though Kab_1 validation metrics were better that the other models, the graphical output, both for the calibration and the validation, shows a weaker performance ( Figure 5, Table 6).
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 33 Chl-a models fitted for the 2017 data were subsequently applied to Landsat 8 OLI datasets from 2016 and 2018. The validation, with chl-a from 4.38 to 20.42 µgr/l, showed very similar results for the three best performing indices ( Figure 5, Table 6). However, all showed a tendency to overestimate chl-a values (bias ranging from 0.94 to 1.13). The lowest root mean square error (RMSE) was obtained with the NDVI index, which also gave the lowest NRMSE and MAPE, very similar to 2BDA, but the highest bias. Despite the good results in the calibration, the validation did not show a good agreement between measured and predicted values. Even though Kab_1 validation metrics were better that the other models, the graphical output, both for the calibration and the validation, shows a weaker performance ( Figure 5, Table 6).  In the calibration of linear models with Sentinel 2 MSI data, the number of significative correlations with the Spectral Band Combinations was less than for OLI data; the relationships were also slightly weaker. Still, 5 of them were strong, with a Pearson correlation coefficient above 0.7 and

Sentinel 2 Imagery
In the calibration of linear models with Sentinel 2 MSI data, the number of significative correlations with the Spectral Band Combinations was less than for OLI data; the relationships were also slightly weaker. Still, 5 of them were strong, with a Pearson correlation coefficient above 0.7 and the highest significance level. In this case, the best fits were obtained with the chl-a value, and not with the Ln transformed data. The chl-a values in this case ranged from 5.09 to 50.37 µgr/L.
As with OLI data, all the linear models with Pearson r > 0.8 and a p value lower than 0.0001 were considered the best performing and were tested in the validation stage. In this case, we also included in the validation the 3BDA algorithm, as it also gave good results and the correlation coefficient was close to 0.8 (Table 6). In case of Sentinel 2 data, not only did the combination of the NIR and red bands fit best with chl-a values but also 3 Spectral Band Combinations with green and blue bands gave similarly good results (Table 7, Figure 6).
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 33 the highest significance level. In this case, the best fits were obtained with the chl-a value, and not with the Ln transformed data. The chl-a values in this case ranged from 5.09 to 50.37 µgr/l. As with OLI data, all the linear models with Pearson r > 0.8 and a p value lower than 0.0001 were considered the best performing and were tested in the validation stage. In this case, we also included in the validation the 3BDA algorithm, as it also gave good results and the correlation coefficient was close to 0.8 (Table 6). In case of Sentinel 2 data, not only did the combination of the NIR and red bands fit best with chl-a values but also 3 Spectral Band Combinations with green and blue bands gave similarly good results (Table 7, Figure 6).
The results obtained in the linear fits with Ln transformed data showed a lower performance in case of Sentinel 2, the Pearson r and R 2 being 0.718 and 0.5167 for 2BDA_2, 0.660 and 0.436 for 3BDA_2, 0.653 and 0.427 for Kab_1 and 0.659 and 0.435 for B3B2, respectively.  The results obtained in the linear fits with Ln transformed data showed a lower performance in case of Sentinel 2, the Pearson r and R 2 being 0.718 and 0.5167 for 2BDA_2, 0.660 and 0.436 for 3BDA_2, 0.653 and 0.427 for Kab_1 and 0.659 and 0.435 for B3B2, respectively.
As can be seen in Figure 6, the models using blue and green bands (Kab_1, B3B2) gave relatively high R 2 values but the correlation was highly dependent on one data point, so they were not considered robust. The results of the calibration/validation exercise are given (Table 8, Figure 7). As can be seen in Figure 6, the models using blue and green bands (Kab_1, B3B2) gave relatively high R 2 values but the correlation was highly dependent on one data point, so they were not considered robust. The results of the calibration/validation exercise are given (Table 8, Figure 7).
The validation of these models was done, as with OLI data, with the valid images for the previous and following years (2016 and 2018) and a chl-a interval from 4.38 to 34.64 µgr/l. The validation metrics showed different performances for these 5 models (Table 8). 2BDA tended to underestimate the chl-a values, giving even negative results (Figure 7), while 3BDA showed a positive bias of 0.25. All indices with good performance in calibration using blue and green bands (B3B2, Kab_1 and GB2) tended to greatly overestimate the chl-a concentration. The best model regarding the validation results for Sentinel 2 MSI was the one using 3BDA algorithm and Band 7 Rrs as the NIR data input, with the lowest values for all metrics, followed by 2BDA, as can also be seen in the validation scatterplots ( Figure 7)  The validation of these models was done, as with OLI data, with the valid images for the previous and following years (2016 and 2018) and a chl-a interval from 4.38 to 34.64 µgr/L. The validation metrics showed different performances for these 5 models (Table 8). 2BDA tended to underestimate the chl-a values, giving even negative results (Figure 7), while 3BDA showed a positive bias of 0.25. All indices with good performance in calibration using blue and green bands (B3B2, Kab_1 and GB2) tended to greatly overestimate the chl-a concentration. The best model regarding the validation results for Sentinel 2 MSI was the one using 3BDA algorithm and Band 7 Rrs as the NIR data input, with the lowest values for all metrics, followed by 2BDA, as can also be seen in the validation scatterplots ( Figure 7) A clear pattern regarding the influence of the time window in the calibration of the algorithms (Figures 4 and 6) could not be found in the graphical analysis for Sentinel 2 but Landsat 8 graphics showed the worst fit for the 3-day time-lapse points. Overall, once the validation was done, the best performing model for retrieving chl-a in the studied reservoir was the one fitted with the 3BDA algorithm applied to Sentinel 2 imagery. The validation results for Landsat 8 were poorer than Sentinel 2.

UAV Imagery. Model Calibration
The two UAV flights (2017 and 2018) were made in late summer-early autumn, when the cyanobacterial bloom events take place in this area. This sampling plan was intended to calibrate and subsequently validate the models; but two opposite conditions of algal bloom development were found (extremely different chl-a values, cf. Table 4). The chl-a interval covered was too broad to validate the models obtained with the 2017 flight with the data of the 2018 flight, so the validation was not possible. Nevertheless, these two opposite situations gave the opportunity to test the multispectral sensor performance in low and high chl-a concentrations, test the coherence of the relationships between flights and to test which of the Spectral Band Combinations were most "portable" to monitor chl-a concentration in waterbodies with different eutrophication levels.
The 2017 flight, done in low chl-a concentration conditions, showed a better performance of the Spectral Band Combinations using NIR and red bands. Even though some positive correlations were found, the only significative one applies a modified version of the 3BDA model (3BDA_ MOD ), for all the three sampling areas tested. (Tables 9-11, Figure 8).
The second flight results, done in high chl-a concentration conditions, confirmed the combination of red and NIR bands as the most suitable for chl-a retrieval with the Rededge Micasense sensor. In this case, the best performing models were always 2BDA and SABI (Figure 8, Tables 12-14).
In another exercise, the data of both flights were analyzed together. We used the 9 available points (5 from 2017 and 4 for 2018) to calibrate linear models that could be adjusted to this broad range of chl-a concentrations. As it was assessed that the results were similar for the 3 buffers tested (Tables 9-14), this exercise was done only with the 0.5-m. buffer data. Even though this uneven distribution of data is associated with a high level of uncertainty, we considered the obtained results promising and worthy of being shown. Table 9. Linear model coefficients and linear fit parameters for all the Spectral Band Combinations tested for the retrieval of chl-a with multispectral data from Rededge Micasense on board the UAV platform for the 2017 flight. The models were done with in situ chl-a and the median of the reflectance values included in a buffer of 0.5 m. around the sampling point.

Index
Interc     The combination of these bands (blue-green and green-red) gave the results shown in Table 15 and Figure 9. It can be seen that this relationship is based upon two separate groups of points, so it is not considered good for directly retrieving absolute chl-a values accurately; however, the results suggested that it can be used to classify the pixels into broad chl-a categories or classes to which apply algorithms specific to this narrower chl-a range inside of each class. This is also a key issue when dealing with cyanobacterial blooms with a spatial patchy distribution in waterbodies, in order to spatially stratify the analyses.    chl-a condition, which corresponds with 2018 flight data. The results correspond to the calculations made with data included in a 0.5, 1.0 and 1.5-m. buffers, which are shown following this order from top to base.
In another exercise, the data of both flights were analyzed together. We used the 9 available points (5 from 2017 and 4 for 2018) to calibrate linear models that could be adjusted to this broad range of chl-a concentrations. As it was assessed that the results were similar for the 3 buffers tested (Tables 9-14), this exercise was done only with the 0.5-m. buffer data. Even though this uneven distribution of data is associated with a high level of uncertainty, we considered the obtained results promising and worthy of being shown.
The only band combinations that fitted in a linear model applied to this broad range of chl-a concentration were those using the first part of the visible spectrum (blue, green and red bands). The combination of these bands (blue-green and green-red) gave the results shown in Table 15 and Figure  9. It can be seen that this relationship is based upon two separate groups of points, so it is not considered good for directly retrieving absolute chl-a values accurately; however, the results suggested that it can be used to classify the pixels into broad chl-a categories or classes to which apply algorithms specific to this narrower chl-a range inside of each class. This is also a key issue when dealing with cyanobacterial blooms with a spatial patchy distribution in waterbodies, in order to spatially stratify the analyses.

Index
Interc   As it is shown in Table 15 and Figure 9, the best results were obtained for a normalized index using blue and green bands (B3B1), which gave the highest correlation coefficient (Pearson r = 0.9907) and R 2 (0.9816) and the highest significance level.
Chl-a = (520.310 * B3B1) − 41.980 (8) Remote Sens. 2020, 12, 1514 20 of 33 With these first results, we designed a preliminary ensemble model that uses the B3B1 algorithm as the first to be applied to the obtained images to classify the pixels into broad ranges of chl-a content. In our case, we could establish only two classes, dividing them in the middle of our 2 datasets: Class 1: chl-a < 50 µgr/L and Class 2: chl-a > 50 µgr/L. Once the chl-a value for each pixel in the images is obtained using this algorithm and classified into these 2 classes, the calibrated algorithms for low and high chl-a is subsequently applied to each pixel. The results of the application of B3B1 algorithm are shown in Figure 10, and the subsequent application for low and high chl-a value algorithms is shown in Section 3.3.
Remote Sens. 2020, 12, x FOR PEER REVIEW 20 of 33 As it is shown in Table 15 and Figure 9, the best results were obtained for a normalized index using blue and green bands (B3B1), which gave the highest correlation coefficient (Pearson r = 0.9907) and R 2 (0.9816) and the highest significance level.
Chl-a = (520.310 * B3B1) -41.980 (8) With these first results, we designed a preliminary ensemble model that uses the B3B1 algorithm as the first to be applied to the obtained images to classify the pixels into broad ranges of chl-a content. In our case, we could establish only two classes, dividing them in the middle of our 2 datasets: Class 1: chl-a < 50 μgr/l and Class 2: chl-a > 50 μgr/l. Once the chl-a value for each pixel in the images is obtained using this algorithm and classified into these 2 classes, the calibrated algorithms for low and high chl-a is subsequently applied to each pixel. The results of the application of B3B1 algorithm are shown in Figure 10, and the subsequent application for low and high chl-a value algorithms is shown in section 3.3. As it is shown in Figure 10, the image of 2017 was classified by the B3B1 algorithm in a low chla value, with almost all the pixels showing chl-a values lower or equal than 50 μgr/l chl-a unless in the northern shore; and most of them below 20 μgr/l chl-a. In this case, the 3BDA_MOD algorithm was subsequently applied to most of the image for a better classification and the SABI algorithm to the pixels identified as > 50 μgr/l chl-a. The image of 2018 was classified by the algorithm into the high chl-a range, with just residual pixels in the south border included in the lower class. In this case, the SABI algorithm for high chl-a concentration was applied to almost all the image. (cf. section 3. 3) The future of this preliminary ensemble model is to validate the use of green and blue bands for the initial classification with intermediate chl-a data and develop a classification algorithm to discriminate intermediate chl-a classes with these two bands. Moreover, specific algorithms have to be developed for intermediate chl-a ranges, as this range was not covered by the data obtained in our study.

3.3.Combined Monitoring Tool
The initial design of the tool contemplates continuous baseline monitoring of chl-a concentration with satellite imagery, and the use of UAV imagery when there is a lack of satellite data due to maintained cloud cover in the area, or a threatening event in the waterbody, so the algorithms found for the UAV can operate as a "stand alone" part of the approach. Nevertheless, if these data can be made fully comparable to satellite data, the integration of UAV images inside of the chain will enhance the overall performance. To account for this, a comparison between the results of both data has been made with the only simultaneous image available for satellite (S2 MSI) and UAV (10/02/2018). The Rrs for each pixel of S2 MSI in bands 2, 3, 4, 5 and 8 were compared to the mean Rrs value of all the pixels of the UAV image inside the same extent of 10 × 10 m., using, in this case, all As it is shown in Figure 10, the image of 2017 was classified by the B3B1 algorithm in a low chl-a value, with almost all the pixels showing chl-a values lower or equal than 50 µgr/L chl-a unless in the northern shore; and most of them below 20 µgr/L chl-a. In this case, the 3BDA_MOD algorithm was subsequently applied to most of the image for a better classification and the SABI algorithm to the pixels identified as > 50 µgr/L chl-a. The image of 2018 was classified by the algorithm into the high chl-a range, with just residual pixels in the south border included in the lower class. In this case, the SABI algorithm for high chl-a concentration was applied to almost all the image. (cf. Section 3. 3) The future of this preliminary ensemble model is to validate the use of green and blue bands for the initial classification with intermediate chl-a data and develop a classification algorithm to discriminate intermediate chl-a classes with these two bands. Moreover, specific algorithms have to be developed for intermediate chl-a ranges, as this range was not covered by the data obtained in our study.

Combined Monitoring Tool
The initial design of the tool contemplates continuous baseline monitoring of chl-a concentration with satellite imagery, and the use of UAV imagery when there is a lack of satellite data due to maintained cloud cover in the area, or a threatening event in the waterbody, so the algorithms found for the UAV can operate as a "stand alone" part of the approach. Nevertheless, if these data can be made fully comparable to satellite data, the integration of UAV images inside of the chain will enhance the overall performance. To account for this, a comparison between the results of both data has been made with the only simultaneous image available for satellite (S2 MSI) and UAV (10/02/2018). The Rrs for each pixel of S2 MSI in bands 2, 3, 4, 5 and 8 were compared to the mean Rrs value of all the pixels of the UAV image inside the same extent of 10 × 10 m., using, in this case, all available bands of the Micasense sensor (1 to 5). The results of this comparison show a good agreement in the shape of the spectral signature for both sensors but a drift in the absolute value of Rrs was found ( Figure 11). These differences are evaluated through RMSE between analogous bands for both sensors as 0.0051, 0.0044, 0.0049, 0.0040 and 0.0046 sr −1 for blue, green, red, rededge and NIR bands, respectively ( Figure 11).
Remote Sens. 2020, 12, x FOR PEER REVIEW 21 of 33 available bands of the Micasense sensor (1 to 5). The results of this comparison show a good agreement in the shape of the spectral signature for both sensors but a drift in the absolute value of Rrs was found ( Figure 11). These differences are evaluated through RMSE between analogous bands for both sensors as 0.0051, 0.0044, 0.0049, 0.0040 and 0.0046 sr -1 for blue, green, red, rededge and NIR bands, respectively ( Figure 11). This comparison was made for the entire reservoir and dividing the pixels into central (45 central pixels) and outer pixels, to extract information about the adjacency effects in the Rrs values. As can be seen in Figure 11, both sensors registered high Rrs in the NIR band when all data are plotted together, and also when only the outer pixels are analyzed, but this effect disappears when only the central pixels data are used. This comparison was made for the entire reservoir and dividing the pixels into central (45 central pixels) and outer pixels, to extract information about the adjacency effects in the Rrs values. As can be seen in Figure 11, both sensors registered high Rrs in the NIR band when all data are plotted together, and also when only the outer pixels are analyzed, but this effect disappears when only the central pixels data are used.
The coherence between the shapes of the curves lead to think that the relationships between bands should be maintained among sensors and indices (specially normalized indices), leading to a potential integration between these sensors. We tried all the indices in Table 2 using these 5 bands, and none of them gave a significative relationship. More data in different trophic conditions are needed to find the proper relationships among these two sensors.
In its current form, the workflow includes the application of the best performing algorithms to valid images from both Landsat 8 and Sentinel 2 satellites (Table 16) to generate chl-a thematic maps on a regular basis. In case of UAV data, the algorithms for low and high chl-a values will be applied once the pixels in the image are classified into these classes by the B3B1 algorithm. The algorithms applied were those developed with the results of the 0.5 buffer sampling area: Table 16. Selected algorithms applied to the images of the reservoir obtained with Landsat 8 OLI, Sentinel 2 MSI and Rededge Micasense on board UAV. The combined monitoring approach gave the graphical results shown in Figure 12. This figure shows the graphical output for September and October of 2017 and 2018, considered the most critical period for cyanobacterial bloom development in this area.

Satellite -Sensor Algorithm
During 2017, there is a good agreement between Landsat 8, Sentinel 2 and in situ water samples (Figure 12), apart from the second half of September. In this period, when the UAV flight was done, the simultaneous in situ water sampling made with the flight did not agree either with the in situ periodic monitoring data. While we found chl-a concentrations ranging from 1.34 to 2.68 µgr/L, the data from the Galician Environmental Information System showed concentrations of 14.09 and 13.51 µgr/L one day earlier, which means a sharp decrease in chl-a concentration in the reservoir during that time. The last figure of September 2017 shows how the point sampling in the area of wind driven concentration of cyanobacterial biomass in control point BP found 89.7 µgr/L of chl-a (dominant cyanobacteria: Mycrocistis sp.), which was not detected by the algorithm applied to Landsat 8 imagery that showed an homogenous chl-a concentration in the reservoir in agreement with the in situ chl-a value found in the second control point (BC).
For the following year (2018), the results showed that the used algorithms worked well in the lower chl-a concentrations, as was expected due to the chl-a concentration values available for the calibration. Nevertheless, in the last figure, corresponding with 10/24/2018 images, we could also see good agreement in high chl-a concentration values between satellites and in situ data, better with Sentinel 2 imagery, with which we could even detect higher chl-a values (around 50 µgr/L). The comparison with the UAV flight of 2018 and the official monitoring in situ data (taken one day earlier) shows good agreement in sampling station BP but not in BC which can be a point circumstance in this station, as, in our sampling of the following day, the in situ data showed high chl-a values in all the sampling points.
Overall, the developed Landsat 8 algorithms worked well in low chl-a values but could not detect high values concentrated in one small area. Given the results obtained in the validation, priority should be given to the results obtained with Sentinel 2 images, which showed an overall better agreement with the ground truth data. Landsat 8 data should only be used to detect general trends in the reservoir in absence of other data source. UAV images can be used in absence of valid satellite images to complement the in situ data with spatially explicit information. During 2017, there is a good agreement between Landsat 8, Sentinel 2 and in situ water samples (Figure 12), apart from the second half of September. In this period, when the UAV flight was done, the simultaneous in situ water sampling made with the flight did not agree either with the in situ periodic monitoring data. While we found chl-a concentrations ranging from 1.34 to 2.68 µgr/l, the data from the Galician Environmental Information System showed concentrations of 14.09 and 13.51 µgr/l one day earlier, which means a sharp decrease in chl-a concentration in the reservoir during

Discussion
Phytoplankton blooms can be highly variable in space and time, cyanobacterial blooms being especially complex due to their buoyancy capacity that enables them to move vertically along the water column [52]. These blooms present a spatial patchy distribution and a fast development of even only some hours [53]. Therefore, their detection and monitoring are difficult tasks that cannot be tackled only by in situ point measurements; the synoptic view offered by satellites and UAV multispectral or hyperspectral imagery is an excellent complement to this routine type of sampling [14,54] The data provided by Sentinel 2 and Landsat 8, when combined in the form of a virtual constellation, opened up new opportunities for capturing limnological dynamics that reflect rates of change at an unprecedent resolution [26,36,54]. The Landsat 8 equatorial repeat cycle (since 2013) is 16 days, while each of the Sentinels 2 is 10 days, which becomes 5 days when the two are combined (Sentinel 2A since 2015 and S2B since 2017). The revisit time of the virtual constellation formed by the Landsat 8 and Sentinel 2 satellites has been estimated at~2.9 days [25], which can be even lower in waterbodies located in latitudes where two overpasses of any of the satellites overlap. Not only the temporal but also the spatial resolution of these satellites (10 m. in case of Sentinel 2 MSI and 30 m. in case of Landsat 8 OLI) has broadened the scope of the inland water quality monitoring through satellite imagery, making it possible to capture much smaller waterbodies than the ones usually studied with dedicated multispectral ocean color sensors on board other satellites (e.g., MERIS or Sentinel 3 OLCI with a 300-m. resolution).
One way to make the most of this combination of platforms and sensors is to validate which models or combination of models are suitable for each type of sensor and type of waterbody and combine the results into a single multi-platform tool with medium to high spatial resolution (10-30 m), which will raise to a very high resolution (8 cm.) in the events of UAV imagery availability. As a pilot experience, we studied a very small waterbody, focusing the efforts on the retrieval of chl-a concentration as a proxy for phytoplankton biomass and trophic status.
The inference process followed in Remote Sensing to obtain information about objects from measurements made from space is always less than perfect. The development of remote sensing products to quantify biophysical parameters is subject to different types and degree of uncertainties due to different sources of error inside of the processing chain. Based on the Image Chain Approach to identify uncertainty in Remote Sensing [55,56], every step in the processing chain from the capture of the image by the sensor to the final product has an associated uncertainty. This uncertainties propagate along the entire chain, but we will focus the analysis of the uncertainties, in our proposed processing scheme, in the weakest points identified, which also have the potential to be improved in the framework of any established monitoring program: the selection of the time window for the match-up points and the performance of the applied atmospheric correction.
The first source of uncertainty is related to the time lapse (± 3 days) between the image acquisition and the in situ ground truth data. Loew et al. [57], in their description of uncertainties in the validation process, state that, regarding the uncertainty due to imperfect spatiotemporal collocation, a compromise should be achieved among two concerns: on the one hand, the measurements should be close to each other relative to the spatiotemporal scale on which the variability of the field becomes comparable to the measurement uncertainties; and on the other hand, the need for sufficient collocated pairs to allow a statistical analysis. In our case, the compromise to minimize the uncertainty of the obtained product had to be achieved among the available data (potential matchups) and the frequency and duration of the algal blooms in the area. Trying to minimize this uncertainty, the data of 2017 were used for algorithm calibration, when the intensity and duration of the algal blooms showed a lower change rate of the 3 studied years (Table 1). Moreover, the extremely high values in sampling point BP were not used in the calibration and validation exercises, being considered exceptional both in time and space and thus with the higher probabilities of not being captured by the satellite images. However, the valid time difference between in situ sampling and satellite overpass used in our work, similar to that used by [51,58,59], is a source of uncertainty that was taken into account when interpreting the results. The uncertainty regarding spatio-temporal collocation can be easily reduced by achieving a compromise with the official agencies in charge of in situ monitoring in which a prioritization of satellite overpassing dates should be established for in situ sampling.
The second source of higher uncertainty identified was the atmospheric correction algorithm applied for each source of data. Lacking vicarious calibration data simultaneous with the satellite overpass to validate the applied AC, we could not give a quantitative estimate of our uncertainty. Nevertheless, this exercise has been done for the Sentinel 2A MSI sensor by Warren et al. [60] in a thorough study evaluating six free source processors for the atmospheric correction over coastal and inland waters. In their comparison against above-water optical in situ measurements, none of the processors performed well over the entire band set tested (444-865 nm). The authors quantify the uncertainties related to the AC algorithms through the value of the Mean Absolute Percentage Difference (ψ), finding the that green (560 nm) and blue (490 nm) bands were best handled by all ACs, but the red, red-edge and NIR bands had often very high uncertainties. More precisely, the uncertainties associated with the atmospheric correction performed by iCOR in Bands 4, 5 and 7 (the most useful in our study to compute 2BDA and 3BDA algorithms) were quantified ψ = 106.05%, ψ = 149.00% and ψ = 681.75%, respectively, while the Root Mean Squared Errors (RMSE) were 0.0029, 0.0027 and 0.0033 sr −1 , respectively.
The performance of atmospheric correction applied to MSI data in our study was also evaluated over inland waters by Pereira-Sandoval et al. [61], showing that other processors were more effective in case of oligotrophic to meso-eutrophic waters, while iCOR showed moderately good results in hypereutrophic environments.
In case of Landsat 8, a similar exercise was performed by Ilori et al. [62] in which four AC processors' performances (including the Landsat 8 Surface Reflectance Code (LaSRC)) were tested against in situ data acquired mainly over coastal waters from the OC component of the Aerosol Robotic Network (AERONET-OC) for Bands 1 to 4 of OLI sensor. Even though AERONET-OC sites are not representative of inland waters, their results showed a poor performance of L8SR product in OLI Bands 1 and 2 but a good performance at a few sites located in nearshore and inland waters. In this assessment, Bands 3 and 4 (used in our final chosen algorithms) showed an RMSE of 0.0030 and 0.0022 sr −1 , respectively, which is similar to what was found for Sentinel 2 MSI by Warren et al. [60] for Sentinel 2 MSI. Even though other valuable comparative exercises have been done regarding the performance and associated uncertainties of the AC processor for OLI data used in our study [63], they were performed over different land covers and so were not taken as a reference for our work.
Another source of uncertainty regarding AC is the adjacency effect related to the distance to shore. This effect cause that the reflected light field from the neighboring land pixels, which are brighter, gets into the path from the target (water pixel) to the sensor through atmospheric scattering [64]. This also makes the reflectance value from the water pixels near the coast larger, and, hence, these datasets are prone to errors for the derivation of water quality parameters. This effect has been shown in the analysis made in our study comparing the Rrs reflectance of the outer and central pixels of the reservoir. Our results show how the Rrs reflectance of outer pixels have a higher value in the NIR band (B8 S2 MSI) than the central ones. This high NIR reflectance over inland and coastal waters due to adjacency effects has been already reported based on airborne and spaceborne data [8]. In the study of Warren et al. [60], it was also shown that iCOR performed poorly in coastal waters but much better in inland waters, due to the land/based methodologies applied to this processor. This result supports the use of iCOR when adjacency effects can be challenging, as in our case.
In spite of this, not only in our study but in some others using different AC processors [59,61] and methodologies [46], promising results were obtained in retrieving chl-a using red and NIR band combinations from Sentinel 2 MSI data. Further developments in the correction of these bands and the specific correction of the adjacency effects, which is only built-in in iCOR of all the publicly available AC processors, it is expected to improve the retrieval of chl-a in inland waters and reduce the associated uncertainties.
The propagation of these errors until the product output is leading to the uncertainties found in our product validation phase. The uncertainty associated with the final output of our Remote Sensing chain was evaluated for every model with different pairwise metrics describing bias and accuracy, as is discussed below.
The development of algorithms for Landsat 8 OLI, with a spatial resolution of 30 m., is hampered by the small size of the reservoir, that makes the adjacency effects an important issue in this case [14]. Moreover, the need to move the match-up points in the reservoir from the shoreline in order to avoid this effect as much as possible can introduce a source of error in the algorithm development. Nevertheless, almost all the algorithms and Spectral Band Combinations tested were positively and significatively correlated with in situ chl-a values. Both the calibration and validation showed that the most appropriate Spectral Band Combinations in Landsat 8 OLI for the retrieval of chl-a were those based on red and NIR bands, specifically in our study: 2BDA, NDVI and SABI. The SABI is an empirical algorithm developed in order to detect water floating biomass that has an NIR response similar to land vegetation, and also includes blue and green spectral bands [48].
Although the correlations in the calibration of the models were good, our validation results did not show a good match between the measured and the predicted values ( Figure 5). Other studies did not obtain good results with OLI data when developing algorithms applying 2BDA [38,51,58]. Watanabe et al. [38] found no relationship with simulated Rrs spectra for OLI/Landsat-8 based on field data in the eutrophic Barra Bonita Reservoir (Brazil), a bigger and deeper waterbody. In this work, also made in a higher chl-a concentration circumstance (mean values ranging from 120.4 to 428.7 µgr/L), the authors suggest that the 2BDA index, as proposed for clear waters, could not explain the chl-a distribution in their reservoir. The best results were obtained in this case with the SLOPE model developed by Mishra and Mishra [65], which uses red and green bands. Similarly, Boucher et al. [51] did not find a suitable model using the 2BDA algorithm with Landsat 8 images Level 1 corrected atmospherically with the dark object subtraction (DOS) method. In this study, made on a regional scale, including different lakes in New Hampshire and Maine (US), the best performing models where Kab1, Kab2 and KIVU, all based in blue and green bands. None of them used the L8SR product for the calibration of the models, which is also a methodological difference. This product has been evaluated for aquatic applications by Bernardo et al. [37]. L8SR was considered by the authors to provide suitable Rrs values for all OLI bands despite achieving the poorest results in Band 5 (NIR). However, our model calibration results with red and NIR bands show that this level 2 product made available by the USGS could be applied successfully for chl-a retrieval even in small waterbodies and should be further investigated with different trophic situations. Not only Red-NIR combinations were found correlated with in situ data in our study, but also algorithms using blue and green bands, as SABI, with very similar values in the validation of statistical metrics, but worst results in the validation graphical output. The worst results were obtained with all algorithms and indices using Band 1 (coastal-aerosol) (Kab2, B3B1 and GB1; cf. Table 3, Table 5 and Figures 4 and 5). Our results confirm that green and blue band combinations are not the most suitable for retrieving chl-a in complex waters but also that the performance improves when Band 2 is used. This shows the difficulty of performing an accurate atmospheric correction for this ultra-blue band for inland aquatic applications. Wang et al. [66] also found big errors when computing band ratios with Band 1 data atmospherically corrected with different AC algorithms in their assessment of Landsat-8 OLI atmospheric correction algorithms for inland waters, even though they did not assess the L8SR product. The authors reported the same situation for the band ratios formed by the NIR band with all AC algorithms assessed. Sentinel 2 data have been extensively tested for chl-a retrieval since its launching in 2015 and 2017, both applying empirical and analytical modelling approaches [14,46,54,59,61,67]. Our results agree with the studies applying empirical modelling approaches [46,59,61], which evaluated the performance of the 3BDA model with good results, even though all of them used what we called here the 3BDA_1 version (with B6 as λ 3 ). Moreover, in agreement with our study are the results obtained by Xu et al. (2018), who also found good relationships with 3BDA, 2BDA and NDCI, although they used a different formulation of the 2BDA algorithm (B5/B4).
To the best of our knowledge, there is only one other study showing the performance of the commercial sensor Rededge Micasense for chl-a retrieval in waterbodies [68], but, in this case, the authors do not use Spectral band Combinations based in known spectral features of the chl-a but fit single and multiple variable linear approaches concluding that the best models in their case where those using green and red bands. We can nevertheless consider the results of our study as the first to test the algorithms developed, be proven effective for satellite remote sensing, adapted to Rededge Micasense sensor, and be tested in different trophic conditions. Again, the lack of field radiometry (as also happens in the study of Arango and Nairn [68]) entails a higher level of uncertainty in the final product. A radiometric calibration made by Cao et al. [69] for Rededge Micasense, showed in laboratory experiments a good relationship between ASD reflectance and Micasense with R 2 all above 0.96 for all bands. In field experiments, the comparison with the Pix4D processed image showed lower accuracy (overestimation) in the rededege and NIR regions, but a consistency in the spectral signature shape which reduced the uncertainty when using normalized indices.
When analyzing separately each flight data, the most suitable algorithms were again different combinations of red and NIR bands, including a modification of the 3BDA which excludes the last part of the equation. In our case, using the first computation involving λ 1 and λ 2 , without removing the effect of the backscattering of the medium with the use of λ 3 , gave the best results with this sensor in low chl-a conditions. The optimal position of λ 3 in 2BDA and 3BDA algorithms were established in practice by Dall'Olmo and Gitelson [70] in 720-740 nm in turbid productive waters. The NIR band of the Rededege Micasense is centered in 840 nm, hence too far from the optimal range, and proved not useful in our study in low chl-a values. As expected by the formulation of the algorithm, when applied in turbid conditions with high chl-a during the second UAV flight, this modification did not perform well (R 2 = 0.4; Pearson r = 0.6) as the effect of the backscattering of a medium with a high load of suspended particles should be removed with the λ 3 spectral region. Therefore, in this case 2BDA and SABI were the best performing algorithm in the presence of a bloom and high chl-a concentration. Again, the NIR spectral band in Rededge Micasense, used in this version of the 2BDA model in coherence with the formulation initially used for 3BDA, is far away from the optimal 720-740 nm range established by Dall'Olmo and Gitelson [70] to be used turbid productive waters, but performed better than the version of this algorithm made with the rededge band (2BDA_2). This band is centered in 717 nm, nearer of the optimum range, but the authors describe the aim of this band λ 3 as to be in a spectral region where the absorption by pigments, tripton and CDOM is negligible, which is something that happens at λ > 740 [21] so the use of longer wavelengths (840 nm in our case) might be feasible. Moreover, as chl-a increases, the λ 3 spectral region of maximal sensitivity to chl-a shifts toward longer wavelengths, and the precise optimal position of λ 1 and λ 3 will depend on the relative importance of the interferences and on the trophic status of the waterbody [70].
When analyzing the data of both flights in a combined dataset, NIR and red bands got lower correlations than the combination of blue-green and red-green. These last combinations were the best to discriminate the differences between very different chl-a concentrations, being this also reported by Arango and Nairn [68] with the same sensor. They also developed the models with data of two opposite trophic states and mean chl-a concentrations (oligotrophic: 8.52 µgr/L and eutrophic: 352.60 µgr/L) finding the green and red bands as the most significant.
This finding was used to develop a preliminary ensemble model for the retrieval of chl-a with UAV-Rededege data. In the ensemble model the pixels of the images were first classified in two chl-a classes, and subsequently a second model calibrated for that chl-a class was applied to every pixel to obtain a more accurate result. This model is still to be validated and further research will be done in next years to obtain data in intermediate and higher trophic states chl-a classes.
Our results show that the most suitable algorithms for the three combinations of platforms and sensors were those using NIR and red bands, and that 3BDA and its special case 2BDA were the most successful in retrieving chl-a with all the sensors tested. Beck et al. [71] obtained similar results when identifying the most "portable" algorithms for the quantification of chl-a. In their study, images from an airborne hyperspectral sensor (CASI) taken simultaneously with dense in situ water quality measurements, and synthetic multispectral datasets, yielded to the conclusion that the NDCI, 2BDA, and 3BDA chl-a algorithms were the most portable, being successful between hyperspectral imagers, WorldView-2 and -3, and Sentinel-2, MERIS-like imagers and, to a lesser degree, MODIS-like imagers.
The overall performance of our approach shows that the use of satellite data as baseline monitoring for eutrophication can be feasible bycombining different sensors, even in very small waterbodies. In this study, we found coherent results when applying the algorithms developed in low to medium chl-a values (1-20 µgr/L) but these failed in retrieving higher concentrations, which could be caused by the "packaging effect" (e.g., [38,72]) but also highlights the need for a further development of new algorithms calibrated with higher chl-a data. As it has been stressed in a number of studies, remote sensing reflectance data in waters with a high level of optical complexity are affected by the concentration of different optically active constituents and, as such, a single algorithm is unlikely to perform adequately in all situations [46,73]. This issue can be solved by using different algorithms in different conditions, which in our case should be algorithms in high chl-a concentration status, to be used alternatively depending on the limits of applicability stated for each model [7,74]. The choice of the appropriate algorithm for each case needs the availability of "a priori data" that could be obtained through the use of autonomous buoys placed in the most representative point of the lake, with real time data acquisition, ingested by the system to apply the suitable algorithm in each case.
The satellite part of our approach has been designed with freely available data, both of global coverage (satellite imagery freely distributed by ESA and NASA) and in situ point measurements (Regional scale. GAIA Portal [28]). If a compromise between the different agencies responsible for water sampling and analysis could be reached, so that the sampling points were equally representative and valid for use in the calibration and validation (cal/val) of RS products in space and time, and made available as FAIR data, the monitoring of these systems would be enhanced by the spatial component offered by RS through cal/val exercises on a global scale. In this regard, the autonomous buoy networks usually deployed by research institutions, where data are collected at high frequency and made freely available through grassroot networks of global coverage like the Global Lake Observatory Network (GLEON), are of particular interest. These buoys are usually located in central areas of water bodies, i.e., in free-water pixels, perfectly valid for its use as calibration and validation match-up points for RS products.

Conclusions
The first objective of this work was to evaluate the combined performance of the chl-a retrieval of three different sensors and remote sensing platforms as a multi-scale approach to be included in the standard monitoring of eutrophication processes in small-sized water bodies, in an area of high cloud cover, such as Galicia (NW Iberian Peninsula), by developing empirical models for the retrieval of chl-a concentration.
Our results show that the best performing models for all sensors and platforms were those using NIR and red bands and that the three platforms can be used in combination to maximize the temporal and spatial resolution offered by any of them independently. The best performing models were first the 3BDA algorithm for Sentinel MSI data and, as a second option in the absence of Sentinel 2 MSI data, the NDVI for Landsat 8 OLI, all of which were validated in a chl-a range of 1-40 µgr/L. Both the improvement in atmospheric corrections above water, especially the accounting for adjacency effects in the AC processors, should be prioritized for the future improvement in the development of this kind of tool for small-and medium-sized water bodies.
The multispectral commercial sensor Rededge Micasense gave the best results in low chl-a values (below 3 µgr/L) with a modified version of the 3BDA algorithm, while, in a high chl-a concentration (85-100 µgr/L), the 2BDA and SABI models were the more successful. However, in the range 1-100 µgr/L, the combinations of blue and green and red and green bands were the most effective in discriminating chl-a values with this sensor. Three of the developed algorithms were combined in a preliminary ensemble model. Further research will be done to validate these models with different ranges of chl-a data to apply them to intermediate and higher chl-a classes.
The developed workflow, in its initial design, can be easily included in the current monitoring scheme of the regional authorities, as the algorithms to be applied to satellite images were calibrated and validated with the data provided by them on a regular basis. With little effort, consisting in the change in the localization of some sampling points from the current location to another more centered in the waterbody and a sampling scheme coincident with satellite overpass (prioritizing Sentinel 2 overpass date and time), the accuracy of these empirical models is expected to improve, reducing the uncertainty related to spatiotemporal collocation and adjacency effects, thus improving their predicting ability. This highlights the value of freely available data sources for capacity building and developing of effective monitoring tools that could maximize the safety of our aquatic resources.