Di ﬀ use Attenuation of Clear Water Tropical Reservoir: A Remote Sensing Semi-Analytical Approach

: The di ﬀ use attenuation coe ﬃ cient of downwelling irradiance ( K d ) is an essential parameter for inland waters research by remotely sensing the water transparency. Lately, K d semi-analytical algorithms substituted the empirical algorithms widely employed. The purpose of this research was to reparametrize a semi-analytical algorithm to estimate K d and then apply it to a Sentinel-2 MSI time-series (2017–2019) for the Tr ê s Marias reservoir, Brazil. The results for the K d semi-analytical reparametrization achieved good accuracies, reaching mean absolute percentage errors (MAPE) for bands B2, B3 and B4 (492, 560 and 665 nm), lower than 21% when derived from in-situ remote sensing reﬂectance ( R rs ), while for MSI Data, a derived MAPE of 12% and 38% for B2 and B3, respectively. After the application of the algorithm to Sentinel-2 images time-series, seasonal patterns were observed in the results, showing high K d values at 492 nm during the rainy periods, mainly in the tributary mouths, possibly due to an increase in the surface runo ﬀ and inﬂows and outﬂow rates in the reservoir watershed.


Introduction
Três Marias (TM) reservoir is an essential tropical inland water body situated in the southeast of Brazil (Minas Gerais State), in the São Francisco River upper valley [1,2]. The main use of TM reservoir is the hydroelectric power generation [3], with an installed capacity of 45 MW. Moreover, TM reservoir provides other ecosystem services such as water supply, irrigation, fishing, and fish-farming [4]. Shifts in rainfall precipitation and surface runoff seasonal variability over the watershed induce significant changes in the water quantity and quality, affecting the management and multiple uses of water and land over the reservoir basin [5][6][7]. With clear waters, TM reservoir might vary from oligo to mesotrophic state waters [8][9][10] depending on the season of the year. Therefore, it is essential to understand the spatiotemporal dynamics of water quality in TM to provide proper management of those waters. Changes in the water optical active constituents (OACs) such as phytoplankton, suspended sediments, and Colored Dissolved Organic Matter (CDOM) modify water transparency, affecting the light attenuation on the water column [11]. As this attenuation is intrinsically connected to  The connection between AOPs and IOPs is the key of understanding how to link K d variations with remote sensing (RS) measurements. This connection is frequently made by empirical or semi-analytical algorithms. For an extended period, the standard method to estimate K d from remotely sensed data has been by empirical approaches using the relationship between K d and remote sensing reflectance Remote Sens. 2020, 12, 2828 3 of 23 (R rs ) [21,22]. However, uncertainties related to the dataset used for reparametrization and validation of empirical approaches are inherent to this method, making it insufficient to provide an accurate estimation of the K d for ranges not used in the reparametrization process [23]. To overcome this issue, semi-analytical approaches based on radiative transfer equations have been proposed [13,[24][25][26], improving the accuracy of K d estimations. The semi-analytical approaches require IOPs as input (absorption and backscattering coefficients), which is done by AOPs inversion. By combining the semi-analytical K d algorithms with existing semi-analytical algorithms for the derivation of the water total absorptions and total backscattering coefficients (water plus constituents), K d can easily be implemented using satellite data [25,27]. Despite the successful use of IOPs inversion algorithms in oceanic and coastal waters, their use in inland waters still requires a reparameterization of the approaches used to provide a high-performance algorithm.
Therefore, the main purpose of this research was to evaluate the performance of Lee et al. [26] semi-analytical algorithm to derive K d from R rs using in-situ and Sentinel-2 MSI data. As the algorithm requires IOPs as input, the Quasi-Analytical Algorithm (QAA) was reparametrized to TM reservoir using IOPs in-situ measurements. The validated K d algorithm was, therefore, applied to a time-series of Sentinel-2 MSI images to understand the spatial and temporal variability of light attenuation in TM reservoir and its relationship with precipitation events. Before that, a brief characterization of TM reservoir waters was done to support the K d results discussions.

Materials and Methods
The development of this study is organized in four steps ( this issue, semi-analytical approaches based on radiative transfer equations have been proposed [13,[24][25][26], improving the accuracy of estimations. The semi-analytical approaches require IOPs as input (absorption and backscattering coefficients), which is done by AOPs inversion. By combining the semi-analytical algorithms with existing semi-analytical algorithms for the derivation of the water total absorptions and total backscattering coefficients (water plus constituents), can easily be implemented using satellite data [25,27]. Despite the successful use of IOPs inversion algorithms in oceanic and coastal waters, their use in inland waters still requires a reparameterization of the approaches used to provide a high-performance algorithm.
Therefore, the main purpose of this research was to evaluate the performance of Lee et al. [26] semi-analytical algorithm to derive from using in-situ and Sentinel-2 MSI data. As the algorithm requires IOPs as input, the Quasi-Analytical Algorithm (QAA) was reparametrized to TM reservoir using IOPs in-situ measurements. The validated algorithm was, therefore, applied to a time-series of Sentinel-2 MSI images to understand the spatial and temporal variability of light attenuation in TM reservoir and its relationship with precipitation events. Before that, a brief characterization of TM reservoir waters was done to support the results discussions.

Materials and Methods
The development of this study is organized in four steps (

Study Area
The study area consists of the Três Marias (TM) hydroelectric reservoir located in the upper valley of the São Francisco River, in the central region of Minas Gerais state ( Figure 2) [1,2]. Constructed in 1957 and operating since 1962, TM reservoir is approximately 100 km long in its main axis along the São Francisco River and has a volume of at least 20 billion m 3 of water [3]. Moreover, the TM reservoir has clear water characteristics, except at the entrance of its tributaries and the banks due to its shallow depth favoring sediment resuspension and bottom effect transparency measurements [28,29]. TM reservoir watershed is characterized by a humid tropical climate with average annual rainfall precipitation of 1400 mm of [30] concentrated between October and March (85%) [31]. Most of its watershed is inside the Cerrado biome (Brazilian Savannah), with some fragments of the Mata Atlântica biome (Atlantic Forest) [32]. The main rivers feeding TM reservoir are São Francisco, Pará, and Paraopeba Rivers ( Figure 2a) [33]. Furthermore, TM watershed main land uses and land cover consist of cattle grazing, patches of the Atlantic forest, urban areas, and a relevant mining industry activity [32].
average annual rainfall precipitation of 1400 mm of [30] concentrated between October and March (85%) [31]. Most of its watershed is inside the Cerrado biome (Brazilian Savannah), with some fragments of the Mata Atlântica biome (Atlantic Forest) [32]. The main rivers feeding TM reservoir are São Francisco, Pará, and Paraopeba Rivers (Figure 2a) [33]. Furthermore, TM watershed main land uses and land cover consist of cattle grazing, patches of the Atlantic forest, urban areas, and a relevant mining industry activity [32].
Most of the mining sites inside TM watershed are located within Paraopeba and Pará Rivers subbasins. Most of them are located along the Paraopeba River, with particular attention to the region between Brumadinho city and Belo Horizonte Metropolitan Region, and Pará River sub-basin ( Figure  2a). This area suffered a huge disaster related to the breaking of a mining tailings dam in Brumadinho city on the 25 th of January 2019, which caused at least 259 deaths [34]. This event highlighted the importance of the Paraopeba river since it drains all the runoff waters from Brumadinho city and part of the Belo Horizonte Metropolitan Region [32,35].
The water transparency and the trophic state of the TM reservoir show seasonal changes with lower concentrations of ℎ and TSS during the dry season [8,10,36]. Generally, the TM reservoir waters are more transparent in the middle of its main channel, having the lowest transparencies at the entrances of its main tributaries. In general, São Francisco River waters that flow into the TM have higher TSS concentrations and lower transparency than those of Paraopeba River [37]. This difference may be explained by the existence of a small-scale reservoir (Retiro Baixo) just upstream of the Paraopeba River entrance that holds part of the sediments which, otherwise would be transported into TM reservoir.  Most of the mining sites inside TM watershed are located within Paraopeba and Pará Rivers sub-basins. Most of them are located along the Paraopeba River, with particular attention to the region between Brumadinho city and Belo Horizonte Metropolitan Region, and Pará River sub-basin ( Figure 2a). This area suffered a huge disaster related to the breaking of a mining tailings dam in Brumadinho city on the 25 January 2019, which caused at least 259 deaths [34]. This event highlighted the importance of the Paraopeba river since it drains all the runoff waters from Brumadinho city and part of the Belo Horizonte Metropolitan Region [32,35].
The water transparency and the trophic state of the TM reservoir show seasonal changes with lower concentrations of chla and TSS during the dry season [8,10,36]. Generally, the TM reservoir waters are more transparent in the middle of its main channel, having the lowest transparencies at the entrances of its main tributaries. In general, São Francisco River waters that flow into the TM have higher TSS concentrations and lower transparency than those of Paraopeba River [37]. This difference may be

Field Data Acquisition
Data collection was carried out at 42 sampling stations (see Figure 2b Radiometric data were acquired using inter-calibrated TriOS-RAMSES radiometers (Rastede, Germany) to measure (θ; ϕ; λ), (θ′; ϕ; λ), (λ) and (z; λ), within 350-900 nm wavelength range. The framework for the four radiometers ( Figure 3) was: (i) (λ) within a nadir angle (θ) of 45° and Sun azimuth angle (ϕ) of 135°; (ii) (λ) within a zenithal angle (θ′) of 45° and azimuthal angle of 135°; (iii) (λ) at 180° in relation to the water surface; and, (iv) (z; λ) positioned in a cage and submerged with the aid of a crane. The above-water sensors were positioned 2 m away from the water surface to avoid shadow and vessel reflections. This acquisition geometry is based on Mobley [38] to avoid glint effects. The total absorption coefficient (λ) was measured with the Wetlabs Spectral Absorption and Attenuation meter (ACS, Bellevue, United States) [39] during the 2013 campaign, while in 2019, it was determined in the laboratory as described in Section 2.3.
Water sampling collection protocols followed the American Public Health Association (APHA) [40]. Water samples were collected at approximately 30 cm of depth using sampling bottles coated with opaque material and then stored inside a cooler with ice until its filtration. After the samples collection, the filtration was conducted to prepare the water samples and the filters for laboratory analysis and determination of ℎ , TSS, (λ), (λ), and its fractions (λ) and (λ). Filter Whatman GF/F was used for chla, (λ), and its fractions (λ) and (λ), and Whatman GF/C for TSS, while water samples filtered using Nylon filters with pore diameter of 0.22 µm were used for CDOM. The chla and TSS samples were frozen inside envelopes at temperatures below 0°C, while (λ) and fractions samples were frozen inside opaque vials in a liquid nitrogen bottle at temperatures below −196 °C. The CDOM water samples were chilled without freezing it.

Optically Active Constituents
The ℎ concentration was determined using the Nush method [41] based on the 665 and 750 nm absorbances read in a spectrophotometer (brand Agilent/Varian, model Cary-50 Conc, Santa Clara, United States). This method extracts the phytoplankton pigments retained in the Whatman GF/F filter by warm ethanol followed by a thermal shock. The extracted pigments were put to rest The total absorption coefficient a(λ) was measured with the Wetlabs Spectral Absorption and Attenuation meter (ACS, Bellevue, United States) [39] during the 2013 campaign, while in 2019, it was determined in the laboratory as described in Section 2.3.
Water sampling collection protocols followed the American Public Health Association (APHA) [40]. Water samples were collected at approximately 30 cm of depth using sampling bottles coated with opaque material and then stored inside a cooler with ice until its filtration. After the samples collection, the filtration was conducted to prepare the water samples and the filters for laboratory analysis and determination of chla, TSS, a CDOM (λ), a p (λ), and its fractions a φ (λ) and a d (λ). Filter Whatman GF/F was used for chla, a p (λ), and its fractions a φ (λ) and a d (λ), and Whatman GF/C for TSS, while water samples filtered using Nylon filters with pore diameter of 0.22 µm were used for CDOM. The chla and TSS samples were frozen inside envelopes at temperatures below 0 • C, while a p (λ) and fractions samples were frozen inside opaque vials in a liquid nitrogen bottle at temperatures below −196 • C. The CDOM water samples were chilled without freezing it.

Optically Active Constituents
The chla concentration was determined using the Nush method [41] based on the 665 and 750 nm absorbances read in a spectrophotometer (brand Agilent/Varian, model Cary-50 Conc, Santa Clara, United States). This method extracts the phytoplankton pigments retained in the Whatman GF/F filter by warm ethanol followed by a thermal shock. The extracted pigments were put to rest for 6 to 24 h inside a fridge protected from light before the measurements. Finally, sample extractions were read in the spectrophotometer using quartz cuvettes with 50 mm of the optical path.
TSS concentrations were computed using pre-calcined and pre-weighed Whatman GF/C filters. Before laboratory analysis, the filters were thawed in the laboratory and dried using a laboratory drying-oven. After that, dried filter samples were weight-measured using a precision scale [42] in a room with controlled air conditions. Finally, the concentrations were calculated using the relation between the filter retained solid weight and the volume of filtered water samples [42].

Apparent Optical Properties
The AOP parameters derived from the radiometric field measurements used in this study were the R rs (λ) and K d (λ). The Remote Sensing Reflectance was calculated using Equation (1).
where ρ(θ; φ) is the reflectance factor at the air-water interface. The values for each field station were obtained from Mobley [38] Look-Up- Table, based on in-situ wind speed and geometry acquisition. For each sample station, around 150 measurements of L t (θ; φ; λ), L sky (θ ; φ; λ) and E s (λ) were obtained. At each field station, the representative spectrum was determined by the closest to the median spectrum, as described in Maciel et al. [43]. The Diffuse Attenuation Coefficient (K d ) was computed from the profiles of downwelling planar From hereafter, K d−measured will represent the field measured K d (λ).

Inherent Optical Properties
CDOM water samples were read using a UV/VIS spectrophotometer (brand Shimadzu, model UV2600, Kyoto, Japan). Spectrophotometer readings were preprocessed to subtract the mean value of the measurements of spectral absorbance (A(λ)) and offset [45,46], and then transformed into a CDOM (λ) [47] using a relationship between A(λ) and the cuvette optical path length (l) (Equation (5)).
The measurements of a(λ) using ACS equipment, during the 2013 campaign, were acquired for the water column profile at each sampling station [29]. The data were smoothed to eliminate the measurement noise. Since the water column was not stratified in most of the sampling stations, the median of the a(λ) profile was used to represent a(λ) for the whole water column.
Finally, to assess a p (λ), and its fractions (a d (λ) and a φ (λ)) for the 2019 campaign, Whatman GF/F filters with retained material were analyzed using the UV/VIS spectrophotometer (Shimadzu UV2600) with dual-beam and an integrating sphere following the methodology proposed by Tassan & Ferrari [48,49]. Measurements of optical density (OD) were determined and corrected for the particulate material to estimate a p (λ) in the first reading [49,50]. After that, the pigments in the filter were leached with sodium chloride and then read again to determine the detritus corrected OD and estimate a d (λ) [49,50]. This approach measures the spectral reflectance and transmittance to estimate the a p (λ) and a d (λ), and then a φ (λ) by subtracting a d (λ) from the a p (λ) [42,48,49].

Sentinel-2 MSI Data
MultSpectral Imager (MSI) images from Sentinel-2 A and B were used in the useful spectral range for assessing water transparency (Table 2) including the visible (B2 to B5) and near-infrared (B6 and B8) bands, and the shortwave infrared region (B11) needed for glint effect correction [51,52]. These images were used to assess temporal changes in K d , resulting from both natural processes and impacts of the Brumadinho disaster on TM reservoir.

Kd(λ) Semi-Analytical Algorithm
In this study, we used the relationship between the diffuse attenuation coefficient (K d ) and IOPs [54,55] to reparametrize and apply a semi-analytical algorithm of K d (λ) (K d−SA (λ)-Equation (6)) to our study site, using a(λ) and b b (λ) data as input [26]. From hereafter, the semi-analytical K d (λ) will be named as K d−SA . where, m 0 , m 1 , m 2, and m 3 are constants independent from wavelength and water type, θ s is the zenithal sun angle, γ a parameter that represents the effect caused by the variation in the scattering constituents in the water, and η w (λ) the parameter that quantifies the contribution of b bw (λ) over the total b b (λ) (Equation (7)).
K d−SA (λ) was used with Lee et al. [26] original values for m 0 , m 1 , m 2 , m 3 and γ (0.005, 4.259, 0.52, 10.8 m and 0.265 respectively), while θ s was calculated using the geographic coordinates, day of the year, and hour of data acquisition.

Deriving Inherent Optical Properties from Remote Sensing Reflectance
To derive a(λ) and b b (λ) as input for equation 6, the empirical steps of the original QAA in its version 6 (QAAv6) [24,56] were reparametrized using λ 0 , the reference wavelength, at 560 nm (although the use of λ 0 = 670 nm in the original QAAv6 (Table 3-steps 2 and 4). 1 The same as the QAAv6, where the r rs (λ) is derived analytically from R rs (λ); u(λ) can be explained from the relation between a(λ) and b b (λ), and derived from r rs (λ), g 0 and g 1 values.
Reparametrizated for MSI bands, choosing the best three bands ratio for this case.

5
The same as the QAAv6.
The same as the QAAv6.
We named this reparametrized QAA as QAA TM. For this reparametrization, only a(λ) derived from ACS [39] measurements carried out during the 2013 campaign (n = 22) were used as no ACS measurements were available for the 2019 campaign.
For empirical step 2, an exponential adjustment was carried out between a(λ 0 = 560), derived from ACS, and a three-band ratio to derive the parameters M and N (Equation (8)). The wavelengths used in adjusting correspond to the central wavelength of MSI B3, B4, and B5 bands [57]. The R rs used in the bands ratio were derived from in situ measurements at each of the 22 sampling stations.
where u(λ) is the ratio of b b to (a + b b ) derived from subsurface remote sensing reflectance (r rs (λ)) using QAA step 1 equation [24,56] (Table 3), and b bw is the backscattering coefficient of pure water [58].
To obtain η from r rs , similarly to step 4 of QAAv6, we adjusted Equation (11), based on the stations η derived from Equation (10) and a two-band ratio using the central wavelength of MSI B4 and B5 bands.
Finally, the parametrized QAA TM was applied using R rs (λ) data to compute a(λ) and b b (λ) used as input for K d−SA .

Semi-Analytical Algorithm Validation
The reparametrized semi-analytical algorithm K d−SA (λ) was applied and validated using two datasets as input: (i) first with the in situ R rs (λ) data; and then, (ii) with the R rs (λ) derived from MSI image bands.
Both K d−SA (λ) derived from in situ R rs (λ) and K d−SA (λ) derived from MSI Sentinel-2 A image, acquired in 1 July 2019 coincident with in situ samplings, were validated with 20 sampling stations of 2019 campaign. The validations were performed from blue to red MSI bands (B1 to B4) for each dataset.
Since the MSI revisit time is 5 days, we used an image covering TM in the middle of the 2019 campaign (1 July). For this reason, the in situ measurements may differ up to 1 day from the satellite image. This difference should not be of great concern because there was no rainfall precipitation during the survey campaign, which could change the OACs.
For each algorithm and each dataset, the estimated K d−SA (λ) was compared with K d−measured (λ) using a linear regression approach to assess algorithm accuracy (Section 2.6). The reparametrized QAA TM algorithm was applied to the MSI image time series, as shown in Section 2.7.

Data Analysis and Accuracy Assessment
For the validation of K d−SA algorithm, the K d−measured was used as a reference and the following statistical metrics were calculate for performance evaluation: correlation coefficient (R 2 -Equation (12)); the mean absolute percentage error (MAPE-Equation (13)); and, the Root-Mean Square Error (RMSE-Equation (14)).
where, y i is the reference field value,ŷ i the predicted value, y i the average of the field reference value, and n the size of the dataset. All the linear regressions and the accuracy assessment for the validations were done based on the MSI bands B1 to B5 wavelengths.

Time Series Kd Algorithm Application
The most accurate algorithm was applied to a 2017-2019 time series of Sentinel-2 A and B imagery for mapping time changes in TM spectral K d−SA (λ). Moreover, images with less than 10% of cloud cover over the scene were used in the query searches to reduce the data volume.
First, the image collection was downloaded and preprocessed (Section 2.7.1). After the preprocessing, the K d−SA (λ) the algorithm was applied to the image dataset and then monthly-average filtered.
Finally, the results for the time series were extracted for three sampling points distributed along reservoir mean inflow and outflow rates for each month of the period. Eventually, months without images following the pre-requirements were not used, causing gaps in the time series.

Sentinel-2 MSI Time Series Preprocessing
First, the image dataset was loaded from the internet and processed using an internal algorithm from the Instrumentation Laboratory of Aquatic Systems (LabISA) that downloads the images at the top of the main axis of TM reservoir (downstream, middle, and upstream) using Quantum GIS software. The extracted results were confronted with the accumulated rainfall precipitation for the watershed and the atmosphere level (L1C) from Google Earth Engine API, then applies AtmosPy atmospheric correction algorithm [59] achieving surface level (L2C). AtmosPy is a processing interface to run the Second Simulation of the Satellite Signal in the Solar Spectrum (6S-6SV) atmospheric correction algorithm [60,61]. In the Atmospy interface, the 6S algorithm receives atmospheric products from Multi-Angle Implementation of Atmospheric Correction (MAIAC) like water vapor, ozone, and aerosol optical depth. After atmospheric correction, the LabISA routine applies the fmask algorithm [62] to produce a water mask without clouds and cloud shades.
After atmospheric correction and masking, the image dataset was standardized to 10 m by resampling bands B5, B6, and B11 using the nearest neighbor method to match VIS bands' spatial resolution. Finally, the glint effect was removed, subtracting the band B11 values from all the remaining five bands [51].

In Situ Limnology Data
The subsurface chla and TSS concentrations, and the Z sd measurements have different ranges for each campaign (Table 4)  The descriptive statistics of chla, TSS, and Z sd variation for each campaign are in Table 4.

Apparent Optical Properties
Field R rs (λ) values showed higher amplitudes for the 2013 campaign (Figure 4a) than those for 2019 (Figure 4b), highlighting P25 and P26 stations during 2013 (Figure 2) with the highest R rs (λ) values near 560 nm (Figure 4a-dotted lines). Both stations presented high TSS values (2.36 and 6.41 g m −3 respectively), and low Z sd values (0.85 and 0.50 m respectively). In general, R rs (λ) peaks have occurred in the green region of the spectrum (between 550 to 580 nm) reaching values from 0.005 up to 0.042 sr −1 , and two other little "shoulders" in the red region, one at 620 nm and the other between 676 and 679 nm, approximately. Moreover, the R rs (λ) in the NIR region (700 to 750 nm) were the lowest ones, reaching null (0) values.  For the calculation, it is necessary to understand the (z; λ) decay as the water column depth increases. Figure    For the K d−measured calculation, it is necessary to understand the E d (z; λ) decay as the water column depth increases. Figure 5 presented one example of E d profile in the water column, in a logarithmic scale, for each campaign for blue (443 nm), green (560 nm), and red (665 nm) wavelengths. For both campaigns, the highest light attenuation occurred at 443 nm, when it reached 1% of the subsurface E d (z; λ) at approximately 6 m depth. The red region followed the blue region with penetration of approximately 7 m of the subsurface E d (z; λ) at 665 nm. Finally, the highest light penetration occurred at 560 nm (green) with approximately 13 m of euphotic depth.  For the calculation, it is necessary to understand the (z; λ) decay as the water column depth increases. Figure   The spectral decay of (z; λ) with depth in the water column ( Figure 5) was used to calculate the (λ) (Figure 6). According to in-situ limnology and data, higher The spectral decay of E d (z; λ) with depth in the water column ( Figure 5) was used to calculate the K d−measured (λ) (Figure 6). According to in-situ limnology and R rs data, higher K d−measured (λ) values were observed for the 2013 campaign. The K d−measured (λ) in the blue region (< 500 nm) reached up to approximately 1.5 m −1 at 443 nm. Moreover, 2013 data also presented a high range of values compared to those of the 2019 campaign, which is clearly due to the high variability in field-measured OACs. It is also important to highlight that stations P25 and P26 (Figure 6-Dotted lines) presented uncertainties in K d−measured (λ) for wavelengths lower than 500 nm, which can be attributed to high CDOM and TSS absorption at those wavelengths [63]. Similar to the (z; λ) decay ( Figure 5), Figure 6 shows that the water column light penetration was higher in the green region of the spectrum (550-600 nm) for all stations, which is represented by (λ) lower values.

Inherent Optical Properties
The measurements of in-situ total absorption coefficient ( (λ)) values ( Figure 7) followed the observed results for and limnological data. Note that for 2013 the (λ) were acquired through ACS measurements whereas for 2019 (λ) data were obtained through a spectrophotometer (See Section 2.2.3). Besides, when both were compared, the results follow the same pattern. All (λ) spectra showed ℎ absorption features at approximately 676 to 679 nm wavelengths. From Figure  7 it is possible to note that the lowest (λ) occurred at 550-580 range, following that observed for . Moreover, one can note that as the wavelength increases the contribution of pure water absorption ( (λ)) increases in the total (λ), (Figure 7-Red-dotted line), with approximately 80% in the red domain (665 nm) and with almost 100% at wavelengths higher than 700.  Similar to the E d (z; λ) decay ( Figure 5), Figure 6 shows that the water column light penetration was higher in the green region of the spectrum (550-600 nm) for all stations, which is represented by K d−measured (λ) lower values.

Inherent Optical Properties
The measurements of in-situ total absorption coefficient (a(λ)) values ( Figure 7) followed the observed results for R rs and limnological data. Note that for 2013 the a(λ) were acquired through ACS measurements whereas for 2019 a(λ) data were obtained through a spectrophotometer (See Section 2.2.). Besides, when both were compared, the results follow the same pattern. All a(λ) spectra showed chla absorption features at approximately 676 to 679 nm wavelengths. From Figure 7 it is possible to note that the lowest a(λ) occurred at 550-580 range, following that observed for K d . Moreover, one can note that as the wavelength increases the contribution of pure water absorption (a w (λ)) increases in the total a(λ), (Figure 7-Red-dotted line), with approximately 80% in the red domain (665 nm) and with almost 100% at wavelengths higher than 700.  Similar to the (z; λ) decay ( Figure 5), Figure 6 shows that the water column light penetration was higher in the green region of the spectrum (550-600 nm) for all stations, which is represented by (λ) lower values.

Inherent Optical Properties
The measurements of in-situ total absorption coefficient ( (λ)) values ( Figure 7) followed the observed results for and limnological data. Note that for 2013 the (λ) were acquired through ACS measurements whereas for 2019 (λ) data were obtained through a spectrophotometer (See Section 2.2.3). Besides, when both were compared, the results follow the same pattern. All (λ) spectra showed ℎ absorption features at approximately 676 to 679 nm wavelengths. From Figure  7 it is possible to note that the lowest (λ) occurred at 550-580 range, following that observed for . Moreover, one can note that as the wavelength increases the contribution of pure water absorption ( (λ)) increases in the total (λ), (Figure 7-Red-dotted line), with approximately 80% in the red domain (665 nm) and with almost 100% at wavelengths higher than 700.  For CDOM measurements, a CDOM (λ) at 443 nm (a CDOM (443)) ( Figure 8) and spectral a CDOM (λ) (Figure 8) achieved higher variation for the 2013 campaign (~0.1 to up to 0.8 m −1 ) compared to the 2019 (~0.1-0.3 m −1 ), mostly due to the stations inside TM main affluent rivers (P13, P23, P25, and P26), represented as blue dotted-lines in Figure 8a. Besides, these stations presented different slopes of CDOM decay, which can be attributed to different sources of organic matter from those tributaries. The a CDOM (443) statistical parameters (mean, median and standard deviation) for each campaign are in Figure 8.

QAATM Parametrization
The exponential adjustment between (560) and the three bands  Nevertheless, the results from the estimation of the η by equation 10 used as preprocessing for step 4 reparametrization were satisfactory, with an average R 2 of 0.77. Table 5 summarizes the empirical steps of QAATM and QAA. The parametrized QAATM as the QAAv6 adapted for MSI bands used to derive (λ) and (λ) from ( ) values are summarized in Table 5. Higher a CDOM (443)values ( Figure 8) were found at the sampling stations with higher R rs (λ) (Figure 4) for wavelengths in the blue to green regions of the spectrum, between approximately 400 and 560 nm.

QAA TM Parametrization
The exponential adjustment between a(560) and the three bands R rs (560)/R rs (665) + R rs (704) (Figure 9a) showed good results, with R 2 of 0.93. Besides that, the linear adjustment between the η and the exponential of the two bands ratio R rs (665)/R rs (704) (Figure 9b

QAATM Parametrization
The exponential adjustment between (560) and the three bands  Nevertheless, the results from the estimation of the η by equation 10 used as preprocessing for step 4 reparametrization were satisfactory, with an average R 2 of 0.77. Table 5 summarizes the empirical steps of QAATM and QAA. The parametrized QAATM as the QAAv6 adapted for MSI bands used to derive (λ) and (λ) from ( ) values are summarized in Table 5. Nevertheless, the results from the estimation of the η by equation 10 used as preprocessing for step 4 reparametrization were satisfactory, with an average R 2 of 0.77. Table 5 summarizes the empirical Remote Sens. 2020, 12, 2828 14 of 23 steps of QAA TM and QAA. The parametrized QAA TM as the QAAv6 adapted for MSI bands used to derive a(λ) and b b (λ) from R rs (λ) values are summarized in Table 5.

In Situ K d−SA (λ) Validation
The QAA TM (Figure 10a) showed better overall results for the spectral regions with the lowest

In Situ ( ) Validation
The QAATM (Figure 10a  The results of (λ) were also statistically assessed from B1 to B4 (Figure 11) to compare QAATM accuracy with the non-parametrized QAAv6 accuracy. For QAATM, the results for bands B2 to B4 showed good accuracies (Figure 11a Regarding the slope and offset of the linear regression, QAA TM showed a steeper slope and slightly similar offset when compared to QAAv6 (Figure 10). The different reference wavelength used in each QAA approach has a visible interference in the slope of the regressions and the accuracy of the results for each chosen band. QAA TM approach had better accuracies for the lower values of K d−SA (λ), agreeing most with the bands at 492, 560 and 665 nm (B2 to B4), while it overestimates higher values of K d (λ) at 443 and 704 nm (B1 and B5). Besides that, QAAv6 had a better agreement with values at 560 and 665 nm (B3 and B4), overestimating the results of K d−SA (λ) at 443 and 492 nm and overestimating results at 704 nm (B5).
The results of K d−SA (λ) were also statistically assessed from B1 to B4 (Figure 11) to compare QAA TM accuracy with the non-parametrized QAAv6 accuracy. For QAA TM , the results for bands B2 to B4 showed good accuracies (Figure 11a), with R 2 of 0.77, 0.80 and 0.27; MAPE of 15% and 9%; and RMSE of 0.11, 0.04 and 0.07 m −1 , for each band, respectively; while B1 showed lower accuracy with R 2 of 0.18, MAPE of 50% and RMSE of 0.45 m −1 . The accuracy results for the non-parametrized QAAv6 (Figure 11b) were worse when compared to those for QAA TM , with R 2 values of 0.27, 0.79, 0.82, and 0.35; MAPE of 16, 23, 19, and 12%; and RMSE of 0.16, 0.14, 0.08, and 0.11 m −1 , for bands B1 to B4 respectively. These statistical results, slopes and offsets for each band linear regression ( Figure 11) agree with the overall results shown in Figure 10 above. While QAA TM showed better results for the band at 560 followed by 665 and 492 nm respectively, the QAAv6 showed better results for the band at 665, followed by 560 nm.
Remote Sens. 2020, 12, x FOR PEER REVIEW 15 of 23 ( Figure 11b) were worse when compared to those for QAATM, with R 2 values of 0.27, 0.79, 0.82, and 0.35; MAPE of 16, 23, 19, and 12%; and RMSE of 0.16, 0.14, 0.08, and 0.11 m −1 , for bands B1 to B4 respectively. These statistical results, slopes and offsets for each band linear regression ( Figure 11) agree with the overall results shown in Figure 10 above. While QAATM showed better results for the band at 560 followed by 665 and 492 nm respectively, the QAAv6 showed better results for the band at 665, followed by 560 nm.

Glint Correction
After the glint effect removal from band B3 to B5 resulted in more similar to those of the in situ (λ). Regarding bands B1 and B2, however, despite the improvement, Sentinel-2 (λ) were overestimated ( Figure 12). Even after the glint removal most of the bands remained with significant errors in its values when compared with in situ measurements. The bands B3 (560 nm) and B4 (665 nm) achieved the best results after the glint removal, with R 2 of 0.90 and 0.82 and MAPE of 11 and 33%, respectively. While the results for bands B1 (443 nm), B2 (492 nm) and B5 (704 nm) achieved the worst results, with MAPE higher than 58% ( Figure 12).

Glint Correction
After the glint effect removal from band B3 to B5 resulted in R rs more similar to those of the in situ R rs (λ). Regarding bands B1 and B2, however, despite the improvement, Sentinel-2 R rs (λ) were overestimated ( Figure 12). Even after the glint removal most of the bands remained with significant errors in its values when compared with in situ measurements. The bands B3 (560 nm) and B4 (665 nm) achieved the best results after the glint removal, with R 2 of 0.90 and 0.82 and MAPE of 11 and 33%, respectively. While the results for bands B1 (443 nm), B2 (492 nm) and B5 (704 nm) achieved the worst results, with MAPE higher than 58% (Figure 12). , and 12%; and RMSE of 0.16, 0.14, 0.08, and 0.11 m −1 , for bands B1 to B4 respectively. These statistical results, slopes and offsets for each band linear regression ( Figure 11) agree with the overall results shown in Figure 10 above. While QAATM showed better results for the band at 560 followed by 665 and 492 nm respectively, the QAAv6 showed better results for the band at 665, followed by 560 nm.

Glint Correction
After the glint effect removal from band B3 to B5 resulted in more similar to those of the in situ (λ). Regarding bands B1 and B2, however, despite the improvement, Sentinel-2 (λ) were overestimated ( Figure 12). Even after the glint removal most of the bands remained with significant errors in its values when compared with in situ measurements. The bands B3 (560 nm) and B4 (665 nm) achieved the best results after the glint removal, with R 2 of 0.90 and 0.82 and MAPE of 11 and 33%, respectively. While the results for bands B1 (443 nm), B2 (492 nm) and B5 (704 nm) achieved the worst results, with MAPE higher than 58% ( Figure 12).

Sentinel -2 MSI (λ) Time-Series
(λ) had the highest accuracies when using (λ) and (λ) derived from QAATM as input. Moreover, only the bands from B1 to B3 produced good statistical results, indicating that the approach combining QAATM plus (λ) should be applied only to these bands regarding the time series production. Figure 14 shows the monthly-averaged variation of (λ) for B2, the reservoir monthly-averaged inflow rate, and monthly-averaged outflow rate at the dam.   Because there were images only from the 2019 campaign, the MSI image results were only validated for the stations measured that year.
The supplementary material of this paper includes the summary for the K d−SA (λ) derived from MSI bands R rs (λ) validations for QAAv6 ( Figure S1).

Sentinel -2 MSI K d−SA (λ) Time-Series
K d−SA (λ) had the highest accuracies when using b b (λ) and a(λ) derived from QAATM as input. Moreover, only the bands from B1 to B3 produced good statistical results, indicating that the approach combining QAATM plus K d−SA (λ) should be applied only to these bands regarding the time series production. Figure 14 shows the monthly-averaged variation of K d−SA (λ) for B2, the reservoir monthly-averaged inflow rate, and monthly-averaged outflow rate at the dam. K d−SA (λ) refers to the Paraopeba River mouth, São Francisco River mouth, reservoir upper, middle and lower stream.

Sentinel -2 MSI (λ) Time-Series
(λ) had the highest accuracies when using (λ) and (λ) derived from QAATM as input. Moreover, only the bands from B1 to B3 produced good statistical results, indicating that the approach combining QAATM plus (λ) should be applied only to these bands regarding the time series production. Figure 14 shows the monthly-averaged variation of (λ) for B2, the reservoir monthly-averaged inflow rate, and monthly-averaged outflow rate at the dam.   Generally K d−SA (492) highest values were found for the rainy period of the year (Figure 15), coinciding with the months with higher accumulated precipitations and higher monthly-averaged flow rates (October to March) [30].
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 23 Generally (492) highest values were found for the rainy period of the year (Figure 15), coinciding with the months with higher accumulated precipitations and higher monthly-averaged flow rates (October to March) [30].   Representative K d−SA (492) results for two consecutive rainy and dry seasons are shown in Figure 16. During the dry season (Figure 16b,d), the K d−SA (492) is more spatially homogeneous, while for the rainy season (Figure 16a,c), higher gradients can be found, mostly at the entrance of the tributaries and in the upper portion of the reservoir. Generally (492) highest values were found for the rainy period of the year (Figure 15), coinciding with the months with higher accumulated precipitations and higher monthly-averaged flow rates (October to March) [30].

Algorithms Results Comparison
The K d−SA (λ) retrieved using IOPs derived from QAA TM as input achieved better performance compared when using IOPs derived from QAAv6 for both in situ ( Figure 11) and MSI ( Figure 13) R rs (λ) applicaiton. The linear regression at 492 nm for in situ data achieved slope factors of 1.21 and 0.96 for the QAA TM and QAAv6 approaches respectively. Besides that, the results in term of MAPE and RMSE were better for the QAA TM approach, achieving MAPE of 15%, and RMSE of 0.11 m −1 (Figure 11b) against MAPE of 23%, and RMSE of 0.14 m −1 for K d−SA (λ) retrieved with QAAv6, (Figure 11f) indicating better K d−SA (λ) accuracy when using QAA TM approach.
For the MSI application, QAA TM approach also achieved better results when compared to QAv6. The K d−SA (λ) derived when using MSI bands as input to retrieve IOPs from QAA TM showed higher accuracies for bands B1 and B2 (Figure 13a These results show that the K d−SA (λ) [25,26] reparametrization combining the QAAv6 and Liu et al. [57] methodology to derive IOPs is a suitable approach for waters optically similar to those of the TM reservoir.
In terms of comparison with the known K d (490) empirical algorithms found in the literature, the K d−SA (λ) algorithm using the QAA TM approach had good results. Table 6 shows the statistical results for 4 K d (490) empirical algorithms developed for different types of water. The MAPE of 14% and RMSE of 0.10 m −1 achieved by the K d−SA (λ) algorithm proposed in this study for the MSI data corroborate with the statement that this algorithm had good results.

Temporal Variation of the Diffuse Attenuation Coefficient of the Downwelling Irradiance
The highest K d−SA (492) values occurred during the rainy seasons, displaying peaks of almost 12 and 6 m −1 at its entrances for December 2017 and 2018, respectively ( Figure 14). The K d−SA (492) spatial distribution ( Figure 16) shows that São Francisco and Paraopeba River mouths display the higher values, followed by the TM upstream region, just before the two rivers confluence.
The raise of K d−SA (492) can be explained by increases in the turbidity [17,68], the chla [18,19], or the stratification caused by heat transfer in the upper layer of the water column [20].
Regarding the Brumadinho disaster, it was not possible to infer if it has affected the TM reservoir water transparency, once it occurred during the 2019-2020 rainy period (November-February), which was not included in the time-series of this study.

Limitations of the Algorithms
The approach applied in this study uses a combination of two semi-analytical algorithms parametrized for TM reservoir. Besides the fact that the algorithms used are semi-analytical, the QAA has some empirical steps that were reparametrized using data acquired in a single date. However, well parametrized semi-analytical algorithms usually do not lose accuracy when applied with input data from other periods of the year.

Semi-Analytical Algorithm Reparametrization
The values for the total a(λ) measured with the ACS were gathered for 2013 campaign, so the reparametrization of QAA TM empirical steps proposed in this study used only 2013 data, while validation used 2019 data.
The results of the reparametrization QAA TM step 2 (a(λ 0 )- Figure 9a) were influenced by two stations with higher values for the a(λ 0 ). These two stations were those with the higher values for all in situ limnology and radiometric data presented in this study (P25 and P26), as presented in the Sections 3.1 and 3.2. These stations were positioned near to the main tributaries entrances (São Francisco River and São Vicente Brook). These measurements are relevant since they represent the upper portion of the reservoir.
The reparametrization of QAA TM step 4 ( Figure 9b) is dependent on the remote sensing reflectance and total absorption coefficient parameters, as they are used for b bp estimates and can represent a source of uncertainty in the estimation. Despite that, this approach was used by other authors [57] and the accuracy are reliable for retrieving b bp values when no in-situ measurements are available. Further, the results of K d−SA (λ) proposed in this study where most affected at wavelengths with higher values (443 and 704 nm), those were MSI bands wavelengths (B1 and B5) that are not so relevant for clear water studies. The band B1 is an aerosol band used for quality assessment with a lower spatial resolution. While band B5 is a band at a wavelength where the total absorption in the water column is almost totally composed by the absorption of the pure water ( Figure 7). Furthermore, the MSI application for 2019 data was carried out for a single date in the middle of the campaign. The campaign occurred from 30 of June to 2 of July, with the image acquired for the 1st of July, which is not a significant limitation of the algorithm, since it did not rain during the campaign. However, the temporal resolution of MSI for Três Marias tile (5 days) and the field operation restrictions do not allow for the use of more images.

Conclusions
This study provides a bio-optical characterization for the TM reservoir for two dates during the 2013 and 2019 dry seasons. This characterization helped us to parametrize and apply the K d−SA (λ) for the water column. It followed international protocols for tropical aquatic systems data collection and processing, representing an important contribution for inland waters remote sensing applications.
The parametrized algorithm achieved more accurate results for K d−SA (λ) derived from field measurement when compared to K d−SA (λ) derived from the MSI image. Besides that, this approach showed good results for most of the MSI bands, making systematic monitoring of the water quality in the reservoir possible, and allowing a 5-day analysis.
To establish a very good spectrally accurate semi-analytical algorithm to estimate K d−SA (λ) for TM reservoir waters we suggest: • Bio optical data collection during rainy periods of the year to reparametrize the algorithm with more variated data; • Establishment of new algorithms for the correlation between OACs and optical properties of the water to better understand the primary production dynamics of the reservoir; • Assessment of other atmospheric correction and glint removal methodologies; • Apply the K d−SA (λ) plus QAA TM approach to a time-series including 2020 rainy period to better understand the Brumadinho disaster effects; • Apply the K d−SA (λ) plus QAA TM approach to other small reservoirs just upstream of TM in the Paraopeba River course, assessing the possibility of the retention of sediments coming from Brumadinho into these reservoirs.