Retrieving Total and Inorganic Suspended Sediments in Amazon Floodplain Lakes: A Multisensor Approach

Abstract

Remote sensing imagery are fundamental to increasing the knowledge about sediment dynamics in the middle-lower Amazon floodplains. Moreover, they can help to understand both how climate change and how land use and land cover changes impact the sediment exchange between the Amazon River and floodplain lakes in this important and complex ecosystem. This study investigates the suitability of Landsat-8 and Sentinel-2 spectral characteristics in retrieving total (TSS) and inorganic (TSI) suspended sediments on a set of Amazon floodplain lakes in the middle-lower Amazon basin using in situ Remote Sensing Reflectance (R_rs) measurements to simulate Landsat 8/OLI (Operational Land Imager) and Sentinel 2/MSI (Multispectral Instrument) bands and to calibrate/validate several TSS and TSI empirical algorithms. The calibration was based on the Monte Carlo Simulation carried out for the following datasets: (1) All-Dataset, consisting of all the data acquired during four field campaigns at five lakes spread over the lower Amazon floodplain (n = 94); (2) Campaign-Dataset including samples acquired in a specific hydrograph phase (season) in all lakes. As sample size varied from one season to the other, n varied from 18 to 31; (3) Lake-Dataset including samples acquired in all seasons at a given lake with n also varying from 17 to 67 for each lake. The calibrated models were, then, applied to OLI and MSI scenes acquired in August 2017. The performance of three atmospheric correction algorithms was also assessed for both OLI (6S, ACOLITE, and L8SR) and MSI (6S, ACOLITE, and Sen2Cor) images. The impact of glint correction on atmosphere-corrected image performance was assessed against in situ glint-corrected R_rs measurements. After glint correction, the L8SR and 6S atmospheric correction performed better with the OLI and MSI sensors, respectively (Mean Absolute Percentage Error (MAPE) = 16.68% and 14.38%) considering the entire set of bands. However, for a given single band, different methods have different performances. The validated TSI and TSS satellite estimates showed that both in situ TSI and TSS algorithms provided reliable estimates, having the best results for the green OLI band (561 nm) and MSI red-edge band (705 nm) (MAPE < 21%). Moreover, the findings indicate that the OLI and MSI models provided similar errors, which support the use of both sensors as a virtual constellation for the TSS and TSI estimate over an Amazon floodplain. These results demonstrate the applicability of the calibration/validation techniques developed for the empirical modeling of suspended sediments in lower Amazon floodplain lakes using medium-resolution sensors.

1. Introduction

Amazon floodplains, mainly those associated to white water rivers [1] such as the Solimões/Amazonas and Madeira Rivers have an important role in the regional-scale biogeochemical processes because they are responsible for almost all of the suspended sediments transported from the Andes to the Atlantic Ocean [2,3]. Particularly, the lower Amazon floodplain is paramount for the understanding of sediment dynamics, as it receives all the sediments transported from the Andes which were not stored on high and medium Amazon floodplains up to the confluence of the Madeira and Solimões Rivers, [2,4,5]. This dynamic is mainly dependent on hydrological and climatic factors such as the magnitude and duration of floods/droughts [6,7,8,9] and on anthropogenic activities, like land use and land cover within the floodplain [10,11,12].

Hydrological and geomorphological ---factors are mainly responsible for the sediment supply to the floodplains, which occurs primarily in high-water seasons through overbank flow and levee breaches [7,13]. Part of these sediment budgets, estimated in approximately 80% [6], is stored in the floodplain through sedimentation and the building of an intricate scroll-bar topography of the floodplain [6,7,13,14,15]. However, the floodplain acts as a temporal storage [6] because, during the receding phase, there is an outflow which transports the remaining suspended sediments as well as those resuspended by waves and currents into the Amazon River [7,16]. Besides natural factors, the floodplain is also affected by anthropogenic factors related to the conversion of a flooded forest into agriculture [10], cattle ranching [12,17], and hydroelectric reservoirs in the Andes [18]. Those factors also may impact the sediment budget in this region, however, with less intensity in relation to the natural factors [2].

Under the pressure of anthropogenic factors [19,20] and climate changes [21,22], the monitoring of the sediment concentration patterns in Amazon floodplains becomes a key indicator of the basin resilience to those impacts. As Total Suspended Solids (TSS) fluxes are fundamental to biogeochemical processes and to the biodiversity of floodplains, their spatial and temporal dynamics have been the object of several studies [6,7,8,13,16,23,24,25,26,27,28,29,30,31]. Moreover, the quantification of sediment concentration contributes to more accurate carbon budget estimates, as suspended organic solids represent an important fraction of the suspended solids (TSS) [12,32] present in the floodplains, which are one of main sources of CH₄ and CO₂ outgassing [33] to the atmosphere [25,34,35,36].

Much of this knowledge, however, is based on a limited number of in situ samples and, therefore, is lacking a spatial portrayal of TSS and TSI (Total Suspended Inorganic Solids) distribution. As they are based on local studies, in situ methods become time- and cost-demanding which restricts their frequent application to regional and large-scale studies [37,38,39]. With this in mind, remote sensing (RS) techniques represent an efficient tool to enhance the spatial representation of TSS estimates through mathematical modelling, as it provides a synoptic view of water bodies with a higher temporal frequency over large areas, favoring the characterization of time changes in the amount of sediment exchange between the Amazon River and its floodplains [39,40,41].

The monitoring of suspended sediments in the Amazon Basin rivers has been carried out with MODIS (Moderate-Resolution Imaging Spectroradiometer) imagery, particularly in the Curuai Floodplain [42,43,44]. However, low spatial resolution sensors prevent the assessment of small channels (<0.25 km²) and river sediments plumes [44,45], which are often overlooked. The new generation of medium-resolution satellite imagery is a suitable alternative to overcome this limitation, increasing the opportunity of quantifying TSS and TSI in a broad range of lakes and rivers comprising the Amazonian aquatic systems. Many studies have been focused on using one sensor (Landsat-8 or Sentinel-2) at a time [19,45,46,47,48,49,50,51,52]. However, recent efforts to combine the data from both satellites have been currently performed by different research groups such as the Harmonized Landsat Sentinel product (HLS), which provides seamless surface reflectance products for both sensors [46] and the assessment of consistency performance between Landsat-8/OLI (Operational Land Imager) and Sentine-2/MSI (Multispectral Instrument)-suspended sediment products provided by Pahlevan et al. [47] which compared images acquired during concurrent satellite overpasses. The combined data from OLI and MSI increases the frequency of image acquisition, creating a robust database for monitoring both the riverine sediment transport and its accumulation in lakes of the Amazonian Floodplains as it increases the probability of cloud-free images [48] since the low frequency of data acquisition is the main constraint for the remote sensing of the Amazon region because of the constant cloud cover in most seasons [49,50].

The remote sensing of TSS is based on the interaction between the electromagnetic radiation and Optically Actives Constituents (OACs) of the aquatic system. This interaction is affected by the physical proprieties of particles [51] such as the index of refraction [52] and Particle Size Distribution (PSD) [53]. Inorganic particles have a higher index of refraction compared to organic particles and, consequently, higher backscattering [53]. The PSD also influences the scattering of light according to the dominant particle within the water body, with an increase in particle backscattering as PSD becomes dominated by smaller particles [54].

The Amazon floodplain waters are characterized by a high spatial and temporal variability in particle proprieties, which are dependent on both TSS origin and hydrograph phase [8,28]. TSS transported from the Amazon River to the floodplains can be both inorganic (TSI) or organic (TSO) depending on the proportion of local to upstream sources. TSI is originated mainly from the erosion processes occurring at the Andes headwaters [2,51,55] modulated by TSS fluxes from the Solimões River (draining the Peruvian Andes) [56] and from the Madeira River (that also drains in part in the Bolivian Andes and in part in the cratonic Brazilian Shields) [2,57], both rivers responding to about 98% of the Amazon River sediments discharge [58]. At the floodplains, TSS originates also from influxes of the Amazon River as well from channel erosion and floodplain lake water resuspension by wind during the low-water period [7,8,28]. TSO is originated from surface runoff that brings organic matter from both Terra-Firme (dry land) and flooded forests [25,30,59] as long as it is from algal blooms [32,59,60].

It is well-documented that Remote Sensing Reflectance (R_rs) and TSS have a high correlation in green, red, and near-infrared (NIR) spectral bands [61], although its magnitude depends on both the proportion and types of OACs in the water column [38,62]. At lower concentrations (up to 50 mgL⁻¹) [38] algorithms based on green or red spectral bands are the most suitable for modeling TSS concentrations [63]. At a higher concentration, above 100 mgL⁻¹, approximately, R_rs in green and red bands saturate and the concentration tends to be underestimated. Therefore, NIR bands become the obvious choice [38,51,62,64]. However, due to the optical complexity and OAC variability at Amazon floodplains [8,28], the development of an overarching algorithm applicable to the different lakes and in different seasons, is still a challenge.

There are several factors which affect the accurate estimate of surface reflectance based on satellite imagery, consequently limiting the accuracy of TSS retrieval. These factors are mainly associated with an electromagnetic interaction with the atmosphere and to the signal reflected at the water surface into the sensor (glint) [53]. The signal registered by satellite images are affected by molecular and aerosol scattering and absorption by optically active gases [65,66]. Besides atmosphere influence, the R_rs signal is also affected by the sky and sun reflected at the water surface (sun glint and sky glint), demanding their correction [19,67,68]. Therefore, the quantification of atmosphere and glint signal is a crucial step for an accurate modeling of TSS and TSI using orbital imagery [63,69]. Several algorithms for atmosphere [65,70,71,72,73] and glint [67,68,74] correction are available in the literature. However, their accuracy is highly variable, according to atmosphere characteristics [71,75,76], water type [66,77,78], and the geometry of data acquisition (e.g., field of view and view angle) [67,79,80].

Therefore, the aim of this paper was to calibrate and to validate for the first time several empirical algorithms based on both in situ OLI and MSI simulated R_rs, which were then tested to OLI and MSI imagery on the lower Amazon floodplain lakes. As atmospheric effects limit TSS and TSI accuracy, this paper also reports the assessment of three different atmosphere correction methods for each sensor as follows: 6S [65], ACOLITE [81], and Landsat 8 Surface Reflectance Product (L8SR) [75] for OLI images and 6S, ACOLITE, and Sen2Cor [73] for MSI images. Both sets of atmospherically corrected images were examined with and without glint correction [68]. The calibrated empirical models were then applied to OLI and MSI scenes in order to spatialize the TSS and TSI distribution.

2. Materials and Methods

To address the objective, the research effort is described as follows: (i) The calibration and validation of empirical TSS and TSI algorithms through OLI- and MSI-simulated R_rs using a Monte Carlo Simulation; (ii) the validation of atmosphere and glint correction for OLI (6S, L8SR, and ACOLITE) and MSI (6S, ACOLITE, and Sen2Cor) using in situ R_rs; (iii) the validation of TSS and TSI models obtained through Monte Carlo simulation in corrected OLI and MSI imagery; (iv) the spatialization of TSS and TSI in the study area.

2.1. Study Area

The study area is located between the Trombeta River mouth to the west (close to Paru Lake) and Monte Alegre city to the east (Figure 1), a region representative of the lower Amazon floodplain which is characterized by large and shallow lakes [8]. The region is exposed to an extreme fragmentation of the floodplain forest, which the cover of has been reduced by 56% in the lasts 40 years (1975–2008) [10], being firstly occupied by jute farming and later by cattle ranching [17,20]. The large lakes in this region create local floodplain systems which the names of are associated with those lakes less affected by the water-level change along the hydrological year.

Among these floodplains systems, the Lago Grande de Curuai (LGC) floodplain stands out. Localized 900 km upstream from the Amazon River mouth, this complex continental aquatic system contains more than 30 interconnected lakes linked to the Amazon River by several permanent and temporary channels. The surface of the LGC open water varies from about 600 km² at the low-water level season (November) to 3500 km² at the high-water season (June) [6]. The LGC hydrological variability is mainly dependent on Amazonas River flow (70–90% of inflow fluxes), which results in about three months of residence time [16] favoring the occurrence of biogeochemical processes [6,8]. According to Barbosa et al. [8], the hydrological dynamic at LGC is characterized by four flow conditions: State 1 (Rising), between January and March; State 2 (High Water), between April and June; State 3 (Receding), between August and October; and State 4 (Low Water), between October and December (See Barbosa [26] and Barbosa et al. [8] for further details). Regarding sediment dynamics at LGC, Barbosa et al. [8] evaluated TSS and TSI concentrations between 2003 and 2004 in the four states and observed higher TSS and TSI values during the Low and Rising states, while lower concentrations were observed during the High and Receding states.

Besides the LGC floodplain, one can note the Paru Lake floodplains system, mainly influenced by the Amazon River flood pulse but less impacted by anthropogenic actions [82]. The floodplains systems Pacoval, Aramanai, and Monte Alegre, however, were influenced by both Amazon floodplains and the Tapajós River flood pulse in some hydrograph phases [83].

2.2. Available Dataset

2.2.1. Limnological Data

The limnological dataset used in this study (TSS, TSI, TSO, and Chlorophyll-a (Chl-a)) were collected during four field campaigns carried out between 2015 and 2017 (Figure 2), sponsored by the project MAS-BNDES 1022114003005. The TSS, TSI, and TSO quantifications followed Wetzel and Likens [84], using a Whatman GF/C (1.2 µm) glass fiber filter for filtration. Chl-a concentration determination followed Nush [85], using a Whatman GF/F (0.7 µm) glass fiber filter for filtration. Both the TSS and Chl-a concentrations were measured in duplicate, giving a concentration value as a simple mean of the duplicates.

2.2.2. Radiometric Data

Three intercalibrated TriOS-RAMES spectroradiometers operating in the 350–950 nm spectral range obtained the radiometric measurements. Each instrument measured, simultaneously, the downwelling irradiance above the water surface (E_d (λ)), Total Water-Leaving Radiance (L_T (λ)), and Sky Radiance (L_sky (λ)). All sensors were set 5 m from the water surface to avoid shadows and boat reflections. All samples were acquired according to Mobley [86] protocols: L_T was acquired with an approximately 45° zenith angle ($\theta$) in relation to the surface and an approximately 135° azimuthal angle ($\varnothing$) in relation to the sun; L_sky were acquired through a 90° (θ′) in relation to L_T and at the same L_T azimuthal angle; and E_s was positioned 180° in relation to the water surface. The R_rs computation followed Equation (1), where the $\rho \left(\theta ,\varnothing \right)$ values were from Mobley [87] and the wind was measured with an anemometer during field campaigns.

$$R_{rs}\left(\theta ,\varnothing ,\lambda \right)=\frac{L_T\left(\theta ,\varnothing ,\lambda \right)}{E_d\left(\lambda \right)}-\frac{\rho \left(\theta ,\varnothing \right)L_{SKY}\left({\theta }^{\prime },\varnothing ,\lambda \right)}{E_d\left(\lambda \right)}·\left(s{r}^{-1}\right)$$

(1)

Around 150 measurements were taken during 30 minutes (approx. 1 spectrum at each 10 s) at each sampling station. The entire set of the spectrum was inspected, and conspicuous outliers were removed. Then, the spectrum with the smaller value of the difference between the median of all the spectra at each wavelength and the remaining spectra measurements (Equation (2)) was used to represent each station.

$$Dif{R}_{rs,i}={{\displaystyle \sum }}_{\lambda =400}^{900}\left|{Rrs}_{i,\lambda }-{Rrs}_{\begin{array}{c}Median,\lambda \end{array}}\right|$$

(2)

where $Dif{R}_{rs,i}$ is the sum of the difference between median R_rs and each in situ-derived R_rs at each wavelength for the i sample, R_rs,i,λ is the R_rs at sample i and wavelength λ, and $Rr{s}_{Median,\lambda }$ is the median R_rs for each station and wavelength (i.e., the median for all spectra measured at each station). The spectrum having the smaller $Dif{R}_{rs,i}$ was then selected. After that, the R_rs values were interpolated for 1 nm interval (originally approx. 3.3 nm) and used to simulate R_rs in spectral bands of OLI [88] and MSI-A [89] sensors based on the Spectral Response Function of each sensor available at https://landsat.gsfc.nasa.gov/wp-content/uploads/2013/06/Ball_BA_RSR.v1.1-1.xlsx for OLI and at https://earth.esa.int/web/sentinel/user-guides/sentinel-2-msi/document-library/-/asset_publisher/Wk0TKajiISaR/content/sentinel-2a-spectral-responses for MSI (Equation (3)).

$$R_{rs,sim\_sensor-band}=\frac{{{\displaystyle \int }}_{{\lambda }_n}^{{\lambda }_m}SRF\left(\lambda \right)\ast R_{r{s}_{m\left(\lambda \right)}}d\lambda }{{{\displaystyle \int }}_{{\lambda }_n}^{{\lambda }_m}SRF\left(\lambda \right)d\lambda }$$

(3)

where $R_{rs,sim\_sensor-band}$ is the simulated R_rs for an ith spectral band of both the OLI and MSI sensors, SRF $\left(\lambda \right)$ is the spectral response function of each sensor at a specific band, and $R_{r{s}_{m\left(\lambda \right)}}$ is the measured R_rs for the interval n -> m that accounts for the OLI and MSI bands. The simulated R_rs will follow the below nomenclature: R_{rs,sim_OBi} (e.g., R_{rs,sim_OB2} for OLI band 2) and R_{rs,sim_MBi} (e.g., R_{rs,sim_MB2} for MSI band 2).

2.2.3. Satellite Data

Satellite imagery from both sensors (OLI/Landsat-8 and MSI/Sentinel-2A) were used. Despite their distinct characteristics (spectral, spatial, and radiometric resolutions) as well as swath width and the Signal-To-Noise Ratio (SNR) (

Table 1. The spectral bands; spectral, spatial, radiometric, and temporal resolutions; swath; and Signal-to-Noise Ratio (SNR) for each OLI and MSI bands. * When both Sentinel-2A and Sentinel-2B are considered.

Sensor/Bands	Spectral Bands Interval (nm)	Spatial Resolution (m)	Radiometric Resolution (bits)	Revisit Time (days)	Swath (km)	SNR at Lref
OLI/B2	452–512	30	16	16	185	321
OLI/B3	533–590	30	16	16	185	223
OLI/B4	636–673	30	16	16	185	113
OLI/B5	851–879	30	16	16	185	45
MSI/B2	447.6–545.6	10	12	5*	290	439
MSI/B3	537.5–582.5	10	12	5*	290	102
MSI/B4	645.5–683.5	10	12	5*	290	79
MSI/B5	694.4–713.4	20	12	5*	290	45
MSI/B6	731.2–749.2	20	12	5*	290	45
MSI/B7	768.5–796.5	20	12	5*	290	34
MSI/B8	762.6–907.6	10	12	5*	290	26
MSI/B8-a	848.3–881.3	20	12	5*	290	20

Source: Barsi et al. [88]; Pahlevan et al. [90].

), several studies have already indicated that their accuracy, regarding TSS concentration estimates, is equivalent [63,70,81,90].

Three atmosphere correction methods were assessed for each sensor. For the Landsat-8/OLI imagery, the following methods were applied: (a) the NASA product Landsat 8 Surface Reflectance (L8SR) in surface reflectance that is based on the 6SV model [75] and are reported as having good performance for inland water studies [64,69,77]; (b) the 6S (Second Simulation Of The Satellite Signal In The Solar Spectrum) [65] radiative transfer code with atmospheric parameters based on Martins et al. [66]. A modified version of Py6S [91] developed at The Instrumentation Laboratory for Aquatic Systems (LabISA: http://www.dpi.inpe.br/labisa/) was used to run 6S (unpublished); and c) ACOLITE [71] (version 20180611.0), an image-based atmosphere correction method that allows for simple and fast processing for coastal and inland waters applications [63,66]. Shortwave Infrared (SWIR) bands were selected for the aerosol atmosphere scattering calculation due to the non-negligible signal at near-infrared (NIR) related to a high sediment concentration. For the Sentinel 2A-MSI imagery, the following methods were applied: both 6S and ACOLITE already tested in OLI imagery and the Sen2Cor (version 2.5.5), another image-based method [73]. After atmospheric correction, the surface reflectance was divided by π to obtain an atmosphere-corrected remote sensing reflectance for both sensors (R_{rs_sat_ac}).

After atmospheric correction and R_{rs_sat_ac} conversion, the Wang and Shi [68] glint correction was applied to both sensors scenes. This glint correction assumes that the remaining signal of SWIR R_{rs_sat_ac}. comes from a water surface specular reflection, as the SWIR signal is negligible even at higher TSS and TSI concentrations [19,68]. Therefore, Equation (4) was used to obtain the $R_{rs,sat\_deglint}\left(VNIR\right)$ in visible and NIR OLI and MSI bands.

$$R_{rs,sat\_deglint}\left(VNIR\right)=R_{rs,sa{t}_{ac}}\left(VNIR\right)-R_{rs,sa{t}_{ac}}\left(SWIR\right)$$

(4)

The OLI and MSI images used in this study were obtained from USGS (United States Geological Service) and correspond to path-row 228/61 (OLI), Aug/10/17 and tile T21MXT (MSI), Aug/08/17. The August 2017 field campaigns were carried out between days 8–12 at Lago Grande de Curuai (LGC). Only this campaign could be used for validation purposes because of the low cloud-cover images acquired concurrently to the field campaign.

2.3. Empirical Models

Three types of empirical algorithms were selected for modeling TSS and TSI: (1) single-band linear algorithms; (2) single-band log-transformed algorithms; and (3) band-ratio algorithms. The simulated R_rs,sim models were assessed for each visible and near-infrared (VNIR) spectral bands of both the OLI and MSI sensors with and without log-transformation. Furthermore, the band-ratio algorithms based on a nonlinear regression analysis (Y = aX^b, where X is the band ratio) [51,92] were also evaluated using ratios between NIR bands and visible (VIS) bands [92,93]. The semi-analytical approach proposed by Nechad et al. [62] was also evaluated. However, this method required a local recalibration for our study area and did not yield accurate results. For these reasons, the Nechad algorithm was not considered further in this study. All statistical analyses were performed using R version 3.5.0 [94].

2.4. Model Validation and Application to Satellite Images

The calibration/validation of TSS and TSI empirical algorithms was based on Monte Carlo (MC) simulations with 10,000 repetitions (Figure 3), similar to the approach adopted by Augusto-Silva et al. [95]. At each repetition, 70% of the dataset was randomly selected for training and 30% was randomly selected for the validation of the algorithm. For each Monte Carlo trial, the Mean Absolute Percentage Error (MAPE: Equation (5)), the Coefficient of Determination (R²), and the Root Mean Square Error (RMSE: Equation (6)) were calculated. For each model, the statistical metrics were obtained through the mode values of the 10,000 statistical metrics obtained from MC simulation.

The following datasets were used in the calibration/validation process: (1) a full dataset (n = 94), including all field campaign and dates, which, from now on, will be named All Data; (2) a data split by field campaigns (See Figure 2 for sample sizes), named from now on as Campaign Data; and (3) the data separated by lakes (Curuai Lake (n = 66) and Monte Alegre Lake (n = 17)) named from now on as Lake Data.

$$MAPE=100\ast \overline{{\displaystyle {\displaystyle \sum }_{i=1}^n}\frac{(y_i-x_{i)}}{xi}}$$

(5)

$$RMSE=\sqrt{\overline{{\left({\displaystyle {\displaystyle \sum }_{i=1}^n}(x_i-y_i\right)}^2}}$$

(6)

where $x_i$ is the in i, $y_i$ is the predicted value for station i, and x_max and x_min are respectively the maximum and minimum TSS or TSI values for each dataset.

MAPE values were used as the primary criterion in the model selection because it provides the absolute percentage of errors, allowing a comparison among different models based on different datasets and concentration ranges [96]. The model coefficients (e.g., slope and intercept) were selected by minimizing the difference between their values and the computed MAPE mode. Figure 3 presents the flowchart of the calibration/validation model process.

The performance of the different atmosphere correction methods with and without glint correction was assessed by comparing R_{rs,sat_AC} and R_{rs,sat_AC_deglint} with R_rs,sim for both sensors. The sampling points subjected to cloud contamination were removed, and the time span between the date of satellite overpass and the in situ data acquisition was kept ≤4 days. Although this time span between satellite overpass and in situ data acquisition is a source of uncertainty in the algorithm calibration, this study assumes that this time span is not a critical constraint to the model validation [97]. The R_{rs,sat_AC} and R_{rs,sat_AC_deglint} from OLI (n = 16) and MSI (n = 14) were the mean value of a 3 × 3 pixels window to reduce the SNR effects [98]. The results were assessed for each spectral band and for all VNIR bands, similarly to Wei et al. [77]. The R² and MAPE statistics were used to evaluate the atmosphere and glint correction for both sensors. The calibrated models for the All Dataset were applied to both OLI- and MSI-corrected imagery, and the R², MAPE, and RMSE statistics were calculated based on in situ TSS and TSI measurements (n = 16 for OLI and n = 14 for MSI) for the 2017 campaign. Therefore, the best models based on the validation statistics for the All Dataset were selected, and the correspondent models were compared with those provided by the Campaign and Lake data sets. Thus, the best model for each sensor/type was selected and applied to the images.

A comparison between the cloud-free images acquired at concurrent overpasses of both sensors were corrected for atmospheric effects using 6S and also by glint effects. Then, TSS-calibrated algorithms were applied to the corrected images to assess the congruence of the TSS estimates [99]. A search for matchups at Curuai Lake from 2015 to 2019 resulted in two scenes for the low-water season (Nov/11/16 and Nov/01/2018), one for the high-water season (May/25/2018) and one for the receding season (Aug/13/2018). Randomly selected 100,000 samples were used to compute R_rs for both sensors. These R_rs were first used to compute the Modified Normalized Water Index [100] (MNDWI: (R_rs,green − R_rs,SWIR)/(R_rs,green + R_rs,SWIR)), setting a threshold of 0.3 to build a water mask. Besides the statistical indexes MAPE and R², the median difference (MD) between OLI and MSI TSS retrievals was also calculated using OLI as a reference [47].

3. Results

3.1. OAC Concentrations and Field R_rs Spectra

The OAC concentration shows a high dependency on the hydrograph phase and water level (

Table 2. The mean, minimum (Min), maximum (Max), and standard deviation (SD) values for TSS, TSI, Total Suspended Organic Solids (TSO), and Chl-a concentration and the mean TSI/TSS ratio for each campaign.

Campaign	Sample Size	Statistics	TSS	TSI	TSO	Chl-a	Mean TSI/TSS
			(mg L−1)	(mg L−1)	(mg L−1)	(µgL−1)	(%)
June/15	31	Mean	16.76	10.08	6.7	3.23	55.93
		Min	5.25	2	3.25	0.35
		Max	33.2	27.6	12.75	8.59
		SD	7.15	6.24	2.58	1.97
March/16	18	Mean	74.46	57.85	16.61	19.76	71.58
		Min	20.5	10.5	10	0.66
		Max	235.5	218.5	27	42.32
		SD	51.16	49.06	4.82	12.75
July/16	27	Mean	32.86	21.93	10.93	21.97	64.40
		Min	17	6.86	3	5.21
		Max	64.67	56.67	21.25	85.05
		SD	11.93	11.8	5.44	20.89
August/17	18	Mean	17.7	11.12	6.59	31.11	59.12
		Min	7	3.38	3.62	9.34
		Max	43.5	33.75	12.12	67.85
		SD	8.52	7.53	2.39	17.24
All Dataset	94	Mean	32.61	22.83	9.79	17.12	62.97
		Min	5.25	2	3	0.35
		Max	235.5	218.5	27	85.05
		SD	31.84	28.64	5.5	17.83

). For June 2015 and August 2017, similar mean values of TSS, TSI, and TSO were observed. However, the Chl-a concentration varied significantly among these campaigns, from 3.23 μgL⁻¹ in June 2015 to 31.11 μgL⁻¹ in August 2017, which presented a mean value 10 times higher. On the other hand, the 2016 field campaigns (March and June) presented higher mean TSS and TSI concentrations in relation to the remaining campaigns. In addition, the amplitude values were much higher in the 2016 campaigns, with a maximum concentration for TSS of about 215 mgL⁻¹ and for TSI of 208 mgL⁻¹ in the March 2016 campaign.

The R_rs assessment (Figure 4), as expected, also indicates a dependency on OAC variability along the hydrological year. The highest R_rs were observed in March 2016 (Figure 4b) due to the highest TSS and TSI concentrations. The R_rs values in July 2016 (Figure 4c) were higher than those of June 2015 (Figure 4a) and August 2017 (Figure 4d) due to the higher TSS and TSI concentrations. Moreover, one can note the Chl-a maximum absorption at approx. 670 nm and the phytoplankton scattering at approx. 700 nm in the July 2016 and August 2017 campaigns in agreement to the high Chl-a concentration (Table 2). The June 2015 R_rs spectra have a more homogeneous shape and amplitude due to the lower amplitude in TSS and Chl-a concentrations.

3.2. Monte Carlo Simulation Results

3.2.1. All Dataset

The results of the Monte Carlo simulation using the All Dataset (

Table 3. The validation results obtained through a Monte Carlo simulation for the linear, log-transformed, and band-ratio models for both OLI and MSI sensors: The table colors refer to the model performance and becomes better towards a greener hue and worse towards red. Number of validation samples: 28.

Sensor	Model	MAPE (%)		R²		RMSE (mgL−1)
		TSS	TSI	TSS	TSI	TSS	TSI
Landsat 8 OLI	ln (Rrs_sim_OB3)	27.79	40.05	0.79	0.75	12.67	12.19
	ln (Rrs,sim_OB4)	24.81	31.93	0.83	0.82	11.72	10.85
	Rrs,sim_OB4	31.86	56.66	0.91	0.88	11.35	11.22
	ln (Rrs,sim_OB5)	26.72	44.1	0.89	0.89	10.84	10.45
	Rrs,sim_OB5/Rrs,sim_OB2	35.67	62.94	0.89	0.9	14.73	13.37
	Rrs,sim_OB5/Rrs,sim_OB3	41.27	67.02	0.82	0.81	15.05	12.94
Sentinel 2 MSI	ln (Rrs,sim_MB3)	27.41	41.61	0.79	0.75	12.85	12.27
	ln (Rrs,sim_MB4)	25.67	32.2	0.84	0.83	11.64	10.72
	ln (Rrs,sim_MB5)	21.55	30.69	0.87	0.86	10.69	9.72
	Rrs,sim_MB6	30.7	57.85	0.91	0.87	9.75	10.63
	ln (Rrs,sim_MB6)	24.15	37.02	0.91	0.91	9.46	9.21
	Rrs,sim_MB7	29.91	55.51	0.91	0.87	9.43	10.33
	ln (Rrs,sim_MB7)	23.77	36.86	0.91	0.9	9.21	9.04
	Rrs,sim_MB8	27.75	57.73	0.91	0.88	9.81	10.3
	ln (Rrs,sim_MB8)	24.86	39.61	0.91	0.91	9.44	9.25
	Rrs,sim_MB8A	33.06	57.22	0.91	0.88	10.75	10.29
	ln (Rrs,sim_MB8A)	30.15	48.13	0.89	0.89	10.75	10.29
	Rrs,sim_MB5/Rrs,sim_MB2	67.25	100.23	0.12	0.12	20.54	17.12
	Rrs,sim_MB5/Rrs,sim_MB3	41.43	54.31	0.37	0.38	20.72	34.37
	Rrs,sim_MB6/Rrs,sim_MB2	31.64	43.3	0.79	0.79	19.49	15.78
	Rrs,sim_MB6/Rrs,sim_MB3	32.91	45.69	0.88	0.89	12.7	11.6
	Rrs,sim_MB7/Rrs,sim_MB2	32.78	43.44	0.82	0.82	18.94	15.32
	Rrs,sim_MB7/Rrs,sim_MB3	35.89	47.73	0.85	0.85	11.45	10.38
	Rrs,sim_MB7/Rrs,sim_MB4	53.54	90.83	0.57	0.45	33.51	15.77
	Rrs,sim_MB8/Rrs,sim_MB2	32.42	47.93	0.85	0.86	16.25	14.55
	Rrs,sim_MB8/Rrs,sim_MB3	37.12	53.38	0.88	0.87	12.29	10.73
	Rrs,sim_MB8A/Rrs,sim_MB3	40.40	69.27	0.83	0.83	15.12	12.53

) indicated that for OLI ln (R_{rs,sim_OB4}) models (MAPE = 24.8% for TSS and 31.93% for TSI) were the best, while the best MSI models were ln (R_{rs,sim_MB5}) (MAPE = 21.55% for TSS and 30.69% for TSI). Despite the similar MAPE results for TSS, ln (R_{rs,sim_OB3}) and ln (R_{rs,sim_OB5}), the latter provided the highest R² (0.91). The performance of green, red, and NIR bands was similar for both sensors, which was expected given the similarity of their spectral response function. The remaining NIR bands (R_{rs,sim_MB6}, R_{rs,sim_MB7}, and R_{rs,sim_MB8}) also produced accurate results for TSS (MAPE < 30%). However, regarding TSI models, those MSI bands did not perform well (MAPE > 30%). In summary, the log-transformed models (MAPE < 30.2% and RMSE < 13 mgL⁻¹) outperformed the band-ratio models (MAPE > 30% and a maximum value of RMSE of 33.51 mgL⁻¹), despite the occasional high R² values (e.g., R_{rs,sim_OB5}/R_{rs,sim_OB2}, R² > 0.88).

3.2.2. Campaign Dataset Models

In comparison with the All Dataset models, the log-transformed models for the Campaign Dataset provided better results for TSS (Figure 5a) and TSI (Figure 5b), indicating that campaign-filtered models perform better than models using data from the All Dataset (Table 3). However, the overall performance depends on both the sediment and Chl-a concentration in a given hydrograph state. At a given season and band, the All Dataset outperformed the Campaign Dataset (e.g., the ln (R_{rs,sim_MB3}) and ln (R_{rs,sim_OB3}) models for the March/16/16 campaign were worse than for the All Dataset). The TSS and TSI log models based on the green spectral band had poor results when using data acquired in March 2016 (Rising Water) (MAPE values > 33%) and in July 2016 (Receding Water) (MAPE > 39%) for both sensors. The TSS models for these campaigns performed better using higher wavelengths due to an R_rs saturation of the visible bands with increasing particle concentration (See Table 2). On the other hand, in June 2015 (High Water) and August 2017 (Receding Water), the campaigns MAPE values tended to be lower for models based on VIS bands (MAPE < 21% and < 33% for TSS and TSI, respectively) than that of NIR band (MAPE > 40% for the ln (R_{rs,sim_OB5}) model). These campaigns presented a lower sediment concentration in relation to the 2016 campaigns, which implies in a low signal in the NIR bands resulting in better results for the visible spectrum part.

In relation to the band-radio models (Figure 6), MAPE values (Figure 6a) lower than 40% were observed for most of the TSS models. However, the band ratios R_{rs,sim_OB5}/R_{rs,sim_OB3,} and R_{rs,sim_MB8A}/R_{rs,sim_MB3} for June 2015 and the ratio R_{rs,sim_MB5}/R_{rs,sim_MB2} for March 2016 resulted in MAPE > 40%. Splitting the data according to the hydrograph states resulted in TSS models more accurately tuned than those provided by the All Dataset for most of the models and campaigns, except for the model based on the 850 nm band using data acquired in June 2015. It also happened in relation to the TSI models (Figure 6b) which presented a MAPE < 50% for the models using the 850 nm band for the June 2015 and August 2017 campaigns. Moreover, the models using data acquired in August 2017 using MSI NIR bands 6 and 7 (740 and 760 nm) also presented MAPE values higher than 50%. The best results for the TSS and TSI band-ratio models were provided by R_{rs,sim_MB5}/ R_{rs,sim_MB3} for both TSS (MAPE < 25%) and TSI (MAPE < 33%) in all the campaigns. Furthermore, for the July 2016 campaign, all the models provided MAPE values smaller than 30% for TSS and TSI with a lower wavelength dependence.

Differently from the All Dataset models, the TSS and TSI models based on the Campaign Dataset showed the lowest MAPE values with a band ratio except for the TSI log model (R_{rs,sim_OB4}, MAPE = 26.89%) for August 2017.

3.2.3. Lake Dataset Models

The Monte Carlo simulation for the Lake Dataset models outperformed the All Dataset models (Figure 7). The Monte Alegre TSS (Figure 7a) models have a MAPE lower than 20% independent of the model type, with the best performance for the red band of both sensors (MAPE < 10%, R² > 0.96). The TSI models (Figure 7b) using the Monte Alegre dataset also resulted in a MAPE < 20%, except for the green band (MAPE approx. 25%). However, the model based on the Curuai Lake dataset resulted in a MAPE higher than that of Monte Alegre, with a similar performance to the All Dataset models. The best model for TSS (ln (R_{rs,sim_MB5})) resulted in MAPE = 21.55%. The 850 nm band, independent of the sensor, resulted in the poorest models for estimating TSI, (MAPE > 40%). TSI errors, however, were much greater than those of the TSS models, with the best performance (MAPE = 28.07%) for TSI by the Curai Lake Dataset model ln (R_{rs,sim_MB5}).

Regarding band ratios, the Lake models provided accurate results for both TSS and TSI when using the Monte Alegre (Figure 8a,b) Dataset (MAPE < 30% for most models). However, for Curuai Lake, the MAPE was always higher than 30%, with the best performance of R_{rs,sim_OB5}/R_{rs,sim_OB2} (MAPE = 32.25%) and R_{rs,sim_MB8}/R_{rs,sim_MB2} (MAPE = 47.13%) for TSS and TSI, respectively. Both the log and band-ratio models provided MAPE results similar to those of the models based on the All Dataset models.

The log-models presented MAPEs smaller than those of the band-ratio models for both the Curai and Monte Alegre lakes. For Curuai, the best model was obtained using ln (R_{rs,sim_MB5}) for TSS (MAPE = 21.6%) and TSI (MAPE = 28.07%), while for Monte Alegre, the best was ln (R_{rs,sim_MB4}) (MAPE = 9.93%) for TSS and ln (R_{rs,sim_MB8}) (MAPE = 13.66%) for TSI.

3.3. Image Application

3.3.1. Atmospheric Correction and Glint Removal

The performance of the atmospheric correction before and after glint correction differed from OLI to MSI sensors (Figure 9). For MSI, the results showed that glint influence was high in all spectral bands, with a MAPE higher than 119% when assessing the entire set of VNIR spectral bands (Figure 9a,c,e). A MAPE reduction of over 100% was observed in all spectral bands after glint correction, independent of the atmosphere correction method. Furthermore, glint correction resulted in increases in R². The best results were observed when using the 6S model (Figure 9a) after glint correction (R_{rs,MSI_6S_deglint}) (MAPE = 14.35% and R² = 0.92) for all VNIR bands. Glint correction (Figure 10) when analyzing individual MSI spectral bands (Figure 10a) produced good results, mainly in the NIR bands, with a MAPE reduction of up to 150%. The Sen2Cor atmospheric correction performed better in the blue (B2), green (B3), and red-edge (B5) bands. However, Sen2Cor had poor performances in the NIR and red bands. The red and NIR bands at 783 and 842 (B4, B7, and B8) presented better results when using the 6S method. ACOLITE presented slightly better results than 6S (a difference of about 0.05% in MAPE values) in the 740 nm band (B6). After glint correction, OLI atmosphere correction (Figure 9b,d,f) produced MAPE < 30%. The MAPE difference between glint-corrected and not corrected images was not higher than 45% for OLI images. The best results for all VNIR bands were obtained by the L8SR (MAPE = 16.68 and R² = 0.93) atmosphere-corrected model. Glint correction also reduced MAPE values in the green (B3), red (B4), and NIR (B5) OLI bands independent of the atmospheric correction methods. However, the 6S atmospheric correction model caused an increase in MAPE after glint correction of the blue band (B2). In the green and red spectral bands, the MAPE was smaller than 15% after glint correction, whereas in the NIR band, 6S and L8SR resulted in MAPE < 23%

3.3.2. Model Application to OLI and MSI Imagery

As described in Section 2.4, the All Dataset models calibrated with in situ R_rs were applied to OLI and MSI imagery acquired on Aug/10/17 and Aug/08/17, respectively over Lago Grande de Curuai (LGC). Only the 2017 imagery were used due to proximity to the 2017 field campaign and the low cloud cover.

Table 4. The MAPE values (%) for empirical models applied to OLI atmosphere-corrected scenes after glint correction (6S, L8SR, and ACOLITE). The R² and RMSE are in the Supplementary Files C. The bold values in each column refer to the best MAPE values (n = 16).

Model	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)
ln (Rrs,OB3_deglint)	6S	18.77	19.87	L8SR	16.81	19.07	ACOLITE	17.23	18.08
ln (Rrs,OB4_deglint)		29.29	31.52		27.95	31.21		25.63	28.24
ln (Rrs,OB5_deglint)		29.55	36.13		28.99	36.08		51.46	62.50
Rrs,OB5_deglint/Rrs,OB2_deglint		494.66	947.42		102.15	161.10		238.99	387.43
Rrs,OB5_deglint/Rrs,OB3_deglint		74.57	116.95		62.02	98.17		140.38	210.00
Rrs,OB5_deglint/Rrs,OB4_deglint		165.31	264.10		166.22	262.35		271.09	456.72

presents the validation results (MAPE) using the equations with the models calibrated using the All Dataset calibrated models: ln (R_{rs,OLI-6S-deglint}), ln (R_{rs,OLI-L8SR-deglint}), and ln (R_{rs,OLI-ACO-deglint}).

The best model was ln (R_{rs,OB3_deglint}) independent of the atmospheric correction method (MAPE < 21%, R² > 0.73 for both TSS and TSI. Moreover, the red and NIR bands present accurate results for OLI. The ACOLITE atmospheric correction had the best performance (MAPE < 28% and R² > 0.74 for both TSS and TSI) for the red band models. For models using NIR bands, both the 6S and L8SR atmospheric corrections presented similar results, being slightly better for 6S (MAPE < 35% and R² > 0.83). The band ratio models, however, presented poor results, with MAPE values higher than 60% and R² < 0.44 independent of the atmosphere correction models and bands. The best TSS model was ln (R_{rs,OB3-L8SR-deglint}) (MAPE = 17.39%, R² = 0.73), while for TSI, it was ln(R_{rs,OB3-ACO-deglint}) (MAPE = 19.22, R² = 0.84).

TSS models using green MSI bands (

Table 5. The MAPE values (%) for empirical models applied to MSI atmosphere-corrected scenes after glint correction (6S, ACOLITE, and Sen2Cor). The R² and RMSE are in the Supplementary Files C. The bold values in each column refer to the best MAPE values (n = 14).

Model	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)	Atm. Cor	MAPE TSS (%)	MAPE TSI (%)
ln (Rrs,MB3_deglint)	6S	16.86	25.87	ACOLITE	18.09	28.88	Sen2Cor	23.88	21.55
ln (Rrs,MB4_deglint)		24.45	23.70		22.45	21.25		21.69	20.30
ln (Rrs,MB5_deglint)		16.67	18.63		19.42	21.83		28.12	31.03
ln (Rrs,MB6_deglint)		16.43	23.22		15.31	24.48		28.65	44.35
ln (Rrs,MB7_deglint)		20.14	25.43		30.10	34.22		31.25	42.85
ln (Rrs,MB8_deglint)		19.89	33.91		23.25	37.17		66.78	76.81
Rrs,MB5_deglint/Rrs,MB2_deglint		65.48	101.17		81.16	133.42		68.15	106.19
Rrs,MB5_deglint/Rrs,MB3_deglint		24.83	39.28		17.12	27.05		71.93	117.00
Rrs,MB5_deglint/Rrs,MB4_deglint		127.97	163.22		118.86	131.20		145.13	230.42
Rrs,MB6_deglint/Rrs,MB2_deglint		492.63	878.88		611.48	1170.11		54.67	60.89
Rrs,MB6_deglint/Rrs,MB3_deglint		30.88	45.51		25.77	42.54		76.08	82.48
Rrs,MB6_deglint/Rrs,MB4_deglint		70.91	137.23		52.55	106.52		69.24	56.20
Rrs,MB7_deglint/Rrs,MB2_deglint		77.30	89.57		219.99	293.71		78.69	86.43
Rrs,MB7_deglint/Rrs,MB3_deglint		55.04	80.17		104.83	150.12		66.25	69.67
Rrs,MB7_deglint/Rrs,MB4_deglint		111.93	136.82		169.23	184.08		64.23	48.26
Rrs,MB8_deglint/Rrs,MB2_deglint		49.19	57.61		81.74	95.00		89.97	96.65
Rrs,MB8_deglint/Rrs,MB3_deglint		42.85	56.61		50.15	66.11		93.07	96.21

) produced MAPE values similar to those of the OLI green band (MAPE < 19% for 6S and ACOLITE for MSI atmosphere-corrected scenes and MAPE < 19% for the three OLI atmosphere-corrected models). However, in the red and NIR bands, the MSI results outperformed OLI (MAPE < 25% at the red band and MAPE < 20% at NIR B8). The best results were produced by the 705 nm band (B5) and 740 nm band (B6) for TSS (MAPE < 20%) and for TSI (MAPE < 25%), using either 6S or ACOLITE atmosphere correction. For the red MSI band, Sen2Cor had the best performance; however, for all NIR bands, the poorer performance of Sen2Cor was related to errors in atmosphere correction (see Figure 10). Despite the lower MAPE produced by ln (R_{rs,MB6-6S-deglint}) and ln (R_{rs,MB6-ACO-deglint}) for TSS and TSI in relation to the ln (R_{rs,MB5-6S-deglint}), the latter had higher R² (0.69 and 0.81 for TSS and TSI, respectively). Therefore the MSI model ln (R_{rs,MB5-6S-deglint}) was selected for spatializing both TSS and TSI.

The best TSS and TSI models provided by the All Dataset were applied to OLI and MSI scenes and compared with models provided by the Campaign (August/17) Dataset and Lake (Curuai) Dataset to assess their performance. Therefore, the following models were applied to OLI images: ln (R_{rs,OB3-L8SR-deglint}) and ln (R_{rs,OB3-ACO-deglint}) for estimating TSS and TSI concentration, respectively. For MSI, the selected model was ln (R_{rs,MB5-6S-deglint}) for both variables (TSS and TSI) (

Table 6. The equations used for model validation in OLI and MSI imagery: All is All Dataset; Curuai is Lake Model for Curuai Lake between 2015–2017; 2017 is Campaign Dataset for 2017 campaign. The subscripts in the model column refers to the model used (All, 2017, or Curuai), and the atmospheric correction method, as well TSS and TSI, refer to the type of sediments.

Sensor	Model	Dataset	Equation
Landsat 8 OLI	ln (Rrs,OB3_All-L8SR) − TSS	All	ln (TSS) = 9.656 + 1.672 × ln (Rrs,OB3_L8SR_deglint)
	ln (Rrs,OB3_Curuai-L8SR) − TSS	Curuai	ln (TSS) = 9.753 + 1.721 × ln (Rrs,OB3_L8SR_deglint)
	ln (Rrs,OB3_2017-L8SR) − TSS	2017	ln (TSS) = 9.833 + 1.781 ×ln (Rrs,OB3_L8SR_deglint)
	ln (Rrs,OB3_All-ACO) − TSI	All	ln (TSI) = 10.73 + 2.08 × ln (Rrs,OB3_ACO_deglint)
	ln (Rrs,OB3_Curuai-ACO) − TSI	Curuai	ln (TSI) = 10.881 + 2.175 × ln (Rrs,OB3_ACO_deglint)
	ln (Rrs,OB3_2017-ACO) − TSI	2017	ln (TSI) = 13.305 + 2.769 × ln (Rrs,OB3_ACO_deglint)
Sentinel 2 MSI	ln (Rrs,MB5_All-6S) − TSS	All	ln (TSS) = 8.318 + 1.336 × ln (Rrs,MB5_6S_deglint)
	ln (Rrs,MB5_Curuai-6S) − TSS	Curuai	ln (TSS) = 8.129 + 1.289 × ln (Rrs,MB5_6S_deglint)
	ln (Rrs,MB5_2017-6S) − TSS	2017	ln (TSS) = 8.605 + 1.413 × ln (Rrs,MB5_6S_deglint)
	ln (Rrs,MB5_All-6S ) − TSI	All	ln (TSI) = 8.447 + 1.511 × ln (Rrs,MB5_6S_deglint)
	ln (Rrs,MB5_Curuai-6S ) − TSI	Curuai	ln (TSI) = 9.59 + 1.799 × ln (Rrs,MB5_6S_deglint)
	ln (Rrs,MB5_2017-6S ) − TSI	2017	ln (TSI) = 8.87 + 1.639 × ln (Rrs,MB5_6S_deglint)

). The results show that, for OLI images, the All Dataset models had the best performance for both TSS and TSI estimation. For the MSI images, however, Curai Lake Dataset had the best performance, whereas, for TSI, the best performance was provided by the All Dataset (Figure 11).

The scatterplot of TSS and TSI (Figure 12) for both sensors indicates a good agreement between the in situ concentration and satellites estimates. Moreover, the high R² (0 > 0.71) and low MAPE (< 21%) (see Table 4 for OLI and Table 5 for MSI values) show that the model can be applied with relatively small errors.

Moreover, the TSS and TSI results show that both models presented similar results (Figure 12). The mean TSS values applied to OLI imagery (Figure 13) at Curuai Lake was 15.47 ± 7.71 mgL⁻¹ (± refers to standard deviation), whereas for MSI imagery, it was 17.23 ± 7.53 mgL⁻¹. These values were very close to the in situ measurements acquired in the August 2017 field campaign (17.7 ± 8.52 mgL⁻¹). The TSI concentration also presents similar results when compared with the TSI measured in the August 2017 field campaign. The mean TSI estimates were 8.11 ± 5.53 mgL⁻¹ for OLI, 9.52 ± 4.93 mgL⁻¹ for MSI, and 11.12 ± 7.53 mgL⁻¹ for the in situ measurements. The difference between the mean in situ and satellite estimates of TSS and TSI were below 20%.

In order to compare the TSS estimates retrieved from the same-day OLI and MSI overpasses, ln (R_rs,OB5) and ln (R_rs,MB8) algorithms were applied to low-water period scenes to avoid saturation. For high- and receding-water periods, the ln (R_rs,OB3) and ln (R_rs,OB5) were applied. Furthermore, the ln (R_rs,OB4) and ln (R_rs,MB4) algorithms were also applied to evaluate the performance using the same spectral band.

At the low-water season, the comparison between TSS values obtained for the algorithms that used the red band (Figure 14a) provided similar median values (15.65 mgL⁻¹ for OLI and 15.47 mgL⁻¹ for MSI) with a MD of −0.7 mgL⁻¹ and MAPE values of 6.66%, indicating the convergence between these two estimates. This convergence was also observed in the comparison between ln (R_rs,OB3) and ln (R_rs,OB5) (Figure 14b), with median values of TSS close to those observed for the red band (13.86 and 14.37 mgL⁻¹ for OLI and MSI, respectively).

Regarding the algorithms applied to the low-water season (Figure 14c), the results demonstrated, again, the convergence between TSS values for both sensors. The MAPE value was smaller than 8%, demonstrating the seamlessness of the TSS product. When all dates were compared, the results, one more time, demonstrated the convergence between the TSS products, with MD values of −2.54 mgL⁻¹ for a broad range of TSS concentration (up to 250 mgL⁻¹, with median values of 28.13 and 30.47 mgL⁻¹ for OLI and MSI, respectively).

4. Discussion

4.1. Variability of OACs and R_rs during Field Campaigns

The TSS, TSI, and Chl-a concentrations and R_rs spectra (Figure 4 and Table 2) showed a high dependence on the hydrograph phase and climatic events, as pointed out in studies in the region [8,16]. Low TSS and TSI concentrations in the June 2015 (High period) and August 2017 (Receding period) field campaigns implied low-intensity R_rs spectra. In contrast, March 2016 (Rising period) presents high TSS and TSI values responsible for high-intensity R_rs spectra. R_rs also merged into two different R_rs spectra clusters (Figure 4b) due to the high variability in amplitude and proportion of OACs in this campaign (See Table 2). Furthermore, the mean TSI/TSS ratio presented high values (approx. 72%), indicating a high concentration of inorganic sediments related to the input of the Amazon River sediments through river channels at Rising Water [7,8]. The inorganic particle inflow increases R_rs values due to its high refraction index [54]_, which increases particle backscattering [101,102]. Higher TSS and TSI values were also observed in the July 2016 (Receding Water level) campaign compared to those measured during the August 2017 receding water level campaign. The interannual differences in concentration might be related to the 2016 drought associated to El Niño anomalies [103]. Moreover, the mean TSI/TSS ratio in July 2016 was higher (approx. 64%) than that observed in August 2017 (approx. 59.12%) (See Table 2). The water level in August 2017 was close to the historical mean, causing TSS and TSI to increase as the water level decreased [8].

In relation to the Chl-a concentration, the values were very small in June 2015 (3.23 µgL⁻¹) and did not imprint its absorption bands in R_rs (Figure 4a). The low Chl-a concentration in this phase is related to the high water level of the Amazon River (See Figure 2) [22] which increases the Amazon inflow into the floodplain lakes. This high inflow increases the turbulence, preventing phytoplankton growth [60,104]. Furthermore, the overbank flow causes the input of Colored Dissolved Organic Material (CDOM), reducing light availability in the blue region within the water column which reduces the rates of phytoplankton growth [7,8,13]. The Chl-a concentration in August 2017 was the highest followed by July 2016. These high values at the beginning of the Receding Water Level Period (July 2016) and in the Receding Period (August 2017) are attributed to the higher nutrient supply to the lakes during rising and high water periods [13,105,106]. Furthermore, in the receding period, the inflow of CDOM and inorganic matter from the Amazon River to the floodplain decreases, resulting in more light availability, which favors phytoplankton growth [60,104].

4.2. In Situ Model Accuracy and Relation to OACs Variability

The Monte Carlo simulation results using the All Dataset (Table 3) provided better results for the TSS log models (MAPE < 28%) than those of the linear and band-ratio models (MAPE > 30% for all models evaluated). This poor performance of band ratios may be related to the concentration ranges of TSS and TSI provided by the All Dataset. The mean TSS and TSI values were 32.61 mgL⁻¹ and 22.83 mgL⁻¹, respectively for the All Dataset. These concentration, associated with a high sediment concentration variability, might affect the accuracy of band-ratio models since the literature reports that they do not work well for TSS values lower than 10 mgL⁻¹ [92,107,108]_. Log TSS models were accurate, with MAPE values below 28% for both sensors (except for ln (R_{rs,sim_MB8A}) (MAPE = 30.15%), with the best results provided by ln (R_{rs,sim_OB4}) (MAPE = 24.81%, R² = 0.83) and ln (R_{rs,sim_MB5}) (MAPE = 21.55%, R² = 0.87) for OLI and MSI, respectively (Table 3). The Ln (R_{rs,sim_MB5}) best results for TSS models could be attributed to the lack of a Chl-a absorption signal at 705 nm and the sensitivity of this spectral band to both phytoplankton and inorganic material backscattering [109]. The green band models displayed a lower accuracy mainly related to R_rs saturation at higher TSS and TSI concentrations at these spectral bands [38,54,62,92]. The best TSI models were also provided by ln (R_{rs,sim_OB4}) (MAPE = 31.93%, R² = 0.82) and ln(R_{rs,sim_MB5}) (MAPE 30.68%, R² = 0.86). These higher errors for TSI models, when compared to TSS models, are related to backscattering by organic particles that, even with less intensity [52,101], affected the results [102,110]. The results show a good convergence in the performance of both sensors in similar bands (See Table 3), which would help to create virtual constellations as they provide similar errors when referred to in situ R_{rs_sim} bands [47,81]. Zheng et al. [111] using R_rs,sim for the OLI NIR band also obtained good results for TSS estimates (MAPE = 18.6%, RMSE = 5.79 mgL⁻¹) at Dongting Lake (China) for a concentration range of OACs as follows: TSS (4–101 mgL⁻¹, mean = 28.7), TSI (2.3–93.3, mean = 27.5), TSO (0.7–24.7, mean = 5.9), and Chl-a (1.9–124.2 μgL⁻¹, mean = 28.7). Liu et al. [112] using R_rs,sim for an MSI sensor reported also accurate validation results for NIR MSI bands for TSS estimates (MAPE < 20%) at Poyang Lake (China). However, as TSS concentration in that study presented higher values (mean = 76.05 mgL⁻¹), their results were slightly different from this present study regarding MAPE and correlation values following the trend of being better at lower wavelengths [38,112].

The performance of the Campaign Dataset models also varied with the hydrograph phases as expected, given the high OACs variability [8] (Table 2). The best TSS and TSI models based in the June 2015 Campaign Dataset were provided for the log models using shorter wavelengths (< 705 nm) (Figure 5). The best band-ratio model (R_{rs,sim_MB5}/R_{rs,sim_MB3} model) (Figure 6) also included short wavelengths due to the low sediment concentration observed during this campaign [38]. In relation to TSS modeling, the MAPE amplitude varied from 18.3% (R_{rs,sim_MB5}/R_{rs,sim_MB3}) to 20.5% (R_{rs,sim_MB3}) (Figure 5 and Figure 6). The TSI log models presented higher MAPE values than the TSS models mainly in the green spectral band for both sensors, with better results in the red band. This effect might be related to the TSO concentration and TSI/TSS (Ratio = 55,93%,Figure 15a) during the June 2015 campaign, probably contributing to the green band backscatter [54]. The MAPE difference between the red and red-edge bands were not higher than 5% due to the lesser influence of TSO and CDOM in these spectral regions [110,113]. Moreover, a low Chl-a concentration does not affect R_rs,sim (See Figure 4), which favors its performance for estimating TSI concentrations [54,110]. At longer wavelengths, the log models showed an increase in MAPE (up to 28.61% for TSS and 40.20% for TSI at 850 nm bands) since a low sediment concentration causes lower particle backscattering [54] relative to the high water absorption coefficient of those spectral bands [53]. The band-ratio models followed a similar relationship with wavelengths as log models with longer wavelengths presented the worst results (e.g., R_{rs,sim_MB8}/R_{rs,sim_MB4}, MAPE = 41.07% for TSS and 52.93% for TSI) due to the low R_rs sensibility with low sediment concentrations at NIR bands [108]. For the August 2017 campaign, the TSS and TSI values were similar to those of June 2015. However, due to a high Chl-a concentration (31.11 µgL⁻¹), the best models were those based on green (for both sensors) and red-edge bands (B5 for MSI with MAPE < 20%). In the August 2017 campaign, both phytoplankton and inorganic particles increased the backscatter at those spectral bands [109,114,115]. The R_{rs,sim_MB5/}R_{rs,sim_MB3} band-ratio model also presents good results (MAPE = 15.4%). For TSI models, the best results were provided by the red band for both sensors (MAPE < 28%) mainly due to the lower backscattering effects of phytoplankton cells [52,54] at the red band, which explains also the higher errors in the green and red-edge (MSI-B5) bands. The band-ratio R_{rs,sim_MB5/}/R_{rs,sim_MB3} also presented the best results for TSI models. Similar to the June 2015 campaign, MAPE tended to increase with wavelength, mainly in campaigns characterized by low TSS and TSI concentrations. Campaigns with high TSS and TSI concentrations (March 2016 and July 2016) showed the opposite trend. High MAPE values were obtained for green and red bands in both sensors due to reflectance saturation in these spectral bands in the presence of a high sediment concentration [63] and a high proportion of TSI in relation to TSS, which increases backscattering due to the higher refraction index of inorganic particles [54,110] (TSI/TSS ratio = 71.48% and 64.40% for March 2016 and July 2016, respectively—See Figure 15a). On these campaigns, the MAPE values tend to decrease with increasing wavelengths for both TSS and TSI as observed in several studies in aquatic systems characterized by a high sediment concentration [62,63,64,112]. In these campaigns, band-ratio TSS models also showed good results, mainly in July 2016, with wavelengths higher than 740 nm (MAPE < 16%).

Regarding the Lake Dataset models for TSS and TSI, Monte Alegre presented the best performance. The reason for that might be the homogeneity of this aquatic system (Chl-a concentration = 6.74 μL⁻¹ and 21.76 μL⁻¹ for Monte Alegre and Curuai, respectively), which minimizes the Chl-a influence on another OAC estimate. Monte Alegre MAPE results were close to each other, with higher values in the green band of both sensors due to the saturation of this spectral band [38]. The band-ratio models for Monte Alegre also presented good results (MAPE < 40%). For Curuai Dataset models, the MAPE values were much higher than those of Monte Alegre due to the high variability in OACs along the hydrograph phases. However, the MAPE values were not higher than 30% for all TSS models, with the R_{rs,sim_MB5} model displaying the best results when using this dataset. Differently from the Monte Alegre Lake models, the MAPE difference between TSS and TSI using models based in Curuai Lake Dataset increases up to 10% in the green-NIR interval due to the influence of organic matter, related to a smaller TSI/TSS ratio (Figure 15b). Band-ratio models also present better results for Monte Alegre than for Curuai, however, with MAPE values higher than those of the log models. For Curuai, the band-ratio MAPE values were close to those of the models based on the All Dataset.

4.3. Influence of Atmosphere and Glint Correction on R_rs Retrieval

Atmospheric and glint correction for both sensors presented an excellent performance (See Figure 9). The results for OLI were better for green and red bands, with MAPE values lower than 15% and R² > 0.83 independently of the atmospheric correction method, provided that glint correction was applied. For NIR band, despite the small MAPE (<23%) for 6S and L8SR, R² was lower than 0.6 for both methods, indicating that the accurate retrieval of R_rs values in high wavelengths [21,56,89] is still a challenge as none of the atmospheric correction methods provide a suitable R_rs estimate in such a band. Good R_rs validation results were also reported by Lobo et al. [19] for the green and red bands of TM (Landsat-5) submitted to 6S atmosphere and glint correction (R² > 0.9) at Tapajós River (Brazil) and by Bernardo et al. [69] using OLI L8SR in an eutrophic Brazilian reservoir without glint correction for green (MAPE = 6.39%) and red (MAPE = 16.7%) bands. Furthermore, Bernardo et al. [69] also obtained poor results for the NIR OLI band using ACOLITE and L8SR (MAPE = 86.99% and 72.8%, respectively). The performance of L8SR in green bands was also reported by Wei et al. [77] at Harbor Bay turbid waters (TSS of approx. 3–8 mgL⁻¹ and Chl-a of approx. 2–5 µgL⁻¹) with MAPE values of 18%.

Atmospherically corrected MSI images, disregarding the method, resulted in a high MAPE (>100%) before glint correction. Two reasons may explain that fact: (i) The orbit issue related to the near-nadir viewing angles of the MSI sensor [67], which contribute to glint signal receiving at the satellite level, and due to the time of satellite overpass (approx. 20 minutes after OLI) [67,90], which implies higher solar elevation angles and an increase in glint effects; (ii) the environmental conditions at MSI passage, as the mean wind speed in August 2008 (MSI satellite) was 3.99 ms⁻¹, larger than that of August 2010 (OLI satellite) that was 1.62 ms⁻¹. This difference in wind speed in shallow lakes [116] might be responsible for the differences in water surface roughness, increasing the probability of sun and sky glint [87]. In that sense, the performance of MSI atmosphere correction was only acceptable after glint correction as MAPE values improved up to 100%. The best performance was obtained with the 6S corrected red band (MAPE = 6.18%). For the remaining visible bands, Sen2Cor and 6S provided results more accurate than ACOLITE. Martins et al. [66] reported similar results for ACOLITE atmospheric correction over lakes of an upper reach of the Amazon Basin (Mamirauá Sustainable Development Reserve). For NIR bands, 6S and ACOLITE produced smaller MAPE (< 30%). However, as observed for the OLI bands, R² values decrease with increasing wavelengths. R² values were acceptable only for the red-edge B5 (approx. 0.75 for all atmospheric correction methods) of NIR bands. These results are different from those of Martins et al. [66], who reported the highest accuracy with Sen2Cor in NIR bands. The poorest results for the NIR MSI bands using Sen2Cor were observed in other studies [117,118], and this difference was related to uncertainties in aerosol optical thickness because Sen2Cor is more sensitive to this atmospheric correction variable [41,117,119].

4.4. Validation of Monte Carlo Calibrated Models on Corrected OLI and MSI Imagery

The accuracy obtained for OLI and MSI sensors provided accurate estimates of TSS and TSI over the Curuai Lake floodplain system. The advantage of the refined spatial resolution even allowed the quantification and distinguishing between sediment plumes coming from narrow entrance channels and overbank flow. Most of the MODIS 250 m algorithms for the Amazon River sediment estimates [42,44,120] are suitable for grasping the low longitudinal variability in sediment concentration (approx. 17mgL⁻¹ at each 100 km) [121,122]. Curuai Lake, however, demands a higher spatial resolution to resolve the small-scale spatial gradients of TSS and TSI surface concentrations related to the high variability of sediment concentration derived from environmental forces such as wind-induced sediments resuspension [7,8], variability in channel depth [123,124], and water influxes from the Amazon River [16]. Previous studies reported the problems related to small-scale spatial heterogeneity in water optical proprieties [125] and in sediment plumes [45,126], which have indicated that coarse spatial resolution underestimates sediment concentrations up to seven times [45,126].

Regarding the three different datasets used (All Dataset, Campaign Dataset, and Lake Dataset models) the best results for both TSS (MAPE = 17.39%, R² = 0.74) and TSI (MAPE = 27.88%, R² = 0.89) were obtained for the All Dataset models (Table 4). These best results for the green OLI band were attributed to the range of OAC concentration in the August 2017 field campaign (e.g., 7–43.5 mgL⁻¹ for TSS, See Table 2). Due to the lower TSS concentration range, there was no saturation in the green band (Pearson coefficient between OLI B3 and OLI B4 = 0.97). Novoa et al. [63] discuss that saturation occurs when a log-modeled curve between green and red reflectance changes the linear behavior and starts to bend when the saturation starts. This effect does not occur in OLI bands, even with concentrations higher than those observed by Novoa et al. [63], probably due to the compensation provided by the absorption by CDOM and suspended organic matter. This factor favors the use of the green band for estimating sediments avoiding Chl-a absorption in the red band. Other authors also found good results for TSS estimation using the green OLI band. Bernardo et al. [69] obtained the best results for TSS (MAPE = 10.06%) estimates using the green OLI band in a Brazilian eutrophic reservoir (Chl-a concentration: 263.2–726.5 µL⁻¹) with similar TSS values within this study. Min et al. [127] uses the green TM (Landsat-5) band in South Korea coastal waters for estimating TSS with a concentration lower than 100 mgL⁻¹ and presents also good validation results (R² = 0.72). The performance of the green OLI band also follows results for ln (R_{rs,sim_OB3}) based on Campaign-Dataset Models for August 2017, with the best results for the green band model (see Figure 5). Beside the best results for the green model, the red and NIR results for TSS and TSI also present satisfactory validation results with MAPE values lower than 35% and R² > 0.82 for 6S and L8SR. As discussed in the R_rs,sim models, the worst results for the red band are attributed to Chl-a absorption. In the NIR band, higher errors are due to the low concentration of suspended sediments [38,64] and errors in atmospheric correction (R² < 0.6 for NIR OLI bands). Higher errors for TSS estimates below 50 mgL⁻¹ using NIR OLI band were also reported by Yepez et al. [64] at Orinoco River (Venezuela). With higher TSS concentrations (>50 mgL⁻¹), NIR bands present a good performance, as shown by Yepez et al. [64] (R² = 0.91 for a range 0 ~ 250 mgL⁻¹). Zheng et al. [111] also demonstrated the good performance of the NIR OLI band for TSS (varying between 4–101 mgL⁻¹) estimation (MAPE = 20.6%) due to a higher TSS concentration when compared with that of the present study. OLI band-ratio models do not present satisfactory validation results (MAPE > 60% for all models), which could be attributed to errors in the atmospheric correction of the NIR band, as discussed before and to the Signal to Noise Ratio (SNR) since band-ratio models tend to amplify noise effects in TSS models, mainly at low concentrations [98]. The results indicated that, when the three evaluated models were compared, the best performance was provided for the All Dataset-based models (MAPE < 20% for TSS and TSI). As the calibration/validation process is performed through a Monte Carlo Simulation, the randomness allows encompassing the entire OACs variability, which makes the models more robust than the remaining validation calibration/validation approaches. Moreover, due to the small sampling size of the dataset acquired in one campaign, water body properties become underrepresented.

The higher spectral resolution provided by an MSI sensor in the NIR (5 bands) in relation to the OLI increases the potential of the MSI sensor for estimating TSS and TSI even with relatively lower sediment concentrations, since the bands centered at 705, 740, and 780 nm are centered in wavelengths where the water absorption coefficient is lower than that of 800 nm [128]. Furthermore, at these bands, the influence of the remaining OACs (mainly at wavelengths > 740 nm) can be assumed negligible [54]. The best result for both sediment types was for ln (R_{rs,MB5_6S}). Regarding the three evaluated models (All, Campaign, and Lake Dataset models), the Lake Dataset performed better for TSS (MAPE = 15.70%, R² = 0.69) and the All Dataset model performed better for TSI (MAPE = 18.63%, R² = 0.81), using 6S after glint correction. These results agree with those obtained with the in situ-simulated R_rs for the August 2017 field campaign, which presented the best performance for the ln (R_{rs,sim_MB5}) model. For band-ratio models, the results for MSI do not present satisfactory R² values (<0.52), and these poorest results could be attributed to the relatively low SNR for the MSI sensor that intensify errors in band-ratio models [98]. However, despite the low accuracy of MSI band-ratio models, these are better than those obtained for the OLI bands which are related to lower wavelengths of MSI NIR bands (e.g., the ratio between MSI band 5 (approx. 705 nm) by MSI band 3 (approx. 550 nm) resulted in MAPE = 17.1%, and R² = 0.52). The MSI results were better than those reported by Liu et al. [112] for a higher TSS amplitude (up to 300 mgL⁻¹). Besides that, the authors found good validation results for MSI B5 (MAPE = 34.56%, R² = 0.81), showing the possibility of using this band with higher concentrations than those of the present study. Olmanson et al. [109] also found good results (R² =0.81) at 705 nm for TSS estimates using a hyperspectral sensor at the Mississippi River (EUA) (TSS varying between 4–80 mgL⁻¹).

The results from the comparison between the concurrent overpasses of both satellites indicate the congruence of their TSS estimates. The median difference in TSS concentration between both sensors was smaller than 0.7 mgL⁻¹ during the high and receding water seasons, with slightly higher values for OLI. Regarding lower concentrations (ranging from 0 to up 8 mgL⁻¹), Pahlevan et al. [47] observed higher TSS estimates when using an OLI sensor compared to MSI and reported a median difference between both sensors of −0.27 mgL⁻¹ over nine lakes around the globe.

During the low-water season in the Amazon floodplain lakes, the median difference was higher than that observed in the high water, which is related to the higher sediment concentrations, which increases TSS values. Besides that, MAPE values still demonstrate the consistency between the estimates (MAPE = 7.61%). These errors were also close to those reported for vicariously calibrated R_rs (<6% in visible bands) by Pahlevan et al. [47].

5. Conclusions

This study presents, for the first time, the performance of empirical algorithms developed using simulated Landsat-8/OLI and Sentinel-2A/MSI in situ R_rs in order to quantify Total Suspended Solids (TSS) and Inorganic Suspended Solids (TSI) concentrations over optically complex waters of lower Amazon floodplain lakes. Through a robust calibration/validation process, this study produced accurate TSS and TSI products in Amazon floodplain Lakes using remote sensing images. Furthermore, the accuracy of atmospheric and glint correction methods for OLI and MSI imagery on the floodplains lakes and their impact on TSS and TSI retrieval were also quantified.

The results support that glint correction improves the accuracy of R_rs retrieval for both sensors. Glint effects are known to be higher in MSI due to the orbital configuration; however, in the present study, environmental effects (wind-induced waves) could also be contributed to the signal degradation. For OLI, however, glint effects had smaller impacts on R_rs retrieval. For example, glint plus 6S atmospheric correction applied to MSI reduced MAPE values up to 100%, whereas for OLI, this reduction was not as high (approximately a 30% reduction in MAPE values) disregarding atmospheric correction adopted.

Empirical models provided satellite-derived TSS and TSI concentrations with a high accuracy (MAPE < 20%) based on log models using OLI-green and MSI-red-edge (B5) bands. This study result indicates that Landsat-8/OLI and Sentinel-2/MSI provide similar errors when compared to in situ TSS and TSI concentrations. The congruence analyses demonstrated that OLI and MSI can estimate TSS concentrations with MAPE values lower than 8% for the high, receding, and low water seasons when near-simultaneous satellite overpasses were compared. However, the unavailability of a high number of satellite near-simultaneous overpasses and in situ data matchups demands creativity in order to assess the accuracy and congruence between OLI and MSI products. The results demonstrated that both sensors could be used together as a virtual constellation in optically complex lake water. The congruence assessment indicated MD values lower than 3 mgL⁻¹ for all dates, throughout the flood pulse, which is sufficient for the systematic monitoring of TSS variability in floodplain lakes. Therefore, this virtual constellation would provide accurate data at a revisiting time of up to 2.9 days. profiting from Sentinel 2A–2B and Landsat 8 [48]. This revisiting time would allow increasing the acquisition probability of cloud-free images mainly during the rainy season, when most of the sediment inputs to the floodplain lakes start in February in the upper Amazon floodplain lakes, which is fundamental for uncovering the river-floodplain interaction processes on a regional scale.