Monitoring water quality in two dammed reservoirs from multispectral satellite data

Providing relatively fine spatial resolution multispectral data, Landsat-8, Landsat-7 (L8 and L7, respectively) and Sentinel-2 (S2) from 2013 to 2018 have been used in this study for enabling high-frequency monitoring of water quality of two small (the smaller with an area of 1.6 km) freshwater dammed reservoirs. Located in Sardinia (Italy) and Crete (Greek), respectively, Mulargia and Aposelemis represent vital resources to supply drinking water in downstream valleys. A total of 400 cloud-free satellite images were turned into information on water quality by using an image processing chain implementing physically based methods for retrieving chlorophyll-a concentration (Chl-a), turbidity, Secchi disk depth (SDD) and surface water temperature. These estimates have been successfully validated (the lower Pearson correlation r was 0.88 for Chl-a) with 23 match-ups of in situ and satellite data. Results of the multitemporal analyses showed a decrease of SDD due to the increase of Chl-a in Aposelemis or an increase of turbidity in Mulargia. For both freshwater reservoirs, the satellite-derived trophic state index assigned both lakes to mesotrophic conditions. The results finally suggested the effectiveness of S2 and Landsat in increasing, for the latest investigated years, the frequency of observations. ARTICLE HISTORY Received 11 February 2019 Revised 26 June 2019 Accepted 28 October 2019


Introduction
Water resources are limited and are facing issues that are caused by over-exploitation, continuous human pressure and also climate changing, which could have serious consequences on water quality (Shevah, 2015), an essential resource for human health, ecosystems and the economy. Degradation of water quality can result in human exposure to harmful diseases and chemicals (Hu & Cheng, 2013), reduced productivity and diversity of ecosystems and damage to aquaculture, agriculture and other water-related industries (Kaiblinger et al., 2009).
With environmental pollution and decrease in water availability becoming an increasingly serious problem, the issue of water quality and quantity has attracted serious attention from the public and the government (Xin, Li, Finlayson, & Yin, 2015). Therefore, the number of artificial reservoirs is rapidly increasing owing to the growth of the world's economy and related energy and water needs (Zhang et al., 2019) such as domestic, agriculture, and energy, despite the costs to build a dam are getting higher every year (Ahmad, El-Shafie, Razali, & Mohamad, 2014). A primary concern of water authorities is to have information on current status of water quality, so different water quality evaluation methods have been developed and extensively applied (e.g. Gao, Birkett, & Lettenmaier, 2012;Gopal, Goel, Sharma, & Trivedy, 1981;Olsen, Chappell, & Loftis, 2012). Some of these studies (e.g. Robert et al., 2016) showed that a successful approach for water quality monitoring might be carried out through measures which include in situ, laboratory, and satellite measurements. While in situ sampling allows the analysis of a wide range of parameters at various depths to be performed, remote sensing technologies support frequent synoptic measurements and extend the ability to study remote waterbodies that cannot be visited regularly in a costeffective way (e.g. Bresciani, Stroppiana, Odermatt, Morabito, & Giardino, 2011;MacKay et al., 2009). In situ and remote detection techniques are complementary and overall support a wide range of possibilities for monitoring water quality, including integration into regional and local models for water quality prediction (CEOS, 2017).
The aim of this study is to investigate the use of Sentinel-2 (S2) and Landsat-7 (L7) and  for monitoring water quality in two freshwater reservoirs. Both satellites have sensors contain overlapping visible, near-infrared (NIR) and shortwave-infrared (SWIR) bands making both instruments capable of monitoring a range of water quality constituents in freshwater ecosystems De Keukelaere et al., 2018;Lehmann, Nguyen, Allan, & van der Woerd, 2018;Pinardi et al., 2018;Toming et al., 2016). In addition, a synergic use of VIS-NIR and TIR bands of L8 has been proved to be crucial in hydrodynamic modelling of mesoscale phenomena (e.g. Brando et al., 2015). Integration of L7, L8 and S2 satellite-based datasets then enables high-frequency revisit times of about 3 days at medium latitudes, where this study focuses on. The study area is in fact represented by two dammed reservoirs, of two islands of the Mediterranean Sea, which are providing drinkable waters to downstream valleys. By observing water quality parameters from 2013 to 2018, the study also aims to evaluate if significant variation of water quality in the study area occurred in recent years.

Study area
Two freshwater dammed reservoirs have been investigated in this study ( Figure 1 and Table 1). Both located in the Mediterranean region, they are Mulargia (Italy) and Aposelemis (Greece).
The Mulargia is a dammed reservoir build on Flumendosa river located in south Sardinia, the second-largest Italian island. With a surface area of 12 km 2 and a capacity of 347 hm 3 (Table 1), Mulargia serves as a drinking water source for the town of Cagliari with hinterland and for about 20 villages around in the greater area, summing up to a population of 700.000 inhabitants. The total annual abstraction for drinking water purposes is estimated to be 100 hm 3 . The Mulargia waters are in mesoeutrophic state as a consequence of loads of phosphorus mostly due to agriculture and zootechnic, with low transparency and high conductivity representing the most important issues.
The Aposelemis reservoir is also artificial; it is located in the north-eastern part of Crete island. With surface area of 166 km 2 and a capacity of 27 hm 3 (Table 1), Aposelemis reservoir serves as a drinking water source for the towns of Heraklion and Agios Nikolaos, as well as local communities in the greater area, summing up to a population of 300.000 inhabitants. The total annual abstraction for drinking water purposes is estimated to 17 hm 3 .

Satellite and field data
For assessing water quality in the study area, imagery acquired from multispectral optical sensors onboard of the polar orbiting satellites L7, L8 and S2 was used. The two S2 twin satellites A and B (S2A/B) carry the MSI (Multispectral Instrument) optical sensor, which acquires in 13 bands, in the spectral region between 442 and 2201 nm, with a spatial resolution of 10-20-60 m depending on bands. The L7 satellite carries onboard the ETM+ sensor (Enhanced Thematic Mapper Plus), with seven spectral bands in the spectral region between 441 and 2345 nm with a spatial resolution of 30 m, and a TIR band with a spatial resolution of 60 m. Finally, L8 equipped with OLI (Operational Land Imager), which acquires in 9 bands in the spectral region between 435 and 2294 nm with a spatial resolution of 30 m. L8 also carries onboard the Thermal Infrared Radiometer Sensor (TIRS) with a spatial resolution of 100 m. Table 2 shows the main characteristics of different satellite sensors used in this work.
For the period 2013-2018, 400 cloud-free satellite images have been processed to obtain water quality products, in particular maps of Chlorophyll-a concentration (Chl-a), turbidity, transparency and trophic status index (based on Carlson, 1977). In particular, the trophic status index has been also implemented by the Copernicus Global Land Service Lake Water Quality (Simis, Stelzer, & Müller, 2018). Furthermore, water surface temperature maps have been produced on TIR measurements. Table 3 shows the frequency of satellite acquisitions available for estimating LSWT and SDD, TUR and Chl-a over the multi-annual study period and for the four seasons. Since TIR measurements are only available on L7 and L8, it is clear that the improved revisiting for observing water quality paraments by using also S2 is not applicable to LSWT.
The physical methods implemented in the Modular Inversion and Processing System (MIP) (Heege & Fischer, 2004;Heege, Hausknecht, & Kobryn, 2007;Heege et al., 2014) were used for retrieval of water constituents from S2 and Landsat. MIP is a sensor independent image processing chain based on radiative transfer-based which couples atmospheric-adjacencysunlight correction images. It solves the Radiative Transfer Equation (RTE) using a Finite Element Model (FEM) in a multi-layer system, including atmosphere, water surface, waterbody and, in case of optically shallow waters, the bottom reflectance (Bulgarelli, Kisselev, & Roberti, 1999;Kisselev & Bulgarelli, 2004). All bidirectional dependencies of scattering and absorption properties in the atmosphere and the waterbody, as well as bidirectional reflectance and transmission of the water surface, are accounted for within one consistent RTE model. The model also considers the full range of possible geometrical conditions between sensor, sun and target   (i.e. the individual water pixels). The software version of MIP is fully based on physics and relates absorption and backscattering of in-water properties as a function of the sensor radiances. The accounted dependencies include variable target level altitude and observer altitude, variable atmospheric aerosol properties amongst other parameters, using the radiative transfer model FEM (Kiselev et al., 2015). The MIP architecture systematically handles the independent properties of sensor parameters and specific optical properties as well as the radiative transfer relationships (at 1nm spectral resolution). The processing includes the acquisition or harvesting of satellite data, radiometric calibration, spatial subsetting of region of interest, land-water-cloud masking (Heege et al., 2014). The MIP finally generates estimation on aerosol, in-water inherent optical properties of water constituent concentrations and creation of quality indicators and quality control. Figure 2 shows the workflow of image processing implemented in MIP.
The LSWT values were retrieved from top-ofatmosphere brightness temperature in the TIRS band 10 (10.9 μm) with the radiative-transfer-based atmospheric correction described in Barsi, Barker, and Schott (2003), Barsi et al. (2014). The model runs with input data on atmospheric profile on air temperature, pressure and relative humidity, preselected according to the latitudes of the target and the season. The method has been extensively adopted in multiple applications, including the retrieval of sea-surface temperature (Brando et al., 2015).
The image-derived products were compared to in situ measurements available for both reservoirs. In particular, a series of field campaigns were organised in coincidence to satellites overpasses in both reservoirs: on 10/11/2016, 05/06/2018 and 22/07/2018 for Mulargia, and 02-03-07/07/2018 for Aposelemis. Water was collected by sampling in the photic layers (measured with Secchi disk) for subsequent laboratory analysis. A total of 23 samples of concentration of Chla were determined spectrophotometrically according to Lorenzen (1967) and by using acetone 90% for the extraction of pigments; 19 samples of concentration of suspended particular matter (SPM) were determined gravimetrically according to Strömbeck and Pierson (2001). According to Braga et al. (2017), SPM data were also used as a proxy of turbidity as their relation might be generally assumed 1 to 1.
Water quality products were analysed to depict temporal trends of the major parameters over the two study sites. Mean and standard deviation for each parameter and date were computed for all valid pixels over the lakes' surface; to this aim, products derived from S2, L7 and L8 were resampled to common spatial resolution of 30 m. Only dates when more than 15% of the lake surface was covered were retained; indeed, if a smaller portion of the lake is mapped due to cloud cover, the output product might not be representative of average water conditions and statistics are affected by a greater rate of noise. Daily product time series were aggregated over the seasons (spring, summer, autumn and winter) based on the date. Given the high number of valid pixels over the season time span, in the boxplot outliers were removed; on average, outlier pixels represent a small proportion of the total number (less than 2%). For sake of clarity, time series data are plotted for LSWT, SDD and Turbidity, while for Chl-a, maps of mean values for the entire study period were computed to depict average conditions over the 2013-2018 time span. Maps of average conditions were accompanied by the coefficient of variation (CV = σ/μ).

Validation
The comparisons of satellite-derived water quality paraments and field data showed a good global agreement as presented in Table 4 with a Pearson correlation coefficient r in the range 0.87-0.94 with the lower value observed for Chl-a. This is likely due to the difference between the time of satellite overpass (instantaneous acquisition) and in situ sampling (typically performed in 3-4 hr). In fact, unlike SDD and turbidity, Chl-a is the parameter showing greater short-term variations due to intra-daily dynamics of phytoplankton. These diurnal periodic changes are particularly evident in summer when solar irradiation  and air temperature change quickly. Even if satellite remote sensing significantly improves the frequency of observation compared to in situ sampling, in highly dynamic water systems, overpass frequency of satellite observations from polar orbiting systems is still limited to capture diurnal or semi-diurnal cycles (Li et al., 2013).

Multi-temporal analysis for mulargia reservoir
Satellite-based products provided long-term information on LSWT, SDD and Turbidity for the entire dataset of Mulargia reservoir, as shown in the three panels of Figure 3. For each panel, the marker represents lake's surface average ± one standard deviation values and filling colour shows the source satellite mission (L7, L8 and S2). Since early 2016 (after the closing of the commission phase), the frequency of observation of water quality has been greatly improved by the availability of S2A data (yellow circles) further enhanced by the launch of S2B in spring 2017. For the year 2017 (full availability of S2 data) a total of 55 images were available for deriving maps of SDD and Turbidity. A separate issue is the estimation of LSWT, which relies on the availability of spectral bands in the thermal-infrared wavelengths, so far available only for Landsat (by considering that a moderate to high spatial resolution is needed in this study). The results clearly show the influence of turbidity levels on Secchi disk values with opposite trends (i.e. highest SDD peaks coincides with lowest turbidity values and vice versa); data shown in Figure 3 highlight maximum SDD and minimum turbidity values occurring in summer with transparent waters. On the contrary, spring and autumn seasons are generally characterized by greater turbidity as a consequence of meteorological extreme rainfall events (Caloiero & Veltri, 2018). Multitemporal data show that significant events of SDD maximum peaks occurred in 2018: in fact, out of a total of 12 dates with SDD > 3 m, six occurred in the first half of 2018. Anomalous values for summer 2018 are confirmed by the analysis of temporal trends by season, as shown in For LSWT, as expected, seasonal trends are evident through the years with highest (lowest) values in summer (winter) season (Figure 3 top panel). In particular, highest summer values were observed in July and August 2017 (average monthly LWST~26°C); for the same months and other years, average LWST was in the range 24-25°C.
Average estimated Chl-a over the study period was 10.8 mg/m 3 , with minima in winter seasons of about 0.8 mg/m 3 and maxima occurring in autumn 2017 (~30 mg/m 3 ). In the same period, a cyanobacteria bloom occurring over the lake surface and observed during in situ surveys (Planktothrix rubescens) determined high Chl-a values that reached local absolute maxima of 77.2 mg/m 3 close to the central area of the lake. Planktothrix is one of the most genus of cyanobacteria presents in the Mulargia (Buscarinu et al., 2003), the depth at which Planktothrix stratifies is therefore explained by buoyancy regulation in relation to the irradiance; in autumn 2017, the high cloudiness has driven migration of Planktothrix to the surface (Walsby, Ng, Dunn, & Davis, 2004). Highest average annual values for Chl-a were observed in 2016 and 2017 (16.7 and 17.6 mg/ m 3 , respectively) compared to previous years when average annual Chl-a was in the range 10.8-13.7 mg/m 3 .
The spatial analysis of Chl-a average distribution and variation (Coefficient of Variation, CV) over the lake's surface ( Figure 5) shows higher values in the northern regions, where Flumendosa river transports nutrient-rich waters into the reservoir. In this region, besides the river water inflow, low bathymetry facilitates resuspension of nutrients that are present in the sediments. Lowest Chl-a values can be observed in the southern areas, close to the Mulargia dam, where a water caption point is located.
The conversion of Chl-a concentration maps into trophic state index highlighted that the highest frequency of the pixels (about 5 × 10 6 ) was found to be of class 6 and the second class of frequency was class 7 with about 1 × 10 6 pixels.

Multi-temporal analysis for aposelemis reservoir
Figures 6 and 7 show the same trend analysis for the Aposelemis reservoir. By inspecting SDD and turbidity panels in Figure 6, it is noticeable how standard deviation (black bars) is much greater than Mulargia's  (cf. Figure 3). Standard deviation decrease with the onset of S2 time series, which is characterized by a greater spatial resolution with respect to Landsat sensors, thus reducing adjacency effects. In Aposelemis, two major significant events of turbidity peak occurred in January 2015 and 2017, with average values of 26 and 15 NTU, respectively. These values are in correspondence with minima in SDD observations, since high turbidity, as excepted, leads to loss of transparency.
Since the end of 2017, the reservoir witnessed a significant increase in turbidity levels with monthly average values above 8 NTU and a maximum of 13.8 NTU in June 2018. This positive trend is likely due to a decrease of water levels; reduced water levels also enhance the effect of wind on the resuspension of fine sediments from the bottom of the lake in the shallower areas of the lake. Lowering of water level is a factor that affects turbidity in many lakes (Giardino, Bresciani, Villa, & Martinelli, 2010;Lisi & Hein, 2019)    . Summary statistics visualized as boxplots by year and season of the three water quality parameters in Aposelemis reservoir. In each boxplot, the median value is the central bold line, the extremes of the rectangle are given by the first (Q1) and third (Q3) quartiles (the 25th and 75th percentiles), whiskers are computed with the following rules: Q1-1.5*IQR and Q3+1.5*IQR, where the inter-quartile range IQR = Q3 -Q1. Outliers were removed representing less than 2%.
period, the temporal trend of SDD and Turbidity shows a slight tendency to loss of transparency.
Looking at seasonal trends shown in Figure 7, the increase through years for turbidity is clearly visible for spring and autumn whereas a more variable behaviour can be observed for the other two seasons. This trend is confirmed for SDD (decreasing trend over spring values) while, as discussed above, more variable and less consistent behaviour is observed for the other seasons.
For what concerns LWST, maximum summer values were observed in 2017; notice that higher values of 2018 in Figure 7 represent estimations only for the month of July hence they are not representative of the entire summer period. In 2014 and 2016, winter statistics show greater variability with respect to the other years of the investigated time series.
Finally, Figure 8 shows average Chl-a and coefficient of variation, CV, as computed over the study period. For Lake Aposelemis, 160 cloud free images show that average Chl-a was 16 mg/m 3 (min 1.7 and max 45.2) with reduced spatial variability within the reservoir (Figure 8). In the central regions, slightly higher average Chl-a values are combined with a high coefficient of variation (CV > 75%). The lowest values were found at the northernmost border close to Aposelemis dam, where the water caption point is located.
The conversion of Chl-a concentration maps into trophic state index highlighted that the highest frequency of the pixels (about 4.4 × 10 4 ) was found to be of class 6 and the second class of frequency was class 7 with about 3.4 × 10 4 pixels.

Conclusions
Landsat-7 & 8 and Sentinel-2A & B have been used to observe a set of key parameters for monitoring water quality. Satellite-derived products, obtained with physically based image processing chains have been providing multi-temporal and multi-scale data over two small freshwater reservoirs from 2013 to the present. The investigated temporal range also clearly showed how, since early 2016, the frequency of observation of water quality (apart for LSWT) has been greatly improved by the availability of S2A data further enhanced by the launch of S2B in spring 2017.
Despite challenges (due to e.g. adjacency effects) for retrieving water quality parameters over small to medium size lakes (i.e. 1.6 km 2 for Aposelemis, which is approximately 10 times smaller than Mulargia), the adopted physics-based models produced accurate results. In particular, satellite and field data show a good correlation. Results of the multi-temporal analyses showed a worsening in water quality due to the combination of increasing Chl-a and turbidity for a decreasing of Secchi disk depth. This was mostly occurring in concomitance with low water levels and intense precipitation events, the latter that foster a high run-off of particulate matter from the basin. For both freshwater reservoirs, the trophic state index evaluated from satellite data assigned both lakes to mesotrophic conditions. Satellite products described in this study were distributed to water managers of Mulargia and Aposelemis basins via web map server for supporting lake water management plans of both freshwater reservoirs.

Disclosure statement
No potential conflict of interest was reported by the authors.