Glint Removal Assessment to Estimate the Remote Sensing Reflectance in Inland Waters with Widely Differing Optical Properties

The quality control of remote sensing reflectance (Rrs) is a challenging task in remote sensing applications, mainly in the retrieval of accurate in situ measurements carried out in optically complex aquatic systems. One of the main challenges is related to glint effect into the in situ measurements. Our study evaluates four different methods to reduce the glint effect from the Rrs spectra collected in cascade reservoirs with widely differing optical properties. The first (i) method adopts a constant coefficient for skylight correction (ρ) for any geometry viewing of in situ measurements and wind speed lower than 5 m·s−1; (ii) the second uses a look-up-table with variable ρ values accordingly to viewing geometry acquisition and wind speed; (iii) the third method is based on hyperspectral optimization to produce a spectral glint correction, and (iv) computes ρ as a function of wind speed. The glint effect corrected Rrs spectra were assessed using HydroLight simulations. The results showed that using the glint correction with spectral ρ achieved the lowest errors, however, in a Colored Dissolved Organic Matter (CDOM) dominated environment with no remarkable chlorophyll-a concentrations, the best method was the second. Besides, the results with spectral glint correction reduced almost 30% of errors.


Introduction
The quality assessment of remote sensing reflectances (R rs ) is an essential task due to the crucial role of R rs products for water quality monitoring [1][2][3].In aquatic systems, targets that drive less than 10% of the radiometric quantities are registered by remote sensors [4].Radiometric quantities to compute R rs should be carefully measured and post processing must be used to retrieve the most reliable optical information [5].The R rs is defined as the ratio of water leaving radiance (L w , W•m −2 •sr −1 •nm −1 ) and the downwelling irradiance (E d , W•m −2 •nm −1 ) that reaches the water surface (z = 0+) for a specific viewing geometry (θ, φ that represents zenithal and azimuthal angles, respectively).Likewise, the R rs can be defined as the ratio of backscattering (b b ) to the b b plus absorption (a), which are inherent optical properties (IOPs) [6].Then, R rs is mathematically described as Equation (1) [7].
Relationships between all radiometric quantities are given by L t (θ, ϕ, λ) = L r (θ, ϕ, λ) + L w (θ, ϕ, λ).Considering L r as the radiance reflected by the water surface straight into the FOV (field-of-view) of the sensor and it is largely derived from L sky into the water [7], L r can be expressed as L r (θ, ϕ, λ) = ρ × L sky (θ, ϕ, λ).Therefore, we can re-write R rs as follow: where ρ (Equation ( 2)) (defined as the ratio of the surface-reflected radiance at the specular direction corresponding to the downwelling sky radiance from one direction) is the water-surface reflected factor used to correct the surface radiance reflection or glint effect-one of the main constraints of above-water measurements [9].The ρ values depend on the viewing geometry of the sensor, environmental conditions (incident irradiance, solar zenith angle-SZA, wind speed), and atmospheric composition (aerosol, molecules, cloud coverage) [7,10,11].Besides, the more accurate estimate of ρ, the more accurate are the R rs datasets for atmospheric correction validations, vicarious calibration of sensors, and bio-optical models [5,12].The glint component did not provide information from an underwater light field, then, it might be removed from R rs to the retrieve optical information about the optical components within a water column [9,13,14].There are some strategies to minimize the surface reflected light, such as: (i) to take measurements below the water surface (upwelling radiance, L u ) and then extrapolate it to above water [15,16]; and (ii) to plug a black cone into the sensor's head to protect the readings from the surface reflection (Skylight-blocked approach-SBA in [17]).Both strategies lead to particular constraints in regards to calculus approximations and the self-shadowing of the cone, respectively.
Regarding the above-water measurements, methods to estimate the ρ values are capable of removing the influence of L sky in the R rs measurements.Simulations of a radiative transfer numerical model with several geometrical views, meteorological conditions, and IOPs showed that for wind speeds lower than 5 m•s −1 , clear sky conditions and specific viewing geometry (θ v , ϕ v = 40 • , 90 • ), it is possible to assume a constant value of ρ = 0.028.However, such an assumption is dependent on the input parameters of the simulation and it did not take into account the spectral influence [7] or other cloudy conditions.Improved simulations considered the polarization ray effects and the water surface roughness [8], however, the spectral dependence was not included in those runs.
In respect to the simulations and field measurements, other glint removal methods includes the wind speed variations and part cloud cover, resulting in ρ = 0.0256 for overcast skies, and for clear to partly cloudy skies the values of ρ and wind speed relations can cause more than two orders of variations [9].Then, the computations of ρ can be made using a proper equation where ρ values are described as a function of wind speed.However, this does not take into account spectral dependence and also suggests an offset correction for residual glint (ε) that can produce negative spectra for overcast sky conditions.
Facing the challenge of the spectral dependence of ρ, Lee et al. (2010) proposed a method that relied on spectral optimization [5].The method was established considering two reflectances: the total remote sensing reflectance (T rs = L t /E d ), and sky remote sensing reflectance (S rs = L sky /E d ).The initial guess of R rs was done by using constant values of ρ, T rs , S rs , and a delta value (∆, average of R rs spectra from 750 to 800 nm) corresponding to a glint correction offset.Then, the initial guess of R rs is optimized by using the hyperspectral optimization method (HOPE, [18,19]).HOPE uses the IOPs (a and b b coefficients) to derive a modeled R rs .The HOPE outputs are used to compute the radiance reflectance just below the surface (r rs ), and it is extrapolated to the above-water surface reflectance (modeled R rs ).The difference between R rs measured (using T rs and S rs ) and R rs modeled (using HOPE outputs) is considered as a bias that is adjusted until reaches the minimal value for producing the final value of ∆, its means, the best correction of glint, and most accurate R rs .
The glint removal methods aim to remove, or at least minimize, the surface-reflected reflectance from R rs .However, all aforementioned methods based their corrections into the simulations of radiative transfer models or specific spectral features of R rs .Simulations of radiative transfer models require IOPs as inputs [7][8][9].In the case of spectral corrections of glint in Lee et al. (2010), absorption and scattering values are required to optimize the R rs spectra and establish ρ and ∆ [5].This method also uses L t , which is the registered radiance from the aquatic system after the incident irradiance penetrates into the water column, interacts with optical significant components (OSC) and water itself, and then re-emerges, passing through the surface until it reaches the sensor's view (sensors are pointed to the water surface to measure L t ).The IOPs are intrinsically related to the concentration and composition of OSCs [20], and L t is the total radiance that reaches the sensor (partially due to surface reflectance and partially from the energy of water column after interacts with OSCs and water itself) [7].The optically complex waters may affect the glint correction results, and consequently the R rs accuracy, because the main steps of glint removal methods are dependent of L t and IOPs that are directly related to concentrations of OSCs within the water column.
The present work aims, therefore, to assess the methods to remove the glint effect in inland water systems with widely differing optical properties.While we sought to trace the glint effect, comparisons between in situ R rs and R rs orbital data that resulted from the Landsat surface reflectance code (LaRSC) of the Operational Land Imager sensor (OLI/Landsat-8) are also outlined, aiming to investigate the reliability of the glint removal method to provide valuable information for quality control of image data.

Tietê River Cascade System of Reservoirs (TCSR)
The Tietê River, located in Brazilian Southeast (Figure 1), contains six reservoirs in cascade-Barra Bonita (BB), Bariri (BAR), Ibitinga (IBI), Promissão, Nova Avanhandava (NAV), and Três Irmãos from up to downstream.The BB, BAR, IBI and NAV reservoirs were chosen for study because they are responsible for producing more than 90% of the hydroelectrical energy from the entire cascade.BAR, IBI, and NAV are run-of-river reservoirs, whereas BB is an accumulation reservoir.(modeled Rrs).The difference between Rrs measured (using Trs and Srs) and Rrs modeled (using HOPE outputs) is considered as a bias that is adjusted until reaches the minimal value for producing the final value of ∆, its means, the best correction of glint, and most accurate Rrs.
The glint removal methods aim to remove, or at least minimize, the surface-reflected reflectance from Rrs.However, all aforementioned methods based their corrections into the simulations of radiative transfer models or specific spectral features of Rrs.Simulations of radiative transfer models require IOPs as inputs [7][8][9].In the case of spectral corrections of glint in Lee et al. (2010), absorption and scattering values are required to optimize the Rrs spectra and establish and Δ [5].This method also uses Lt, which is the registered radiance from the aquatic system after the incident irradiance penetrates into the water column, interacts with optical significant components (OSC) and water itself, and then re-emerges, passing through the surface until it reaches the sensor's view (sensors are pointed to the water surface to measure Lt).The IOPs are intrinsically related to the concentration and composition of OSCs [20], and Lt is the total radiance that reaches the sensor (partially due to surface reflectance and partially from the energy of water column after interacts with OSCs and water itself) [7].The optically complex waters may affect the glint correction results, and consequently the Rrs accuracy, because the main steps of glint removal methods are dependent of Lt and IOPs that are directly related to concentrations of OSCs within the water column.
The present work aims, therefore, to assess the methods to remove the glint effect in inland water systems with widely differing optical properties.While we sought to trace the glint effect, comparisons between in situ Rrs and Rrs orbital data that resulted from the Landsat surface reflectance code (LaRSC) of the Operational Land Imager sensor (OLI/Landsat-8) are also outlined, aiming to investigate the reliability of the glint removal method to provide valuable information for quality control of image data.

Tietê River Cascade System of Reservoirs (TCSR)
The Tietê River, located in Brazilian Southeast (Figure 1), contains six reservoirs in cascade-Barra Bonita (BB), Bariri (BAR), Ibitinga (IBI), Promissão, Nova Avanhandava (NAV), and Três Irmãos from up to downstream.The BB, BAR, IBI and NAV reservoirs were chosen for study because they are responsible for producing more than 90% of the hydroelectrical energy from the entire cascade.BAR, IBI, and NAV are run-of-river reservoirs, whereas BB is an accumulation reservoir.  , respectively.The time residence of water ranged from 25-41 days in IBI, and 46 days on average in NAV [25].IBI is considered a meso-to-eutrophic environment [26], and NAV is considered an oligo-to-mesotrophic waters [27].
Two field campaigns were conducted in each reservoir (see details in Table 1), and the sample geographic locations were selected following the design plan steps described by Rodrigues et al. (2016) [28].Measurements of radiometric quantities and water sample collection were taken at each sampling spot.The following water quality parameters were also collected: Secchi Disk Depth (Z SD , m), suspended particulate matter (SPM, mg•L −1 ) and its respective organic and inorganic fractions (POM and PIM, mg•L −1 ), chlorophyll-a (Chl-a, mg•m −3 ), wind speed (m•s −1 ) and depth (m).Turbidity (NTU), wind speed (m•s −1 ), Z SD (m), and bottom depths (Z, m) were taken in triplicate to reduce systematic errors.Water samples, in duplicate, were collected in each spot, and were filtrated at low vacuum by using Whatman GF/F filters (47 mm of diameter and 0.7 µm pore size), to establish the SPM, PIM, POM, and Chl-a concentrations.The SPM concentrations, and their organic and inorganic fractions, were obtained following the APHA protocol (1998) [29].The Chl-a concentrations were obtained by the acetone extraction method and spectrophotometer analysis [30].Replicates were made in order to avoid systematic errors.

Optical Data
Optical measurements were made in situ, to gather radiance and irradiance measurements to compute R rs and apply glint correction methods.In the laboratory the water samples collected at each sampling location established the absorption coefficients from the phytoplankton, non-algae particles, and colored dissolved organic matter, in addition to the total absorption coefficient.
The in situ measurements were made between 400 and 900 nm using TriOS Hyperspectral radiometers (RAMSES TriOS ® ; TriOS, Germany).L t and L sky , both in W•m −2 •sr −1 •nm −1 , were measured with ARC sensors that contained a 7 • field of view, and 3.3 nm of the sampling interval.The downwelling irradiance incident onto the water surface (E s , in W•m −2 •nm −1 ) and the downwelling irradiance below the water surface (E d , in W•m −2 •nm −1 ) were measured by an ACC sensor equipped with a cosine collector for registering the light within the hemisphere.Errors due to slight wavelength differences were avoided by linearly interpolating all the measurements at 1 nm spectral resolution before computing R rs .
All field campaigns assumed the same geometric view and were set up to minimize the sun glint effect.Measurements of L sky and L t were made using the zenith solar angle (θ v ), which equaled 40 • from zenith and the sensors were pointed towards nadir.An azimuthal angle (ϕ v ) of 90 • (solar plan was considered the origin) was also adopted to avoid boat shadows and to minimize the specular radiance effects [7,31].

Inherent Optical Properties
Laboratory measurements were made to obtain the OSC's absorption features.The particulate absorption (a p = a nap + a phy , sum of non-algae particle and phytoplankton absorptions) was computed by using the Trasmittance-Reflectance (T-R) method [32,33].A know volume of sampled water was filtered through Whatman GF/F filters (47 mm of diameter and 0.7 µm of pore size) at low vacuum.The filters were submitted to spectral readings of absorbance taken by a dual-beam UV-2600 UV-VIS spectrophotometer (SHIMADZU, Japan) with 1 nm of spectral resolution, ranging from 280 to 800 nm and using a blank filter dried with Mili-q water drops as a reference.The absorbance readings were converted into a p , and the filters were bleached of the organic fraction of the particulate by using hypoclorite solution (NaClO at 10%).After removing the solution with died pigments, the bleached filter pads were re-submitted into new readings of absorbance to retrieve a nap .Then, a phy was obtained by subtracting a nap from a p .
The Colored Dissolved Organic Matter (CDOM) absorption coefficient (a cdom ) was obtained by filtrating the water samples through a nylon membrane with a diameter of 47 mm and 0.22 µm pore size.The filtrate was placed in quartz cuvettes with an optical path of 10 cm (L = 0.1 m).Then, the spectral absorbance readings were taken to derived a cdom (λ) after the filtrated being acclimates to room temperature.Besides, a baseline correction was applied by subtracting residual absorbance (mean absorption values between 700-750 nm) [34].

Remote Sensing Reflectance
The water color radiometry can be defined in terms of R rs , which is strongly affected by the surface-reflected component since L r is affected by the meteorological conditions, cloud coverage, and wind speed.Each facet of a wave is capable of reflecting the light into the sensor view [5].In order to correct the glint, the factor acts as the probability factor to estimate the proportion of energy from the sky that will reach the water surface and reflect into the sensor.The sun glint effect can be disregarded due to the low latitudes [35].
Four approaches of glint removal of R rs were evaluated to understand the impacts of "ρ factor" in different optical compositions.The first approach, hereafter M1, adopted ρ = 0.028 for all sampling spots [7]; the second approach (M2) considered a ρ factor that changes according to the viewing geometry, acquisition time, and illumination geometry [8].The third approach (M3) used a ρ factor that is spectrally dependent [5]; and the fourth approach (M4) used the ρ factor as a function of wind speed [9].
The M1 approach adopted ρ = 0.028 because it represents the maximum value of L sky that will be reflected into the sensor.When the wind speed in the field was not higher than 5 m•s −1 , the geometry view was (θ v ,ϕ v ) = 40 • , and 90 • was maintained, the measurements were made under clear sky conditions.In addition, it is recommended that measurements be taken during the solar hour between 10 a.m.-2 p.m. [7].
The M2 approach was calculated based on a the LUT (Look-Up-Table ) of ρ values, which requires as inputs the wind speed (WS, in m•s −1 ), sun zenithal angle (SZA), and viewing geometry adopted from the sensors represented by ϕ v and θ v [8].The LUT of ρ results from several HydroLight simulations with several wind speeds, and geometries; although it was not analyzed, the values of the simulations are expressed in terms of quads that are greater than the areas comprised by the sensor's FOV.The differences between M1 [7] and M2 [8] are the methods used to simulate the water surface wave (Cox-Munk or Fourier transformation) and consider the polarization of light.Moreover, M2 takes into account the wave shadowing effect.The inputs for LUT were wind speed measured in situ by using a portable anemometer, the viewing geometry adopted as θ v , ϕ v = 40 • , 90 • , and the SZA calculated by Equation ( 3): where δ is latitude from each sampling spot collected by Garmin GPS, γ is the solar declination calculated via the NOAA website (http://www.esrl.noaa.gov/gmd/grad/solcalc/),H is the sidereal hour retrieved by the difference between the sidereal hour at 0 h in the UTC system and the right ascension (α).H and α values were available from the Astronomic Annularly of the Brazilian National Observatory (2014, 2016, and 2017) and the values were interpolated considering the time zone for Brazil in each sampling spot.The values in LUT were used to compute R rs using Equation ( 2).The M3 approach considers a spectral ρ variation [5].As a premise, there is an initial guess of the reflectance value calculated by Equation ( 4): where R Fresnel is the radiance reflectance computed with ρ = 0.022 and ∆ is the R mean(750-800nm) , i.e., an offset calculated as the mean of R Fresnel between 750 and 800 nm, that supposed to be set as zero for Case-1 waters and extended to longer wavelength intervals if the measurements longer than 800 nm were made for Case-2 waters.The step forward consists of using the hyperspectral optimization algorithm (HOPE) [17,36] for solving the inversion problem and modeling the total scattering (b), total absorption coefficients (a t ), and specific coefficients (a w , a dg , a phy ).These values of b, a t , a w , a dg , a phy , are used as inputs to compute the subsurface reflectance (r rs ) that is converted into remote sensing reflectance R rs using Equation (5).
The r rs values are modeled on the HOPE solution.Considering the R rs and R guess , the called bias [5] represents the surface correction value, and can be computed as an objective function (Equation ( 6)).
With R guess calculated from (4) and R rs is calculated from (5), 675 400 R guess − R rs 2 is the square of mean reflectance from 400 to 675 nm, and 800 750 R guess − R rs 2 is the square of mean reflectance from 750 to 800nm.The Error ( 6) is computed for each adjustment until reaches the minimal value for producing the final value of ∆, initially assumed as ∆ = R mean (750-800 nm).
It is important to highlight that if we used the assumptions of HOPE, which uses a 0 and a 1 coefficients to model the phytoplankton absorption coefficient, the interactions did not achieve the minimum value of Error, leading to high values of delta and the shape of spectra was negative.Then, we adjusted the method to use specific absorption coefficients (a * phy ) instead of the modeled coefficient, then the values are closer to the real measurements and the interactions are converted into a minimum Error, resulting in corrected R rs spectra, also named optimize reflectance values (R opt -Equation ( 7)).
where the terms in Equation ( 7) are the total and sky radiances (L t and L sky ), downwelling irradiance (E d ) and proportional coefficient of the sky reflected light (ρ).∆ is related to a residual glint correction, which means sun glint, whitecaps, and foam reflected light [5].The M4 approach relies on sun sky conditions as well as features at 750 nm, where L sky (750)/E d (750) were lower than 0.05 (from clear to partially cloud covered sky), the value of ρ is calculated as a function of wind speed (W, m•s −1 ), as described in Equation ( 8), if the ratio is higher than 0.05, then ρ = 0.0256 and represented by an overcast sky.The main assumption is that reflectance is dominated by pure water absorption in near infrared (700-900 nm), therefore, a registered signal in this range results from glint (assuming no other sources of errors) [9].
where W is the wind speed (m•s −1 ) in each sampling location, and the coefficients in Equation ( 8) were obtained according to radiative transfer model simulations [7].Then, the values retrieved from Equation ( 8) were used to compute R rs (Equation ( 2)).
The benefits and disadvantages from each glint removal method used in this work are accessed in Table 2.

M1 [7]
• Minimizes the skylight effects using Cox-Munk simulations for waves and unpolarized ray • Spectral dependence of skylight and wind speed variations are not considered.

•
Simulations were developed for sea optical compositions.

M2 [8]
• Minimizes the skylight effects using wave spectral variance and considered the ray-polarization.
• Spectral dependence of skylight is not considered.

•
Desires additional information collected in fieldwork, such as wind speed and SZA computations.

M3 [5]
• Minimizes the sky and sunlight effects by using HOPE-semi-analytical scheme for retrieving IOPs.
• Needs additional information collected in fieldwork, such as the absorption and scattering coefficients.

M4 [9]
• Minimizes the skylight effect accordingly to the wind speed and was developed for turbid environments.

•
Requires additional information about wind speed quantities.

HydroLight Simulations
HydroLight (version 5.2; Sequoia Scientific, Inc., Bellevue, WA, USA; [37]) is a radiative transfer numerical model capable of computing the radiance distribution, which relies on the inputs (IOPs and the OSCs concentrations) and contour, such as the nature of the wind-blown water surface and the sky conditions.HydroLight produces an in-water light field, with outputs split into radiance components that allows us to compute the reference dataset of R rs from 400-800 nm via Equation (1), removing the unwanted contributions, it means surface reflected light.
The most representative R rs curves from each fieldwork were selected by applying the k-means, and then clustering results were compared to the sample curves by using a Spectral Angle Mapper (SAM [38]).Overall, a total of 24 of 175 R rs curves (at least 2 curves from each field) were selected and 24 sampling spots with different optical features were used in HydroLight simulations.
HydroLight experiments were performed using a Case-2 model, where it is required to define concentrations and specific absorptions from four OSCs (pure water, chlorophyll-bearing particles, CDOM, and mineral particle).We used the default IOPs of water that were relied on by Pope and Fry (1997) for absorption [39], and Smith and Baker (1981) for scattering [40].Specific coefficients of Chl-a and mineral particles were user-supplied (files with specific absorptions between 350-800 nm were set into the HydroLight), as well as the concentrations (in mg•m −3 and g•m −3 ) that were considered as constant with depth (we did not collect measurements along the water column).The CDOM component required the value of Slope (S cdom ) and the absorption at λ 0 , considered at 440 nm.
The scattering model and phase function were kept as the default setup that represents the GAM model parameters and Petzold's average-particle phase function [6], respectively.Some sampling stations were simulated, including chlorophyll fluorescence.The inelastic scattering was not included in the runs.Besides, the geographical location used to compute the sun zenith angles, wind speed, water temperature, and acquisition dates was informed as an input of each sampling spot.For all simulations, we assumed a clear sky in RADTRANX, by using cloud cover equal to 0, and no atmospheric parameters were modified because we did not take measurements from the atmospheric composition.Infinite depth (no bottom contributions) and a spectral range of 400-800 nm (∆λ = 1 nm) were defined.All inputs are summarized in Table 3.The printout returned information about all of the parameters and inputs used during the runs.The outputs, written into ASCII files, were formatted to be read by spreadsheets in Excel (macros).Considering this option, we extracted the R rs dataset from the outputs of HydroLight.Then, reference HydroLight R rs (reference dataset) was compared to the in situ R rs corrected by the four methods aforementioned.

Accuracy Assessment of R rs
Considering several approaches to evaluate accuracy on R rs , a better understanding is provided when we evaluated the optical closure between in situ measurements to the simulated data [41]; therefore, comparisons between R rs calculated by M1, M2, M3, and M4 to the simulated R rs assessed which method retrieved the most reliable optical information.The matching was evaluated using statistical errors-Root Mean Squared Difference, RMSD in sr −1 , normalized Root Mean Squared Difference, nRMSD in %, and Mean Absolute Percentage Error (MAPE, in %) (Equations ( 9)-( 11)).
× 100 (10) where x e,i and x r,i are the values estimated and measured R rs from M1, M2, or M3 and simulated R rs (reference data); max and min are maximum and minimum values.RMSD can be used for establishing random errors, such as environmental issues [42], whilst MAPE is used as a scatter factor of measurement, indicating systematic errors such as radiometric or platform errors [42].

OLI Assessment
OLI applications for inland water systems were recently investigated and demonstrated great results [43][44][45][46] due to their spatial resolutions that are able to represent significant details in hydrologic process in such environments.Besides, the Landsat dataset has demonstrated high suitability for inland and coastal studies [47,48].Assessment of matching between the best glint removal of R rs and the reflectance product of the Operational Land Imager (LaSRC) onboard Landsat-8 was performed.R rs -glint removal based methods were spectrally convoluted (Equation ( 12)) with the spectral response function (SFR) of each band of OLI centered at 593 nm, 482 nm, 561 nm, and 665 nm.The last band of OLI was not considered in our analysis due to the limitation of HydroLight simulations in the 400-800 nm range.
where R s rs is the convoluted in situ R rs for OLI bands (centered at 443 nm; 482 nm; 561 nm and 655 nm), SFR is the spectral function response for each band [49], and λ max and λ min are the longer and shorter wavelength in each OLI band.Then, the convoluted values of R rs were compared to the R rs from OLI images to assess if our glint removal dataset matches the OLI dataset.Considering that LaSRC retrieved the irradiance reflectance, the images were converted into R rs by dividing them to π, and applying the scale factor of 10,000 [46].The images were also visually assessed in order to identify the presence of glint.None of them presented visual glint effects in the water surface.
The matching between in situ samples and image samples was made considering only in situ measurements that were taken at the same day that OLI overpass, corresponding to 34 sampling locations.A summarized flowchart of methodology is shown in Figure 2.

Water Quality Data and Optical Characterization
A wide variability of water quality parameters was observed in TCSR's dataset.The Chl-a concentration ranged from 1.37 mg•m −3 to 797.8 mg•m −3 , SPM ranged from 0.10 mg•L −1 to 44 mg•L −1 (Table 4), and most of the samples presented a high composition of particulate organic matter

Water Quality Data and Optical Characterization
A wide variability of water quality parameters was observed in TCSR's dataset.The Chl-a concentration ranged from 1.37 mg•m −3 to 797.8 mg•m −3 , SPM ranged from 0.10 mg•L −1 to 44 mg•L −1 (Table 4), and most of the samples presented a high composition of particulate organic matter (POM), excluding the NAV's samples.Furthermore, Z SD comprised between 0.37 and 4.80 m, whilst turbidity comprised between 1.01 and 80.9 NTUs (Table 4), leading to the most variable water quality parameter (CV = 88.7%).Except BB, the other fieldworks reached wind speeds higher than 5 m•s −1 .The optical contributions from each OSC to the total absorption disregarding water itself (a t-w (λ)) were plotted in ternary diagrams (Figure 3) using the central wavelength of OLI bands.Each reservoir demonstrated specific OSCs dominance.Furthermore, the light absorption budgets for each reservoir did not change along the wavelengths, it means, the OSCs dominance were not spectrally dependent.

Rrs from M1, M2, M3, and M4 Approaches
Results from all glint removal methods showed a marked reduction of magnitudes on Rrs curves (Figure 4).Overall, the most variability was found in the NIR range (700-800 nm), where we observed that M3 showed the lowest magnitudes of CV (not valid for BB reservoir), it means that in the NIR range the Rrs_M3 provided the lowest variation among sampled curves, or rather, the most clustered curves in the NIR interval.Furthermore, NAV showed a marked reduction in all spectra Where a ϕ is related to phytoplankton absorption, a NAP to non-algae particle absorption and a CDOM is related to the absorption from colored dissolved organic matter.

R rs from M1, M2, M3, and M4 Approaches
Results from all glint removal methods showed a marked reduction of magnitudes on R rs curves (Figure 4).Overall, the most variability was found in the NIR range (700-800 nm), where we observed that M3 showed the lowest magnitudes of CV (not valid for BB reservoir), it means that in the NIR range the R rs_M3 provided the lowest variation among sampled curves, or rather, the most clustered curves in the NIR interval.Furthermore, NAV showed a marked reduction in all spectra magnitudes by M2 and M3 approaches (Figure 4n,o).The same occurred in BAR curves when we used the M3 glint removal method (Figure 4g).It is possible to observe that Rrs curves presented several features along the spectra that changed for each fieldwork.Regardless of approach, the spectral features are conserved along 400-800 nm, varying only in magnitude.In order to quantify such variations along the spectra, we computed the CV for each 20 nm of range and Figure 5 represent them in graphs.
Similar results were obtained from M1 and M2 approaches in terms of CV and standard deviations (SD) for almost all fieldworks (Figure 5).The highest variations of Rrs where found in BB fieldworks for M3 in NIR ranges (700-800 nm), and NAV2 (Figure 5a,b,h), which are the opposite results compared to the other reservoirs (BAR1, BAR2, IBI1, IBI2, NAV1) that presented M3 as the approach with minimums CVs.It is possible to observe that R rs curves presented several features along the spectra that changed for each fieldwork.Regardless of approach, the spectral features are conserved along 400-800 nm, varying only in magnitude.In order to quantify such variations along the spectra, we computed the CV for each 20 nm of range and Figure 5 represent them in graphs.
Similar results were obtained from M1 and M2 approaches in terms of CV and standard deviations (SD) for almost all fieldworks (Figure 5).The highest variations of R rs where found in BB fieldworks for M3 in NIR ranges (700-800 nm), and NAV2 (Figure 5a,b,h), which are the opposite results compared to the other reservoirs (BAR1, BAR2, IBI1, IBI2, NAV1) that presented M3 as the approach with minimums CVs.The assessment of glint methods was made comparing the achieved results to the HydroLight simulations.However, to simulate the entire dataset composed by 175 spectral curves would imply a lot of time consumption and computational processing.Aiming to optimize the simulations, a group of representative samples was selected by using k-means to organize the entire dataset in clusters.Then, the most similar in situ spectra selected with the clustered data by using SAM.A total of 24 spectra were selected for HydroLight simulations/comparisons.We take care to select at least three samplings from each field-one sample near to the highest concentration of SPM, the other one near to the lowest concentration, and another one near to the average.
The plots from specific central wavelengths of OLI are shown in Figure 6, as well as the error assessments of glint removed-Rrs are shown in Table 5.The OLI images were chosen because they were suitable for inland water studies due to their spatial and temporal resolutions [50].The assessment of glint methods was made comparing the achieved results to the HydroLight simulations.However, to simulate the entire dataset composed by 175 spectral curves would imply a lot of time consumption and computational processing.Aiming to optimize the simulations, a group of representative samples was selected by using k-means to organize the entire dataset in clusters.Then, the most similar in situ spectra selected with the clustered data by using SAM.A total of 24 spectra were selected for HydroLight simulations/comparisons.We take care to select at least three samplings from each field-one sample near to the highest concentration of SPM, the other one near to the lowest concentration, and another one near to the average.
The plots from specific central wavelengths of OLI are shown in Figure 6, as well as the error assessments of glint removed-R rs are shown in Table 5.The OLI images were chosen because they were suitable for inland water studies due to their spatial and temporal resolutions [50].Overall, the nRMSD ranged from 7.6% to 54.1% over the entire TCSR's dataset.Excepting BB1 and IBI1, the M3 approach retrieved the lowest errors (Table 5).It is important to highlight that IBI2 showed a meaningful variation of errors comparing all corrections, ranging from 16.6% to 40.2%.Meaningful differences between glint removal methods were also found in BB2-the errors comprising between 18.8% and 54.1%, reducing the error by about 35%.Regarding all glint removal methods tested in this work, the most suitable Rrs dataset was provided by M3, it means, the spectral ρ correction (Lee et al., 2010).Overall, the nRMSD ranged from 7.6% to 54.1% over the entire TCSR's dataset.Excepting BB1 and IBI1, the M3 approach retrieved the lowest errors (Table 5).It is important to highlight that IBI2 showed a meaningful variation of errors comparing all corrections, ranging from 16.6% to 40.2%.Meaningful differences between glint removal methods were also found in BB2-the errors comprising between 18.8% and 54.1%, reducing the error by about 35%.Regarding all glint removal methods tested in this work, the most suitable R rs dataset was provided by M3, it means, the spectral ρ correction (Lee et al., 2010).
Besides the analysis for each fieldwork, we evaluated the errors based on OLI central wavebands (Table 6).Therefore, glint-removed R rs and Hydrolight R rs were spectrally convoluted (Equation ( 12)) and comparisons were made by wavelength.Relying on the results, the lowest errors retrieved from M3 approach again (20.34-29.8% of nRMSD and 13.9-39.3% of MAPE), whereas the highest errors retrieved from M4 approach (30.1-64.2% of nRMSD and 26.2-86.6% of MAPE).Regarding the wavelengths, the highest errors were found in 443 nm, whereas the lowest errors were found in 665 nm.

OLI R rs Versus Glint Removal R rs_M3
Plots of matching pairs between the LaRSC product of OLI bands and the in situ R rs retrieved from the M3 glint removal method are shown in Figure 7.For such an analysis, we used solely the sampling spots (n = 34) that matched with the Landsat-8 overpass in order to minimize the temporal effects over the results.It is important to highlight that if the image for NAV1 presented an overcast sky, then, this fieldwork was not considered in our analysis.
Besides the analysis for each fieldwork, we evaluated the errors based on OLI central wavebands (Table 6).Therefore, glint-removed Rrs and Hydrolight Rrs were spectrally convoluted (Equation 12) and comparisons were made by wavelength.Relying on the results, the lowest errors retrieved from M3 approach again (20.34-29.8% of nRMSD and 13.9-39.3% of MAPE), whereas the highest errors retrieved from M4 approach (30.1-64.2% of nRMSD and 26.2-86.6% of MAPE).Regarding the wavelengths, the highest errors were found in 443 nm, whereas the lowest errors were found in 665 nm.

OLI Rrs Versus Glint Removal Rrs_M3
Plots of matching pairs between the LaRSC product of OLI bands and the in situ Rrs retrieved from the M3 glint removal method are shown in Figure 7.For such an analysis, we used solely the sampling spots (n = 34) that matched with the Landsat-8 overpass in order to minimize the temporal effects over the results.It is important to highlight that if the image for NAV1 presented an overcast sky, then, this fieldwork was not considered in our analysis.The OLI3 (561 nm) and OLI4 (655 nm) plotted in Figure 7 are well distributed along the 1:1 line and showed r = 0.60 (p-value < 0.05).OLI1 and OLI2 did not show relevant correlation coefficients (r = 0.002 and 0.05, respectively), which implied higher MAPE in OLI1 and OLI2, 43.6% and 29.6% when compared to OLI3 and OLI4-MAPE of 18.4% and 28%.

Discussion
Computing suitable Rrs spectra is not an easy task.The glint removal methods are widely used in a broader range of oceanic and freshwater systems; however, there is no consensus of which method is the optimum way to provide accurate datasets [5].The TCSR presented a wide range of The OLI3 (561 nm) and OLI4 (655 nm) plotted in Figure 7 are well distributed along the 1:1 line and showed r = 0.60 (p-value < 0.05).OLI1 and OLI2 did not show relevant correlation coefficients (r = 0.002 and 0.05, respectively), which implied higher MAPE in OLI1 and OLI2, 43.6% and 29.6% when compared to OLI3 and OLI4-MAPE of 18.4% and 28%.

Discussion
Computing suitable R rs spectra is not an easy task.The glint removal methods are widely used in a broader range of oceanic and freshwater systems; however, there is no consensus of which method is the optimum way to provide accurate datasets [5].The TCSR presented a wide range of water quality parameters, from turbid to clearer waters (turbidity ranged from 1.01-80.9NTU), with algae blooms, and scums.In some cases there is a CDOM dominance in some reservoirs with brown color, which was observed during in situ measurements (BAR2 and IBI2).
Water quality parameter variations from up-to-downstream inside a unique continuum system (Tietê River) are dependent on the time residence of the water, inputs from tributaries along the cascading reservoirs, auto-depuration, the availability of solar irradiation, and the activities developed within the sub-basin areas that contribute to wastewater discharges.Barra Bonita, for instance, is near to São Paulo (capital of São Paulo State) and it is the second most populous city in Brazil [51].The Bariri reservoir receives discharge from the river that flows throughout Bauru city, an urban area that recently adopted wastewater treatments, and Ibitinga is near Environmental Protection Areas, which contribute to the organic matter as a result of the vegetation coverage [26].Besides, the continuum river cascade concept (CRCC) in TCSR contributes to meaningful changes throughout a sequence of impoundments within the same river-retention of particles as turbidity decreases and a higher amount of energy along the water column in downstream accumulates [25].
Several water types can be found in TCSR identified by ternary diagrams (Figure 3).With absorptions displaced towards the CDOM axes in BAR and IBI, BB showed absorption coefficients displaced towards the phytoplankton axes, and NAV presented absorption coefficients well distributed inside the ternary diagram with slightly higher contributions of NAP [27].The wide variation in OSC's absorption affects the spectral features in the R rs curves (Figure 4).We observed absorption features of Chl-a at 440 nm and 675 nm in all reservoirs [52], and the phycocyanin absorption near 620 nm [51] can be realized in BB and BAR fieldworks (Figure 4a-d).The iron absorption features (450-550 nm) and high CDOM absorption coefficients (<500 nm) resulted in flat Rrs spectra at shorter wavelengths (<550 nm), mainly in curves obtained in BAR, IBI, and NAV fieldworks [53][54][55].In case of black waters, CDOM absorption can be higher than water absorption even at 700 nm [56].The increase in magnitudes at 550 nm is related to several SPM concentrations and compositions, with inorganic particles causing light scattering in the green spectral region [57,58].A minor reflectance feature near 760 nm, in almost all spectral curves, can be associated to oxygen absorption [14] and normally related to the glint effect.Computing the dip of oxygen absorption, the results from all approaches were close to zero (D ≤ 0.001), which means the glint effect was removed [59].
After applying the glint correction methods, the features in the R rs spectra only changed in magnitude.The main differences between glint removal methods are presented in Table 2.The adoption of specific geometry views aimed to rule out problems with sun glint, shadows, reflections, foam, and whitecaps [7].During in situ measurements, we adopted the sensor's viewing of 40 • from nadir and from zenith to measure L t and L sky as recommended by Mobley (1999) [7], nevertheless, the wind speed overcame 5 m•s −1 in BAR, IBI and NAV fieldworks, which makes the adoption of ρ = 0.028 unfeasible.Besides, the effects on several viewing geometry sets under radiometric measurements are described in [41].The wind speed variations caused roughened surfaces, where each wave facet caused the solar light to be reflected into the sensor's viewing and resulted in different L sky contributions [3].Facing such a challenge, the surface radiance needs to be modeled and removed from R rs by using the sea-surface reflectance factor-ρ.Then, ρ models the glint effects caused by the waves reflection of radiance.When better ρ is estimated, more accurate is R rs computations are achieved [5].
Accurate estimates of ρ values depend on the sky conditions, sea state, viewing geometry, and environmental conditions.In order to control the variable geometry viewings, our sensors were coupled on steel frames positioned in boat platforms, so the instability of the boat caused by waves was taken into account by the measurements that were not commonly taken in a perfect viewing geometry.Apart from geometry, some measurements were taken under unstable conditions beyond human control, such as the variability of environmental conditions (wind speed, cloud coverage variations, and atmospheric composition with particles due to sugarcane burned areas near reservoirs [60]).
The glint removal results demonstrated similar patterns between the M1 and M2 approaches for almost all reservoirs (NAV2 was the exception), with a very similar coefficient of variation of R rs (Figure 5).A previous study developed by Garaba et al. (2015) indicates that some glint removal methods are overlapped, resulting in similar corrections due to the same optical assumptions.68% of ρ values from M2 were comprised between 0.025 and 0.028, which demonstrated a relative overlap of ρ values with M1, and consequently, very similar results [2].The remainder values of M2 were equally divided into two main categories: 0.020-0.025;and 0.028-0.037.The highest value of ρ was also accomplished with wind speed near 5 m•s −1 .
M4 did not show similar results from M1, M2, and M3 approaches in environments where inorganic loads are an expressive component of SPM (around 40% on average).This conclusion can be obtained by observing IBI1 and NAV2 (Figure 5e,h) fieldworks, which can indicate that M4 showed different behavior when applied to environments with marked inorganic composition, although they were developed for specific turbid areas.
CDOM-rich environments (BAR2, IBI1 and IBI2) did not present wide variations in the R rs visible range (400-700 nm), on the contrary, the NIR range showed a higher variation as depicted in Figure 5 (CV of M4 in BAR2 ranged from 40% to 80%).Considering BAR1 as a CDOM environment with remarkable phytoplankton concentrations, the results did not follow the same pattern as BAR2, IBI1, or IBI2.BAR1 presented variability very similar to the BB1 and BB2, which were fieldworks with high SPM concentrations and high levels of Chl-a, with samples reaching more than 250 mg•L −1 .The M3 approach depicted the most distinguishable CVs among all of the tested approaches.The variations in M3 results can be a consequence from the ∆ corrections that took into account, as an initial guest, the average of R rs between 750-800 nm (in some cases, considering even longer wavelengths for such computation).This spectral region was influenced by a marked peak at 760 nm in some R rs spectra, then, the initial guess of ∆ can be highly variable, resulting in widely spectral optimizations.
Facing the challenges of glint removal, we simulated 24 spectra from all fieldworks using HydroLight to obtain a reference dataset.Then, we spectrally resampled the in situ and HydroLight R rs to compute the errors.Error assessments of each fieldwork retrieved MAPE ranging from 7.6% to 54.1%, with M4 retrieving the highest errors.The errors from M4 can be caused by the wind speed equation based on optical simulations (Equation ( 8)) that probably did not correctly estimate the spectral-reflectance proportional factor, and did not retrieve the most accurate R rs datasets.
We relied on the errors assessment demonstrated in Table 4. Almost all glint-removed methods did not work for CDOM-rich environments, such as IBI1 and IBI2, where the errors achieved were higher than 20%.The most suitable glint removal method was M2 (with errors around 8% and 17%) for IBI1 and IBI2.Regarding the M3 results for IBI2, we improved from 20.8% in M3 to 11% when we did not take into account the sampling with lowest Chl-a and SPM concentrations.This result indicates that the M3 approach is highly sensible for measurements carried out in environments with remarkable CDOM amounts and less expressive Chl-a/SPM concentrations.
Modeling glint effects from higher Chl-a concentrations, and consequently high SPM concentrations, was also a challenge.Removing the samples with the highest concentrations for BB1 and IBI1, the MAPE changes from 18.4% to 11.5%, and 16.5% to 9.0%, respectively, suggesting that M3 is a good approach overall the methods, however, it was very sensible when the dataset contained a few concentrations of outliers.Overall, the M3 approach showed the best results in BB2, BAR (both fields), and NAV (both fields).Considering that M3 approach corrects the sun and sky glint effects, and also takes into account the spectral variation of ρ (in our case, the values varied in a factor of 8 such as the results showed in [5], the best results were expected from the glint removal method.Overall, the M3 approach was the most suitable method for glint correction in high turbid waters, as found in BAR and BB, and clearest waters, as found in NAV.
Considering several studies about bio-optical modeling for inland waters [27,51,61] and the potential of OLI for its aquatic monitoring [62,63], we evaluated the matching between our glint corrected R rs and LaSRC from OLI.The OLI dataset has been widely applied to inland water systems due to its spatial resolution and the detailed richness of such environments [64][65][66][67].
In the relied error analysis (Table 6), it was observed that the coastal band retrieved errors near 43% of nRMSD (and 54% of MAPE), which can be related to atmospheric issues, such as high influence of Rayleigh scattering and aerosol concentrations [64].Besides, the green and red band retrieved lower errors (nRMSD < 26%) even with some sampling spots displaced towards higher values (R rs_M3 > 0.020 sr −1 ) for the green region (Figure 7).These points related to high levels of Chl-a concentrations (>250 mg•L −1 ).No wind and haze days meant the sky was partly clouded and this increased L sky , making the glint correction more challenging, which probably caused the mismatching.
OLI errors from Table 6 indicate the suitability of OLI bands in the green and red bands, however, the short wavelengths represented by the Coastal (443 nm) and Blue (482 nm) bands showed significant errors.In 443 nm, the MAPE was higher than the other wavelengths, which can be caused by two main problems: atmospheric issues in the blue range or high Chl-a concentrations that are additional to the higher CDOM absorptions.A deep decrease in the R rs in this range complicates atmospheric correction.It is worth mentioning that the spectral mixture and spatial resolution impacts the acquisition of images, and in situ measurements were not evaluated in this paper, yet they have implications for matching.Based on our results, and considering the OLI applications for inland environments, it will be possible to use OLI green and red bands for monitoring SPM.For instance, for CDOM estimates [68], if 443 nm is used for CDOM absorptions [69,70], it might be necessary to achieve better evaluations of the atmospheric correction methods, as well as, the impacts of glint-removal.

Conclusions
Computing reliable R rs datasets is crucial to obtaining accurate estimates, describing optical features, developing bio-optical modelling, and establishing environmental monitoring by satellite products.The main challenges for in situ measurements are to reduce the glint effects from the sky and sun light into the water-leaving radiance, or R rs .Besides, other aspects such as viewing and illumination geometry, environmental conditions (wind speed, clouds, solar incident radiation), and vertical mixing of water components produced uncertainties in R rs .
The present study evaluates four methods of glint removal and identified that the spectral method (aforementioned as M3) [5] provided the most accurate data when compared to synthetic datasets obtained from HydroLight simulations.The main advantage of computing the spectral glint correction is considered to be the spectral influence of the surface reflectance into the remote sensing signal.In addition, using a delta to remove glint residuals fits very well for almost all optical compositions evaluated in this paper.One exception needs to be highlighted: the spectral approach is too sensible when the dataset contains some outliers of Chl-a concentration.Our results also support that CDOM-rich environments are a challenging task, since the best glint removal for this case was variable ρ, which was different from all other fieldworks Our results pointed out that there is no general value of ρ factor to be adopted in different inland water conditions.On the contrary, the findings demonstrated that the most suitable methodology is the spectral ρ.That in turn means, for very unstable platforms (boats) and inland water environments (optical complex systems), non-flat surfaces considerably affect the sensor's measurements, and at some moment, the sensor is registering the sun glint effect, which needs to be spectrally corrected.All other tested methods provided less accurate R rs datasets.
Comparisons to the LaSRC product demonstrated the same trend: M3 is the best approach to provide R rs for spectral matching and the quality control of the optical information.Coastal and blue bands (centered at 442 nm and 481 nm) showed the highest errors.The error sources can be attributed to the high levels of CDOM and Chl-a concentrations that increase the absorptions in this spectral region and make the atmospheric correction a challenge.Other reasons for mismatching in lower wavelengths are related to Rayleigh scattering from atmospheric components, or adjacency effects (spectral mixture contamination), which are beyond the scope of our study.Special attention has been given to CDOM-rich environments, since the glint removal methods were sensitive to this optical composition.Regarding the other optical compositions (NAP or Chl-a dominance), the spectral glint correction when compared to the other glint-removal methods, retrieved suitable R rs and reduced the errors by almost 30%.

Figure 1 .Figure 1 .
Figure 1.Brazil with São Paulo State location (a); Map of the Tietê River Cascade System of Reservoirs (TCSR) and respective Land cover classes of Tietê basin (water body, forest, shrubland, Figure 1.Brazil with São Paulo State location (a); Map of the Tietê River Cascade System of Reservoirs (TCSR) and respective Land cover classes of Tietê basin (water body, forest, shrubland, bare soil, umid areas and urban areas by CPLA, 2010) (b); and sampling locations were detailed for each reservoir Barra Bonita (c); Nova Avanhandava (d); Bariri (e); Ibitinga (f).

Table 4 .
Descriptive statistics from all field campaigns carried out in TCSR.The notations represent: SPM-Suspended Particle Matter, PIM-Particle Inorganic Matter, POM-Particle Organic Matter, Aver-average, SD-Standard Deviation, Min-Max-minimum-maximum, CV-Coefficient of Variation.

Figure 7 .
Figure 7.Comparison between OLI derived Rrs and in situ Rrs using 34 sampling spots.

Figure 7 .
Figure 7.Comparison between OLI derived R rs and in situ R rs using 34 sampling spots.

Table 1 .
Sampling descriptions from each fieldwork in Barra Bonita, Bariri, Ibitinga and Nova Avanhandava reservoirs from TCSR, and respective water quality parameters and radiometric measurements.L t is the total radiance, L sky is the incident sky radiance (both in W•m −2 •sr −1 •nm −1 ), E s is the downwelling irradiance incident onto the water surface and E d is the downwelling irradiance (both in W•m −2 •nm −1 ).

Table 3 .
Inputs of Hydrolight of simulations using case-2 model.

Table 5 .
nRMSDs (%) for in situ R rs retrieved from M1, M2, M3 and M4 methods using HydroLight R rs as reference data.

Table 6 .
Error assessment (nRMSD) from comparisons between in situ R rs (retrieved from glint removal corrections-M1, M2, M3, M4) and simulated R rs .NRMSD (Normalized Root Mean Squared Difference, in %) and MAPE (Mean Absolute Percentage Error (in %) resulted from simulated R rs and M1, M2, M3, and M4 comparisons for the entire spectrum.