1. Introduction
Suspended particulate matter (SPM—see
Table 1 for symbols and acronyms) is a major component of the aquatic environment, composed by organic and inorganic fractions. It plays an important role in the hydrophysical functioning and biogeochemical cycles of inland waters [
1]. It essentially controls, through the absorption and backscattering of light and the turbidity and transparency of the water column, which can affect the total available energy for photosynthetic activities [
2]. Furthermore, SPM controls the transport of materials and contaminants in aquatic systems, representing an index of general uses of water resources [
3]. Thus, mapping the distribution of SPM concentration is considered critical in water resource management. High levels of SPM may alter the nutrient composition available in water and decreases the water clarity, affecting the light penetration through the water column [
4]. The gradient of available energy underwater largely determines the biogeochemical cycles and biodiversity of aquatic organisms [
5,
6,
7].
Although the SPM in inland water systems has been monitored using several different approaches [
8], remote sensing can be considered the most promising and efficient way to map the large-scale spatiotemporal dynamics of SPM [
9,
10]. Traditionally, SPM could be estimated from inherent (IOPs) or apparent (AOPs) optical properties using empirical or analytical models [
11]. Using the IOPs, SPM concentrations can be obtained by employing a spectral absorption index (SAI [
12]) or by applying empirical regressions with backscattering coefficient (b
b) [
13]. Using the AOPs, SPM can be derived from the remote sensing reflectance (R
rs) from unique or band ratios [
14,
15]. Establishing a reliable model to estimate SPM concentration using R
rs in inland waters, however still remains a challenge, because R
rs saturates at certain levels of SPM [
16]. Further, the chlorophyll-a (Chl-a), colored dissolved organic matter (CDOM), and organic and inorganic fractions of SPM respond differently to the incident energy, resulting in a wide range of magnitudes and shapes of R
rs spectra [
17,
18]. The Tietê River is a representative case that presents widely ranging IOPs, where the phytoplankton absorption coefficient (a
∅) at 443 nm ranges from 0.02 to 10.9 m
−1 [
19].
The Tietê River is the longest river in São Paulo State, running more than 1000 km before meeting the Parana River, which later becomes the Plata River reaching the Atlantic Ocean. The Tiete River is an important water resource for the communities at the local scale, providing valuable drinking water, food sources, irrigation, water for industrial use, transportation and recreation [
20]; thus, the water quality of the river is considered a critical issue in this region. A series of six hydroelectrical reservoirs along the Tietê River are constantly filtering the water, and intensive anthropogenic activities occur within the catchments (e.g., agriculture, dredging, industrial production and fishing activities). Therefore, different types of water draining the surrounding catchments are impacting directly the dynamics of the SPM in the Tiete River [
21,
22], resulting in a wide range of magnitude concentration and varying composition of the SPM. Several attempts were made to assess the spatiotemporal patterns of the SPM in the Tietê River Cascade System (TRCS); however, they were not successful in retrieving a robust result. For instance, empirical algorithms have been used to estimate the SPM over a specific [
23] or combined reservoirs [
18]. However a universal model that accounts for both a large variability (and spatially heterogeneous) in SPM concentrations and its varying biogeochemical composition still does not exist by considering K
d, retrieved from an semianalytical scheme, as the main predictor of SPM in wide ranges.
In this paper we used diffuse attenuation coefficient (K
d) as a key parameter for the SPM retrieval. K
d depends on several factors, such as depth (z), incident light (E
d), IOPs (a and b
b) and optically significant constituents (OSCs—Chl-a, CDOM and SPM). In the inland water systems such as the TRCS, SPM concentrations are the main OSC that control the K
d values [
24]. As an AOP, K
d is largely determined by the IOPs (a and b
b) and secondarily on the light field geometry [
25]. Therefore, K
d is considered more stable than the R
rs in estimating SPM, because it is directly derived by summing the IOPs [
25], while R
rs is a function of the ratio of IOP, due to the analytical configuration of absorption and backscattering. In addition, K
d is also relatively easier to validate compared to SPM models based on b
b [
26,
27]. Thus, K
d can be used efficiently to map SPM dynamics over inland waters where the SPM concentration varies widely, which has only been used so far in some coastal waters [
27].
We aim to develop a semi-analytical scheme to estimate SPM in the TRCS that is sensitive enough to capture the widely varying effects of OSCs. For this, K
d is estimated based on the methods by Lee et al. [
28], derived from the IOPs (a and b
b) estimated using a re-parameterized version of a quasi-analytical algorithm (QAA [
29]). The QAA was tested for the TRCS [
20,
30] and the results suggested that further re-parameterizations of the QAA are required for more accurate results. The K
d model from Lee et al. (2013 [
28], 2015 [
31]) was also evaluated for other inland waters [
22,
30], which provided the lowest error when compared with other published algorithms to retrieved K
d. However, the suitability of Lee’s model was never tested for the entire TRCS where the optical composition of surface water varies over reservoirs. Therefore, the specific objectives of this paper are to (i) re-parameterize the QAA considering the optical composition and the available spectral bands onboard Landsat-8, (ii) compare our results with other versions of QAA tuned for inland waters, (iii) evaluate the suitability OLI bands to estimate Kd for the first time, and (iv) evaluate the influence of different water compositions on SPM estimates.
4. Discussion
QAA
TRCS provides the most accurate results, mainly for the a
t estimates. The improvement of its performance is related to the four main modifications we made in deriving the QAA
TRCS. The first modification was adapting the method from Wang et al. [
55] to compute α(λ) and β(λ) coefficients instead of using 0.52 and 1.7 values to retrieve r
rs. The second change was shifting λ
0 towards longer wavelengths. Integrating the OLI/Landsat-8 bands and our absorption measurements (400–800 nm), we selected two wavelengths for this—561 nm and 665 nm. Nevertheless, when we assessed the total absorption against water absorption contribution, we found that a
t-w(655) complies with Lee et al.’s (2002) [
29] requirements for choosing the reference wavelength.
Another change was identifying the α coefficient (see
Table S1) that relevantly estimates a and b
b. For this, we tested α=2 and α = 5 [
29,
64], which presented huge discrepancies in their errors, reaching almost 20% in some wavelengths (results of our tests not shown in this paper). In regard to the performances, we used α = 2 for reservoirs with higher CDOM contributions (BAR2 and IBI2) and α = 5 for the others. The final change was modifying the bands to compute the η parameter, since the two reservoirs (BAR and IBI) presented a high contribution of CDOM into the absorption conditions (
Figure S2). We used a 561/655 band ratio to account for the CDOM effects [
56]. Among all tested bands, the 561/655 band ratio also provided the best correlation with SPM concentrations for the TRCS’s dataset (r = 0.65). As a result, the QAA
TRCS significantly improved its performance in estimating a
t and b
b by 14% and 30% respectively from QAA
OMW and QAA
V5. The main reason for the relatively poor performance of a
t estimates in BB2 (
Table S3) is related to the high levels of Chl-a. Comparing the results with eutrophic aquatic systems by Watanabe et al. (2016) [
20] and Mishra et al. (2014) [
46], the poorest performance is likely to be caused by the λ
0, since OLI did not have the 709 nm (or near) band to be used as a reference wavelength. It is important to highlight that QAA
TRCS presented the lowest average error of at estimates (
Table 7) with a MAPE = 30.7%, when compared to the lowest performance that retrieved a MAPE = 43.9%, which probably can be caused by the QAAOMW that was developed for the inorganic environment and failed when tested in a more eutrophic environment such as BBHR (as demonstrated in 18), BAR or IBI. QAA
CDOM and QAA
TRCS retrieved nRMSE of 17.7% and 16.8%, respectively for average a
t estimates. Considering that QAA
CDOM was developed for environments with high levels of CDOM such as Itumbiara, it was expected that QAA
CDOM would develop a good performance in BB, BAR and IBI, which presented high a
cdom coefficients.
Regarding the backscattering, the estimated values agreed with in situ b
b, except in BB1. The backscattering measurements in BB1 and NAV1 were conducted using HydroScat-6P (HOBI-labs Inc.2008), which was originally designed for ocean waters [
65]. Therefore, when the sensor is used in turbid water with high scattering and absorption properties, the measurements are susceptible to the signal losses of path length and saturation [
66]. BB1 is optically more complex than NAV1, where the measurements can be relatively unstable due to limitations of the equipment [
66,
67]. An additional source of errors in BB1 can be related to the post-processing of the HydroScat-6P data, which includes the corrections of power losses due to the sensor’s path length, also known as the Sigma correction. Even when processed according to the manufacturer’s instructions, it may still contain some unexpected variations within the blue-green spectral range when the surveying environment does not completely meet the desirable usage conditions defined for the HydroScat [
67].
It is important to highlight that QAA
OMW also presented adequate performances in estimating b
b during field campaigns, mainly for non-turbid waters (turbidity <6 NTU). Overall, differences in b
b estimates for QAA versions were less expressive, with QAA
TRCS retrieving 39.5% of MAPE and QAA
CDOM retrieving 39.5%, which can be considered as statistically equal results. When we considered the errors retrieved in the TRCS (
Table S3), we observed that QAA
TRCS retrieved the lowest error (nRMSE = 18.70%) when compared to QAA
OMW (nRMSE = 19.70%); however, when we evaluated each fieldwork, we verified the lowest errors in QAAOMW, with the exception of BB and BAR, which are more eutrophic environments and failed the QAA
OMW estimates for b
b. Differences in MAPE of the b
b estimates between QAA
OMW and QAA
TRCS in non-turbid waters ranged from 5.9% (IBI2) to 16.6% (NAV2). A possible source of this difference could be the band ratio used to compute η. Rodrigues et al. [
18] used 655/754 nm, which are the bands not available in OLI. We tested all band combinations to provide η values close to the ones retrieved from QAA
OMW, and the best result was achieved with 561/655 nm. Despite these differences, the magnitude of estimated b
b via QAA
TRCS did not affect the final K
d estimates (unlike the case for the a
t estimates), since η values are mostly influential over shorter wavelengths [
49].
Overall, QAACDOM and QAATRCS also presented similar performances for at and bb estimates. Such results confirmed that the inclusion of 561 nm was important in improving the QAA performance. Additionally, computing rrs using spectral coefficients instead of fixed values and using an interchangeable value of α also contributed to improving the IOPs estimates.
Outputs of QAA
TRCS were used in the K
d equation (Equation (1)). The model published by Lee et al. (2013) derives K
d from IOPs and is considered a semi-analytical model which provides reliable estimates for inland waters according to Gomes et al. (2018) [
22]. When comparing K
d_QAA and K
d_r, resulting errors were less than 25% for the entire TRCS’s dataset, with the minimum at 561 nm (19.2%). The poorest performances were found in BB2, BAR1 and BAR2, which could be affected by the complexity of water types—the presence of relatively high concentrations of OSCs variably impacts the light attenuation and consequently produces higher values of K
d that were not considered by Lee et al. (2013). This was clear when we observed the higher errors in shorter wavelengths for BB2, BAR1 and BAR2, which indicates an additional effect of CDOM absorption into K
d.
The highest modeled K
d value reached 8.5 m
−1, whilst the highest measured value was 11.9 m
−1, which is almost 30% less than the maximum reference value. The same error, about 30%, was found for the minimum value of K
d. Overall, the average error of about 21% is a satisfactory result (
Table 8). Values of K
d decrease from upstream to downstream (as observed in
Figure 4), which are directly related to the SPM concentrations, which, with isolated peaks of SPM that imply to K
d peaks. Higher values of K
d(443) and K
d(482) also confirmed that absorption from CDOM and phytoplankton in shorter wavelengths highly contributes to the K
d values, attenuating the light field inside water. Another observable trend in
Figure 5 is that the accuracy of K
d estimates depends on spectral zones and that K
d in longer wavelengths are more precise than in shorter wavelengths (also presented in
Table 8).
Once 655 nm was determined as the most suitable wavelength, we used this band of Landsat 8/OLI to estimate a
t, b
b, K
d and then finally to derive SPM concentrations (
Figure 6). R
rs(655) values were higher in the BB2 image, whilst the values were the lowest in NAV. Given that all Landsat 8 images used in this study were atmospherically corrected (LASRc product [
60]), we conclude that the variances observed in each image arise from the widely varying OSCs.
In the BB, although it is an accumulation reservoir, distribution of SPM is rather homogeneous, and a typical decreasing pattern of SPM concentration downstream is not observed. We consider that at this time the dam was in operation releasing the water because October is in the middle of wet season. The SPM map of NAV collected in May, which is close to the peak of the hydrograph, also did not show any longitudinal trend. In contrast, a clear longitudinal pattern of decreasing SPM toward the dam is detected in the BAR. Although this is a run-of-river dam, since the image was acquired during the dry season (August), the reservoir seems to be storing water and thus temporally trapping sediment. The same logic is applied to the IBI, which shows a decreasing gradient towards the dam.
The tributaries of each reservoir presented higher SPM concentrations than the main channel of reservoirs, which indicates the SPM contribution from the tributaries of the Tietê River. It is noteworthy that the OLI images we used to estimate SPM were acquired on the same day (or near) with our field campaigns, enabling direct validation of SPM estimates. Now that we have a field-validated semi-analytical model, it is important to reconstruct a time series SPM map for continuous monitoring of SPM dynamics and to build standards of SPM concentrations in the entire cascade using Kd as a predicting parameter.
5. Conclusions
SPM directly impact the biological aquatic process due to several effects, such as adsorbing contaminants, increasing the temperature by absorbing heat, and affecting the penetrability of light within the water column. Regarding the light attenuation caused by SPM within the water, we observed throughout our experiments that Kd, which represented the light attenuation, was a suitable single explanatory variable to estimate SPM concentrations for inland waters. Due to the optically complex characteristics of such systems, taking into account the OSCs variation is an important issue to develop an accurate optical model, especially when we use analytical models to derive IOPs, which are described as a function of the concentration of OSCs. We used the most applicable analytical model, the QAA, and adjusted it from existing schemes to be suitable for OLI sensors, and also to be applicable to inland waters with widely varying OSC concentrations.
Our changes in the QAA consist of using an interchangeable parameter for the CDOM environment, as well as adapting the spectral methodology to compute r
rs instead of using fixed values as originally proposed in Lee et al. (2002) [
29]. The main reason to use spectral coefficients for computing r
rs is its spectral dependence. Another relevant modification was including a specific band to incorporate and mitigate the CDOM absorption effects, mainly for brownish waters such as the BAR and IBI reservoirs in the TRCS. The adoption of a 561 nm band retrieved the lowest errors in all tested versions of the QAA.
Our re-parameterized model, the QAATRCS, improved the IOP estimates, yielding better accuracies for at, bb and consequently Kd estimates. Kd was responsible for explaining over 74% of SPM variation among the widely varying SPM concentrations in the TRCS. Then, the predictive SPM values were retrieved by using Kd(655) and a power fitting, capable of providing estimates with errors less than 30%.
Changes that were made in bands for some parameters and spectral optimization also implied some constraints. One of them was that OLI bands were not sensitive enough to estimate in more eutrophic environments; however, the relative error did not affect the SPM concentration estimates by Kd. The 655 nm band was the most suitable band of the OLI sensor to derive Kd, since the at coefficient was dominated by water, at least for 70% of the TRCS’s dataset. Using other sensors might improve the performance in SPM estimates when SPM > 30 mg.L−1—the highest estimated SPM value using QAATRCS, which is indicative for further research using other satellites, such as Sentinel-2A, which presented spectral bands near 700 nm.
Throughout the cascade, Kd showed a decreasing gradient from upstream to downstream (along with the SPM variation). In addition, the highest values of Kd(443) confirmed that CDOM and phytoplankton absorption are markedly representative in Kd values. Spatial distribution of SPM is homogeneous in the downstream reservoir, while in some intermediate reservoirs in the TRCS, a clear gradient towards to the dams was presented.
In conclusion, QAATRCS was capable of deriving Kd in inorganic, organic and CDOM dominant aquatic systems and providing reliable SPM estimates for the entire TRCS. Further investigations are needed to assess the suitability of using a sensor that has a 700-nm band and adequate spatial resolution to capture moderate to high SPM concentrations. Once the model was validated with in situ measurements, time series SPM could be reconstructed to identify the environmental standards of SPM. Future investigations can apply the QAATRCS as an analytical model to compute Kd and consequently estimate SPM concentrations over the entire cascade, aiming to identify the SPM standards and eventual drivers to extreme SPM values.