The Parana River is the second longest river in South America, running through around 4000 km [36]. It is the natural boundary of multiple provinces in Argentina, reaching the Río de la Plata into its exit in the Atlantic Ocean. Paraguay River, with 2550 km [36], is a tributary of Parana River in its middle basin. The Bermejo River, an Andean tributary [37], is the primary sediment source in the Parana-Paraguay confluence. Due to the high solids presence in the Paraguay River, its discharge into the Parana River alters the characteristics of its composition, creating two distinct regions of high (West) and low (East) sediment concentration [38].

The Metropolitan Area of Gran Resistencia (MAGR) is an urban region in Chaco Province, North-East Argentina. It comprises four cities, including Resistencia, the capital city of Chaco. According to the last census, MAGR has a population of 423,000 inhabitants [39]. Parana River significantly impacts the region's fishing industry, tourism, recreational activities of the local communities, and transportation routes [40]. The water source for the MAGR potabilisation plant is located in the Barranqueras River (an arm of the Parana River), which is connected to two main rivers in the metropolitan area, the Black and Tragadero Rivers.

Figure 1 shows a map of the region of interest. The inset image corresponds to Argentina with Chaco province (pink), MAGR location (white dot) and Parana River extension (blue line). The main image is a real-colour satellite view of the study area, with the potabilisation treatment plant (yellow triangle) in Barranqueras city and the main rivers.

Figure 1.

Region of study, indicating main rivers, water plant location and sample site; in the inset, Chaco Province in Argentina, MAGR location and Parana River extension

The sample point (red star), located at 58°54'23''W 27°28'20''S, was selected over the Barranqueras River at the plant's inlet position. The water from the inlet point is pumped into a chamber and distributed to the different plant sections. Samples collected from the chamber are delivered to the in-site laboratory to measure a series of parameters, mainly turbidity, pH, electrical conductivity, and alkalinity.

Daily measurements were performed by the in-site laboratory at the water treatment plant [41] in Barranqueras City from 2017–01–01 to 2021–09–03. In this period, 1732 observations were recorded. The parameters and their units were: pH; electrical conductivity in micro siemens per centimetre [μS/cm], alkalinity in parts per million of calcium carbonate [ppm CaCO₃], and turbidity in nephelometric turbidity units [NTU]. Alongside these data, supplementary water samples were taken to assess more sediments-related parameters, such as total suspended matter (TSM), total dissolved matter (TDM), and total matter (TM). These parameters, measured in parts per million [ppm], are related since TM is the sum of TSM and TDM. One-litre water samples from the distribution chamber were stored in dark glass bottles. In total, 28 complementary samples were collected from 2021–08–24 to 2022–12–07.

The physicochemical methods applied to measure pH, conductivity, alkalinity, turbidity, TSM and TM were 4500 H⁺, 2510–B, 2320–B, 2130 B, 2540 D and 2540 B, respectively, according to Standard Methods techniques [42]. TDM was calculated as the difference between TM and TSM.

Satellite spectral R_S data were obtained from S2-MSI. Table 1 resumes the characteristics of two sensors since platforms' S2A and S2B products were used, featuring a maximum spatial resolution of 10 m (when available), 5 days revisit time, and 11 spectral bands. Bands B09 and B10, at 945 nm and 1373 nm, respectively, were discarded since no surface measurement was done at those wavelengths.

Copernicus Data Space Ecosystem provides complete, open, and free access to S2 products. A set of 382 images was acquired for the sample collection period. S2-MSI data at the processing level L2A are atmospherically corrected by the Sen2Cor processor [43]. The images disturbed by clouds were discarded using the quality assessment band QA60, acquired from the Sentinel-2 dataset in the Google Earth Engine platform [44]. This simple method was preferred over more complex approaches [45] since the QA60 band is a coded bit mask detecting clear skies, dark clouds, and cirrus clouds.

After removing the data from the days with clouds over the study area, 181 satellite products remained to continue the analysis. The products were cropped around the area of interest and then resampled to the uniform spatial resolution of 10 m. The R_S values were extracted using a 3×3 pixel window around the point near the plant water entrance on Barranqueras River (Figure 1). The final pixel value was the mean of the individual values in the grid.

As a preliminary step, the relationship between turbidity and R_S per band was studied by evaluating the potential impact of individual bands on the turbidity value. The target parameter in the modelling process was turbidity as a mathematical regression problem, with the spectral bands as predictors. Two main modelling methods were used: linear, with algebraic relationships between the predictors, and a tree-based machine learning RF approach. The linear modelling utilised several spectral bands and the normalised difference turbidity index NDTI, obtained by the red and green bands, B04 and B03 [46]. This index was used for water quality assessment because it is proportional to turbidity [47], according to eq. (1).

Machine learning techniques were applied to improve traditional methods for parameter retrieval [48]. RF operates by an ensemble of decision trees [49], each trained by a subset of the whole data. RF can manage many predictor variables and maintain low levels of overfitting [50], which is a negative aspect of modelling.

RF modelling used all spectral bands available (Table 1) since this method is suitable for finding non-linear relationships between multiple predictors. A tuning step improving the performance of RF was applied to obtain the best arguments, the hyperparameters, required in this model. The tuned hyperparameters were the minimum number of samples taken from the dataset to form a node in a decision tree (min_n) and the number of predictors that will be sampled (mtry). The 'trees' hyperparameter was fixed at 1000 units. Both steps of sample observations and predictor selection are random across all trees. The final turbidity estimation was an average value of all tree's estimations.

The tuning process used the racing technique [51]. It evaluated the model in a subset of resamples, obtaining the performance metrics and continuing only with the hyperparameters that showed promising results. Usually, racing techniques are faster to compute than traditional methods, such as grid search [52].

Pearson's coefficient of determination (R²) and root mean squared error (RMSE) were calculated to measure the model's performance. Equations (2) and (3) show the mathematical expressions for these metrics.

Where i represents each measurement from a total of n samples; y_i denotes the real turbidity value, x_i – the estimated value from the correspondent model, and the mean value of all y_i. A preferred model consists of a high value of R², closer to 1, and a low RSME.

Developing a linear or RF model included splitting the dataset into two parts: 75% of the samples for training and tuning the corresponding model, and the remaining 25% only for testing and finalising the model, that is, getting the last version of the model specification. This methodology followed the best practices for data modelling [53]. Since the original dataset corresponded to a time series of turbidity values, the testing split corresponded to the most recent dates. The training dataset was resampled using a 10-fold cross-validation.

The aim of performance metrics calculation in the training step was to select the best model. Following the selection, the model was evaluated using the testing dataset, with new and later observations, to calculate the final performance metrics. The estimated turbidity values were compared with the validation values in a time series plot.

Applying the selected model to S2-MSI products in the Barranqueras River for four dates enabled several maps showcasing the spatial distribution of turbidity. The spectral index MNDWI (Modified Normalised Difference Water Index) was used to mask the water from the scene [54]. An automatic method allowed the identification of the MNDWI threshold value [55].

Figure 2 shows the parameters measured in the water treatment plant as a time series plot. The number of samples (n) is shown in the top right corner of each panel. Scarce rains and drought in 2019 caused a historic low level of water in the Parana River [56]. It can be seen in the steady increase in turbidity (Figure 2a) and conductivity (Figure 2d). Conductivity from 2019 and forward started to be more dispersed than in previous years. Turbidity presented yearly cycles, with high values at the beginning of each year, between January and April-May, then followed by a low-turbidity period. The main statistical values per parameter are summarised in Table 2, which shows the mean, median, standard deviation (SD), initial and final sampling date, and number of samples (n).

Figure 2.

Time series of measured water parameters: turbidity, alkalinity, pH, and conductivity

Water properties are heavily influenced by dams' operations [57] in the North basin (south Brazil) being Itaipú dam (Paraguay-Brazil) and Yaciretá reservoir (Paraguay-Argentina), the closest to the study area. The primary source of sediments in the Parana River is the Bermejo River, a Paraguay River tributary, creating a high turbidity imbalance [58]. Due to heavy rains in Bermejo River headwaters between October and April, a high sediment concentration occurred in the Parana River between December and May [57]. Said period corresponded with the turbidity cycles observed in Figure 2a.

In a regional scope, the Black River is in a meandering area with a low surface slope, causing low soil erosion that carries sediments to the Tragadero River (Figure 1). MAGR has flood risk, and heavy rains (1500 mm a year [59]) can cause hydric emergency [60], altering the water properties in the treatment plant inlet and increasing the Parana River flow rate.

The highest water flow in Parana River occurs between February and March, with values over 30,000 m³/s, with a mean of 17,000 m³/s [61]. The change in water flow, land use and hydroelectric development (by dam constructions) alters the hydrologic characteristics of Parana waters, thus affecting water treatment operations in the last decades [61].

To estimate turbidity from R_S it was necessary to inspect the relationship between these quantities. Figure 3 illustrates the spectral signatures of all the observations in light grey lines. Grouping the data observations by turbidity ranges, the obtained mean spectral signatures are shown in black lines. Lower turbidity (<150 NTU) presented the lowest R_S . As the turbidity range increased, the spectral signature response rose until values higher than 1050 NTU.

Figure 3.

Mean spectral signatures per turbidity range for S2-MSI bands

Bands B01, B11 and B12 (Figure 3) are not sensitive to turbidity change since the points for different turbidity ranges remained in the same position. Bands B05, B06 and B07 presented the highest changes. These bands were related to algorithms for turbidity estimation [5].

Several models were tested to estimate the inlet water turbidity for the treatment plant, with predictor variables selected according to the model. For RF, S2-MSI bands shown in Table 1 were used as predictors. Traditional linear models took account of the following variables: an interaction between B06 and B07; individual bands B05, B06 and B08; and the spectral index NDTI. The aforementioned spectral bands were selected according to the results from Figure 3.

The characteristics and performance metrics for all proposed models are resumed in Table 3. Following model selection based on these metrics, the model is finalised using the preserved testing dataset, with observations not used in the training. The best results were achieved by the RF model. R² values for individual bands B05, B06, and B08 were 0.693, 0.732 and 0.736, respectively. In a similar work [5], also in turbid lakes, R² values for the same bands were 0.83, 0.66 and 0.63, respectively.

RF model was selected in the following analyses since it presented the best performance metrics, with the lowest deviations (RMSE) and highest correlation (R²). This machine learning algorithm is superior to the proposed usual regression models in capturing the complex relationships between satellite spectral data and turbidity values. An RF model with proper modifications performed better than other machine-learning options for turbidity estimation [62]. The interaction model between B06 and B07 was the second-best model combining bands in the red edge related to sediments in water [63]. In comparison, the linear model using a single band performed poorly. Unlike other studies, the NDTI index presented the lowest R² and the highest RMSE [64].

Table 4 shows the tuned hyperparameter values and the main characteristics of the final RF model.

After the model selection, the last fitting was performed. For the RF model, the final performance metrics were obtained by the testing dataset. These observations were kept apart, so they have no influence on the modelling. The metrics obtained with 39 data points were RMSE = 143.9 NTU and R² = 0.913. These values are different from Table 3 data derived from the training data set and used only for model selection.

The comparison between measured and estimated observations in the testing split is shown in Figure 4. The solid line represents the linear relationship between estimated and measured turbidity, with a dashed line at 45°. Lower estimated turbidity values are closer to the real values. The deviations increase for higher turbidity, with estimates being lower than measured values, as indicated by the solid line below the dashed line. The outcome variable presented a wide range, with many observations under 100 NTU and some measurements as high as 1100 NTU.

Figure 4.

Estimated and measured turbidity using the validation dataset

The measured and estimated turbidity values for the validation dataset (Figure 4) are shown as a time series plot in Figure 5. The crosses represent the estimations made by the RF model, while the turbidity measurements are plotted as a solid line. The number of samples in the testing dataset is shown in the top right corner.

Figure 5.

Time series of estimated and measured turbidity in the validation dataset

The most significant differences between estimated and measured turbidity in Figure 5 are within the larger values, equivalent to Figure 4. The estimations followed the same trend seen in the time series in Figure 2, with high turbidity at the beginning of the year and lower later.

The complexity of the RF model is difficult to explain since the explicit form is not as clean as a more straightforward linear model. The explanatory technique of global feature importance can assist in understanding the driving predictors of RF aggregated in all training observations.

The results of the global feature importance technique for the obtained RF model, analysing each spectral band, are shown in Figure 6. The technique employs the notion of the overall change in the model due to the perturbation of a specific variable [65]. A permutation-based approach is a valuable tool for model explanation since after the permutation of said variable, the model performance is expected to decrease [49].

Figure 6.

Global feature importance of S2-MSI spectral bands in RF model

Spectral band B07 presents the most effect in the model, according to Figure 6, since the boxplot had the highest RMSE (104.9 NTU). Close to B07 were B06, B08, and B05. Spectral bands B05 [66] and B08 [13] have been reported to be related to turbidity. The lowest effects were given by B01, B02, and B03 since the perturbation of these bands had a much lesser impact on the overall model. The vertical dashed line represents the base RMSE.

The most influential bands ranged from 704 nm to 830 nm, with the least influential between 440 nm and 500 nm. For comparison, the same method in the research on the North Tyrrhenian Sea [67] gave a similar result, but the band importance order was B05 (the highest), followed by B07 and B06. Turbidity estimations were most successful at the wavelength between 700 nm and 800 nm for surface water [68]; this range includes B05, B06 and B07.

Partial dependencies profiles allowed to show the change in the expected value of a model estimate alongside a single explanatory variable [65]. According to the global feature importance technique, B07 was the spectral band with the highest effect on the RF model. Figure 7a shows the partial dependencies profile obtained for this band. The thin grey lines in Figure 7 correspond to 100 randomly selected observations from the training dataset. The black line indicates the mean. The effect of B07 (Figure 7a) on turbidity estimates was constant until R_s = 0.12, then started to increase until its highest effect at R_s = 0.2. In this range of surface reflectance, the turbidity changed from 218.7 NTU to 353.8 NTU. For comparison, Figure 7b corresponds to an identical analysis for B01, the band with the lowest feature importance (see Figure 6). The partial dependency profile of this band was constant, meaning that the turbidity presented no change in the entire range of R_S from B01 values. This result was consistent with Figure 3 and Figure 6.

Figure 7.

Partial dependencies profiles for spectral bands B07 and B01

The obtained RF model was applied to the spectral values from the Barranqueras River to evaluate the spatial turbidity distribution. Figure 8 shows the maps for four different dates from 2020. The yellow triangle on the top centre of each panel represents the water plant location. A water mask was applied to the region of interest to extract only pixel values from the Barranqueras River. Figure 8a (2020–01–07) and Figure 8b (2020–12–22) indicate relatively low turbidity with a wide dispersion of values. Figure 8b (2020–04–11) and Figure 8c (2020–08–24) present a narrower turbidity dispersion, thus the colour homogeneity. For the former date, the estimated turbidity values were high; for the latter, turbidity values were lower.

Figure 8.

Barranqueras River turbidity maps for four different dates

To better understand the turbidity spatial distribution, histograms were plotted to showcase the estimations of dispersion alongside the Barranqueras River, Figure 9. The bin width was set to 10 units. Figure 9a (2020–01–07) presented relatively low values with a wide dispersion; the median turbidity was 117 NTU. High values and a narrow dispersion, with an 889 NTU median, were observed in Figure 9b (2020–04–11). The lowest turbidity distribution was obtained in Figure 9c (2020–08–24), presenting a single peak at 69 NTU. Finally, the values increased in Figure 9d (2020–12–22) until a 128 NTU median and a wider dispersion.

At extreme values, the turbidity dispersion was low, as seen in Figure 9b and Figure 9c. The estimations observed in Figure 8 maps and Figure 9 histograms followed the same measured turbidity trends as Figure 2a for 2020.

Figure 9.

Turbidity histograms per the exact dates as seen in Figure 8