Chlorophyll-a Estimation in 149 Tropical Semi-Arid Reservoirs Using Remote Sensing and Machine Learning

1. Introduction

Eutrophication, an essential manifestation of water pollution in freshwater reservoirs, results from the excessive nutrient loads that cause algal blooms [1,2]. The excessive algae and aquatic plant growth result in oxygen depletion, ultimately causing severe impacts on aquatic ecosystems and considerably increasing the costs related to water treatment [3,4,5,6].

Chlorophyll-a (Chla), a photosynthetic pigment in major algae groups, is widely used as a critical indicator of phytoplankton presence [7,8]. Because the abundance of algae can reflect the state of eutrophication, Chla is one of the most important parameters for evaluating the trophic condition of water bodies [9,10]. While Chla concentration has long served as a critical parameter in monitoring harmful blooms, accurately predicting Chla in reservoirs has proven to be a persistent challenge [11,12]. This difficulty is mainly due to the non-linear and non-stationary characteristics of Chla concentration, which are influenced by anthropogenic and hydrometeorological factors [2].

In regions characterized by semi-arid climates, such as Northeastern Brazil, establishing an extensive network of multi-purpose artificial reservoirs has emerged as a reliable solution to confront the water scarcity challenges imposed by environmental constraints [13,14]. These reservoirs are notably susceptible to eutrophication due to hydroclimatic characteristics that combined favor photosynthesis and biodegradation [15,16], such as interannual variability of precipitation and stored volume [17], high temperatures and evaporation rates [18], and prolonged hydraulic retention time [19]. Moreover, this susceptibility is further aggravated by continued anthropogenic pressure on water bodies due to internal enrichment from aquaculture practice [3,20], inadequate coverage of sanitation systems [21], and a dense reservoir network [22,23].

Therefore, regularly monitoring Chla concentrations is crucial for implementing effective water quality management strategies to prevent further deterioration [24]. However, traditional sampling methods are expensive, time-consuming, and impractical for many reservoirs [25,26]. As an attractive option, satellite remote sensing and machine learning (ML) techniques offer a cost-effective approach for monitoring Chla concentrations and their spatiotemporal variations, providing large-scale data on complex environmental systems [27]. The Sentinel-2 constellation has been proved to be a valuable asset in monitoring inland and coastal waters once it has improved spatial resolution compared to other freely available sensing systems data, such as Landsat 8 [28,29].

The ML approach retrieves complex non-linear relationships within satellite data, capturing the underlying structure bonding the satellite data and the desired target variable [30,31]. Combining ML architectures with remote sensing data has been able to provide top-notch results in a plethora of scientific fields, such as solar irradiance forecasting [32], mapping of mineral extraction sites [33], forest fire mapping [34], and crop water stress evaluation [35].

ML and satellite data performance has also been explored for Chla monitoring inland and ocean waters. For ocean water, including coastal waters, specific Chla ML forecasting models, namely OLCI Neural Network Swarm (ONNS) and Ocean Colour 4 for MERIS (OC4ME), using Sentinel-3 satellite data, showed good performance for such a task [36,37]. However, other possible ML paradigms can be found for inland and coastal waters. A random forest (RF) based model was developed in [38]. In their study, the authors used Sentinel-2 imagery to retrieve Chla concentrations for lake Chagan in China. Their proposed model provided good performance in determining the Chla concentrations and complyingwith the biological mechanism in lakes, offering robust results to seasonal changes.

In Cao et al. [9], the authors used Landsat-8 remote sensing data together with extreme gradient boosting tree model (XGBoost) ML to determine Chla in lakes located in China. Their approach was implemented to analyze the spatiotemporal data from 2013 to 2018, and it demonstrated satisfactory performance in identifying the Chla behavior in the study location. In Hu et al. [39], the authors developed methodologies to mitigate spectral noise in satellite data to improve the performance of ML models to estimate Chla in global oceans using remote sensing from several satellites. Their results proved that the support vector regression (SVR) was the best-performing ML approach, surpassing the traditional band-ratio models and providing reduced image noise.

The Group Method of Data Handling (GMDH) has also been applied to hydrological scenarios, including Chla estimation [40], water quality prediction [41], and image classification for plant diseases [42]. However, there is a gap in the knowledge regarding the usage of this approachin modeling Chla concentration using satellite data, which the present study aims to fulfill. Additionally, the present study seeks to provide a deep insight into the performance of ML paradigms when applied to a vast area containing heterogeneous reservoirs.

Given this scenario, the potential to increase algae blooms and further degradation of these aquatic systems has increased the need to study often poorly monitored reservoirs in semi-arid regions. Therefore, this study aims to estimate Chla concentration through combined remote sensing and machine learning techniques. For this, the following specific objectives were pursued in this research:

• A comprehensive investigation on several input parameters for Chla modeling, including all the 13 bands of the MSI onboard of Sentinel-2 constellation and 16 different spectral indexes.
• A comprehensive analysis and characterization of all the 149 tropical reservoirs extensively spread across the state of Ceará, located in the Brazilian semi-arid region.
• The usage of stepwise approach on parameters selection.
• The investigation of different machine learning paradigms for modeling Chla values in heterogeneous reservoirs distributed over a vast region.
• To fulfill the gap in the knowledge regarding the usage of the GMDH ML model for Chla modeling using remote sensed data and spectral indexes.

2. Materials and Methods

2.1. Study site location

The Brazilian semi-arid region is one of the most populated semi-arid regions in the world, with approximately 28 million inhabitants and an area that occupies 12% of the national territory [43]. The state of Ceará, inhabited by around nine million people, covers an area comparable in size to England of 150,000 km², where 98.6% of its territory is in the semi-arid region [44]. The region is characterized by extreme climatic events such as recurrent droughts, sporadic flooding and high inter-annual variability [45,46]. According to the Köppen classification [47], the climate is 'BSh', characterized by a mean precipitation of 750 mm per year, a potential evaporation rate of 2000 mm per year, a mean annual temperature of 31°C and negative water balances for most of the year [13,48]. Moreover, Ceará's state is influenced by two distinct seasons: a rainy season from January to April, when 80% of the total precipitation occurs, and a dry season during the rest of the year [49].

Since Ceará's water supply relies predominantly on artificial reservoirs, the state has a dense reservoir network, serving as the primary source for over 90% of the region’s water resources [20] and with a storage capacity of approximately 18.6 billion cubic meters. Notably, the three largest reservoirs in the state, Castanhão (6,700 hm³), Orós (1,940 hm³), and Banabuiú (1,600 hm³) collectively represent approximately 55% of the total storage capacity [50].

The study area includes data from 149 monitored reservoirs distributed across Ceará in twelve watersheds. The longitude, latitude, and basic information of the reservoirs are listed in Table S1 in the supplementary material. These reservoirs are mainly used for human water supply, aquaculture, fish farming, and irrigation. As a result, pollution is predominantly attributed to nonpoint sources such as livestock, agriculture, and soil erosion, and point sources originating from sewage and fish farming [51]. The trophic condition of these reservoirs fluctuates with seasonal changes and flooding or drought events [52]. Figure 1 illustrates the geographical location of Ceará within the Brazilian territory and in the semi-arid region, as well as the location of the reservoirs distributed across its area.

2.2. Water quality data

This study used data from 149 spatially distributed reservoirs (Figure 1), covering the years 2015 to 2021, including 1,399 Chla samples. This information was collected from the database provided by the Portal Hidrológico do Ceará platform, a system in which monitoring has been carried out by the Water Resources Management Company (COGERH) [50], a public agency maintained by the local government, presenting consistent time series of hydrological and water quality parameters. The company carries out monitoring campaigns in all seasons of the year through quarterly sampling campaigns and in situ measurements of parameters.

The Chla concentration data provided by the COGERH was obtained through in-situ water samples. The samples were gathered at each reservoir sampling point at 0.3 meters from the surface, collected into a dark flask to avoid exposure to light. Later, the sample’s pigments were extracted using a solution of 90% acetone to cold extract them. Finally, the pigments were used to assess the Chla concentration for each sample were analyzed by accredited laboratories according to a standardized protocol APHA 10200 H spectrometric method [13,53].

2.3. Sentinel-2 satellite data

2.3.1. Band data

The Sentinel-2 mission is an effort of the European Space Agency (ESA) to monitor the Earth’s environment. This mission comprises a constellation of two polar-orbiting satellites, namely Sentinel-2A and Sentinel-2B. Both satellites are placed in the same sun-synchronous orbit, aiming to monitor the Earth’s environmental changes [54,55]. The mission started with the launch of the first Sentinel-2A satellite in June 2015. The latter deployment of the second satellite, Sentinel-2B, in March 2017 reduced the revisiting time from 10 days to 5 days [56].

Besides the brief revisiting time, the Sentinel-2 mission is equipped with Multispectral Instruments (MSI), which provide information for different spatial resolutions, ranging from 10 m to 60 m, given the state-of-the-art anastigmatic telescope within it [55,57,58]. The MSI can register data from 13 spectral bands, varying from visible to near-infrared (NIR) and short-wave infrared (SWIR), providing high-resolution data for both inland and coastal areas [29,55,59].

The remote-sensed data used in the present work undergoes Level-1C processing. For this level of data processing the Earth’s spatial region is partitioned into tiles based on the UTM/WGS84 projection, with each tile separated by a distance of 100 km. Radiometric values for each tile’s pixel are determined for the top-of-atmosphere reflectances, which are later converted to radiance [60,61]. The quality of Level-1C processing is ensured by the ESA monthly, following rigorous quality standards [62].

The combination of reduced revisiting time, high spatial resolution, and wide range of spectral bands deem the Sentinel-2 an important mission with agriculture applications and forest monitoring [28]. Table 1 compiles the information for each spectral band.

Table 1 informs the wavelength for each band captured by the MSI. It is possible to notice that the wavelengths around 700 nm (NIR) suggest that the Sentinel-2 constellation is suited for capturing the phytoplankton spectral characteristics, including Chla, as the microscopic organisms cause a surge in spectral reflectance around the 700 nm mark [13,29,64].

Furthermore, given the satellite resolution, the water samples from the reservoirs were taken at a minimum distance of 60 meters from the shoreline. This precaution ensures that the pixel, corresponding to the smallest spatial resolution is entirely encompassed within the reservoir boundaries. Also, given the satellite revisiting time, most of the remotely sensed data had an interval of 1 to 2 days. Given that the Chla concentration does not vary significantly for the given period, this time-range is valid for our study application.

2.3.2. Satellite spectral indexes

Satellite spectral indexes are derived from mathematical equations combining two or more spectra of the satellite bands. These indices are a helpful approach to extracting information from the spectral bands pixelwise to model terrestrial processes and features, like vegetation, water, urban development, and agriculture [63,65,66]. A comprehensive investigation of 16 different indices and their impact on the model’s result was performed in the present study. Their mathematical formulations are displayed in equations 1 to 16.

$$NDVI = \frac{Band 8 - Band 4}{Band 8 + Band 4}$$

(1)

In Equation 1, we have the normalized difference vegetation index (NDVI) formulation. This index is vastly used in remote sensing, being primarily used for the evaluation of green areas and the changes related to them. It is a valuable input on different remote sensing applications [67,68].

$$GNDVI =\frac{Band 8 - Band 3}{Band 8 + Band 3}$$

(2)

Equation 2 shows the formulation for the Green Normalized Difference Vegetation Index (GNDVI), an adapted version of the NDVI index aimed explicitly at detecting Chla in the vegetation [69,70].

$$EVI = 2.5\cdot \frac{\left(Band 8 - Band 4\right)}{\left(Band 8 + 6\cdot Band 4 -7.5\cdot Band 2 + 1\right)}$$

(3)

The enhanced vegetation index (EVI) is similar to the previous NDVI but removes the impacts of the atmosphere and the soil over the vegetation signal [71,72].

$$SAVI = \frac{Band 8-Band 4 }{Band 8 + Band 4 + 0.428}\cdot 1.428$$

(4)

The soil-adjusted vegetation index (SAVI) improves the NDVI index by considering the soil effects due to multiple scattering of soil [72,73].

$$NDMI=\frac{Band 8 - Band 11}{Band 8 + Band 11}$$

(5)

The normalized difference moisture index (NDMI) is used to verify changes in vegetation physiology by determining its water content [74,75].

$$MSI = \frac{Band 11}{Band 8}$$

(6)

The moisture stress index (MSI) is used to evaluate changes in the water content in the vegetation via canopy stress analysis. It is also used to indicate water concentration in the soil [76,77].

$$GCI = \frac{Band 9}{Band 3}-1$$

(7)

The green chlorophyll vegetation Index (GCI), as its name implies, is applied to the remote sensed data to estimate chlorophyll concentration in the vegetation and, consequently, determine the health of the analyzed vegetation [78,79].

$$NBRI =\frac{ Band 8 - Band 12}{Band 8 + Band 12}$$

(8)

The normalized burned ratio index (NBRI) seeks to identify the occurrence and severity of natural or man-made fires in vegetation areas [80,81].

$$BSI =\frac{ \left(Band 11 + Band 4\right) - \left(Band 8 +Band 2\right)}{\left(Band 11 + Band 4\right) + \left(Band 8 +Band 2\right)}$$

(9)

The bare soil index (BSI) formulated in Equation 9 retrieves information from the vegetation in cases where its coverage is less than half of the assessed area. This index allows us to determine the vegetation health concerning the exposed soil area [82,83].

$$NDWI = \frac{Band 3 - Band 8}{Band 3 + Band 8}$$

(10)

The normalized difference water index (NDWI) retrieves information from water bodies effectively from remote sensing data [84,85].

$$NDSI = \frac{Band 3 - Band 11}{Band 3 + Band 11}$$

(11)

The Normalized Difference Snow Index (NDSI) is a tool that detects snow cover in a specific area by analyzing the light reflection properties of ice. This index retrieves information by distinguishing snow coverage from other surfaces and adjusting for atmospheric and terrain effects [86,87,88].

$$NDGI = \frac{Band 3 - Band 4}{Band 3 + Band 4}$$

(12)

The normalized difference glacier index (NDGI), similar to NDSI, identifies glacier coverage in a region mainly composed of snow, ice, and debris [89,90].

$$ARVI = \frac{Band 8 - 2\cdot Band 4+Band 2}{Band 8 + 2\cdot Band 4}$$

(13)

The atmospherically resistant vegetation index (ARVI) is an improvement over the NDVI index by implementing atmospheric corrections. The ARVI is especially useful for regions under dense aerosol coverage [91,92].

$$SIPI = \frac{Band 8 - Band 2}{Band 8 - Band 4}$$

(14)

The structure-insensitive pigment index (SIPI) was initially proposed to identify vegetation stress via the ratio between carotenoid and chlorophyll in vegetation. It is also useful for analyzing vegetation structures with different canopy configurations [93,94].

$$SWM = \frac{Band 2 + Band 3}{Band 8 + Band 11}$$

(15)

The sentinel water mask (SWM) specifically seeks to detect water data from the Sentinel-2 constellation [95].

$$AWEI = 4\cdot \left(Band 3 - Band 11\right)-\left(0.25\cdot Band 8 + 2.75\cdot Band 12\right)$$

(16)

The automated water extraction index (AWEI) is used to accurately detect water given various environmental interferences [96,97].

A first analysis of the indexes may result in the conclusion that some of them, e.g., NDSI and NDGI, may not be suitable for the study location due to the semi-arid characterization of the Ceara state. However, a careful examination of Equations 11 and 12, together with Table 1, elucidates that the mentioned indexes result from bands within the spectral range that favour Chla identification and are of high spatial resolution [13,29,64]. Therefore, these indexes may still carry relevant spatiotemporal information that may help uncover the latent relationship between the input parameters and Chla, thus improving the performance of the ML models. This potential improvement provides a compelling reason for further investigation into the participation of these indices in Chla modeling. Figure 2 and Figure 3 show the correlation between the satellite bands and Chla and the indexes and Chla, respectively. In both images, the lighter the color, the more correlated the variables are.

Figure 2 shows that the bands share a significant correlation between themselves, but band 10. In contrast, in Figure 3, the indexes show less correlation with each other, indicating potential to be used by ML models due to their low collinearity. Regarding the correlation with the Chla attribute, the bands present significantly lower values, up to ten times smaller than the correlation between the indexes and Chla.

The elevated correlation within the bands’ data does not necessarily mean a positive impact on the performance of the forecasting model. Highly correlated variables, when used as inputs on a predictive model, may add noise to the dataset, increasing the model’s variance and thus reducing its performance [96,97]. However, discarding highly and low correlated variables altogether may also harm the model’s performance since they may still carry relevant spatiotemporal information relating to the inputs' and outputs' attributes, leading to improved forecasting performance [30].

2.4. Machine learning models

Chla was estimated using satellite data and machine learning models. The first implemented approach was the random forest (RF). The RF is a model composed of trees trained using the bagging method of resampling considering just a subset of predictors (Figure 2 and Figure 3) [100,101], deeming it an ensemble model. The trained trees are low correlated, reducing the ensemble model variance and improving performance [98,99,100]. Different published works have investigated the RF methodology [101,102,103].

XGBoost [104] was also implemented as a standalone model and used as benchmark for the models' performance. This tree-based approach is an ensemble model and an extreme improvement over the random forest. It consists of bag sampling smaller tree models, which are combined into a larger and more robust tree-model, reducing the model variance while improving its generalization and reducing the model's tendency to overfit [98,105]. The XGBoost model can handle missing data and manipulate increasing dataset size, keeping its generalization. This approach reached excellent results when applied to different time-series forecasting tasks [106,107,108].

Another ML used in Chla prediction in this work was the k-nearest neighbors (k-NN). The k-NN algorithm is a supervised ML model, which uses non-parametric vectors to determine an unknow point [109]. However, given that it is also based on the distance between points distributed on a possible multi-dimensional space, a distance metric must be implemented, often the Euclidean distance [109,110,111]. The k-NN architecture can be applied to data classification and regression model cases [109,112]. For the former application, the model classifies an unknow datapoint considering its neighbors by a simple voting system, assigning the unknow point the label of the most common class around it. The latter approach has a similar implementation, but this time, a continuous value is assigned to an unknow point given the average value of the target attribute of its neighbors [109,112,113]. Despite its simplicity, the regression performed by the k-NN approach offers competitive results within the ML field and has been explored for different scenarios in previous studies [114,115,116].

The Support Vector Machine (SVM) was also implemented in this study. The SVM is a flexible ML approach with diverse applications for classification, regression, and outlier detection [117]. Initially proposed by Cortes and Vapnik in 1995 [118], with previous formulations dating back to the 1960s, the SVM is a generalization of the maximal margin classifier [112]. The SVM learns a boundary function discriminating the dataset when implemented for classification. At the same time, regression provides a best-fitting function describing the data behavior by taking into consideration its extreme attribute vectors. A unique feature of SVM is the use of kernel functions that allow the dataset to be manipulated into higher-dimensional spaces, making it possible for the model to learn complex non-linear relationships by applying convex optimization without being computationally expensive, reducing the training error [109,119,120]. However, one drawback of the SVM approach is that it does not handle large datasets efficiently, requiring extended computational time to be trained [121,122].

The least absolute shrinkage and selection operator (LASSO) regression [123] was another implemented ML methodology. This approach is a more straightforward ML methodology, which seeks to implement the best linear regression line in the dataset. Besides that, the LASSO paradigm is also a regularization and parameter selection approach, making its results more interpretable than other traditional ML models [98,123].

Lastly, using satellite data, we investigated the GMDH ML model for the Chla modeling. It was first proposed by Ivakhnenko [124] as an alternative to address the challenges of linear dependency and equation complexity for higher dimensional problems and small data sequences [40,125]. This methodology is a feedforward unidirectional ML model, similar to a multilayer perceptron [126,127]. It is also a self-organizing model, indicating that its parameters are selected automatically, not needing parameter tuning [128]. The resulting value by the GMDH model is a quadratic approximation, using pair combinations of the input variables [129,130], which are used to model the relationship between input and the output parameters [131]. Unlike other artificial neural network paradigms, the GMDH does not need vast training data, as the estimation of its parameters is automatically determined without recursion [40]. The performance of the GMDH model was verified in previous studies for time-series challenges, including hydrological applications [132,133,134].

2.5. Evaluation metrics

To assess the forecasts by the proposed models, we opted to calculate the metrics root mean squared error (RMSE), normalized RMSE (nRMSE), mean absolute error (MAE), mean absolute percentage error (MAPE), mean bias error (MBE), and coefficient of determination (R²). Their equations can be verified in [135] for R², and [136] for the remaining metrics.

2.6. Dataset preprocessing and attribute selection.

The dataset standardization seeks to restrain the dataset to within the same dimensionless scale, setting the mean value equal to zero and the standard deviation to one, often improving the model’s performance [32]. Besides the data standardization, in this study, the Yeo-Johnson transformation [137] was implemented for some of the ML models. This transformation is an improvement over the Box-Cox approach, which is restricted to handling positive numbers only. The Yeo-Johnson transformation is based on the power transformation with different parameters to positive and negative values and is presented in Equation 18 [137].

$$\psi \left(y,\mathsf{\lambda }\right)=\left\{\begin{array}{ll}\frac{(y+1{)}^{\lambda }+1}{\lambda }& y\ge 0, \lambda \ne 0\\ \log \left(y+1\right)& y\ge 0, \lambda =0 \\ -\frac{(-y+1{)}^{2-\lambda }-1}{2-\lambda }& y<0, \lambda \ne 2\\ -\log \left(-y+1\right)& y<0, \lambda =2\end{array}\right.$$

(18)

Where the transformed value ψ is a function of the original attribute value, $y$, and a parameter $\lambda$, which is determined via maximum likelihood [101,102]. This transformation seeks to reduce the data skewness, approximating the original dataset distribution to a normal distribution as $\psi \left(y,\lambda \right) \sim N\left(\mu , {\sigma }^2\right)$ [138,139].

The best parameters were selected among the bands and indexes. To this end, using a step-by-step approach [112], we started investigating the influence of the bands only. First, we investigated the performance of just one band at a time, from band 1 to 12, and we selected the best one that returned the best R², using the XGBoost model as the benchmarking approach. After that, we tested the combination of the already selected bands, with the remaining ones, again, one by one, and kept the band that returned the highest coefficient of determination value (forward variable selection). This process was repeated until all the bands were investigated and their best combination was kept [112].

The same procedure was repeated for the selection of the indexes. This time, we started the investigation with the best band selection previously set. Then, we added one index at a time to the input parameters, selecting the index which yielded the highest R². Again, the process was repeated until all the indexes were assessed. The selection of the best bands and indexes are presented in Figure 4 and Figure 5, respectively.

As depicted in Figure 4, the modeling for Chla substantially improved with the combination of bands 6, 4, 10, 7, 8a, 9, 5, 12, and 3, resulting in an R² of 0.36. Figure 5 shows that the addition of the indexes also managed to increase the coefficient of determination to 0.50 after including NBRI, NDSI, AWEI, SAVI, NDVI, MSI, GCI, GNDVI, and ARVI in the selected input data. Interestingly, the NDSI, an index related to snow coverage, promoted improvements over the Chla identification, as indicated in Figure 5, proving that it conveyed spatiotemporal information for the modeling approach. Therefore, the results discussed in the subsequent section are based on grouping the aforementioned spectral bands and indexes.

3. Results

3.1. Limnological behavior

The mean value for Chlorophyll varied widely among reservoirs and according to seasonality. Considering the complete dataset, Chla values ranged from a minimum of 1 μg/L to a maximum of 1001.78 μg/L, with a mean value of 39.62 μg/L and a standard deviation of 65.78 μg/L. The boxplot in Figure 6-A shows the inter-annual and intra-annual distribution of the dataset. Mean values decreased throughout the years (from 81.34 μg/L in 2015 to 27.84 μg/L in 2021), while maximum values did not reflect this tendency, with recorded values increasing, for example, from 707.06 in 2015 to 1001.78 in 2016, indicating the occurrence of concentration blooms. Regarding the intra-annual distribution of Chla, higher concentrations were detected in the rainy seasons for most of the studied years, except in 2017.

The state’s three largest reservoirs were also analyzed (Figure 7-B). The mean recorded values for Orós were the highest (31.89 μg/L), between 9.20 and 60.94 μg/L, while Chla for Banabuiú was found in lower mean concentrations (21.64 μg/L) between 1.82 and 54.55 μg/L. For Castanhão, values ranged from 5.65 and 78.72 μg/L with mean Chla value of 29.03 μg/L.

3.2. Results for Chla estimation by the ML models

The proposed models were built using the variables presented in Figure 4 and Figure 5 as input parameters. We used data from the Sentinel-2 satellite and on-site measurements from the 149 reservoirs to train and test these models. We adopted a 70/30 split for the training and validation datasets. These datasets were randomly assembled using data spanning from 2015 to 2021. It is important to notice that given the random separation for the training and validation data, not all the reservoirs may have been included in the validation dataset. However, this does not mean that the validation dataset is not significant and statistical characteristics similar to the training dataset.

The hyperparameters of the used ML models, except for the GMDH self-organizing approach, were selected using the GridSearch tool for Scikit-Learn, with 5-fold cross-validation [140,141]. The outcomes for Chla modeling by each one of the assessed models are presented in Table 2.

Table 2 shows that the ML approaches reached similar results regarding the RMSE metric, but for the GMDH, which surpassed all of the assessed paradigms, achieving RMSE of 20.38 μg/L. The LASSO model achieved the highest error value for RMSE, with 89.87 μg/L, followed by the k-NN model, with RMSE value of 61.82 μg/L. The tree-based models XGBoost and RF achieved similar RMSE values of 55.60 μg/L and 56.75 μg/L, respectively. The SVR output scored an RMSE value equal to 50.64 μg/L. It is essential to disclose that an analysis of the RMSE alone may be misleading to assess the performance of the models. The coefficient of determination is another crucial factor to consider in this scenario. This parameter indicates the total variance of the dependent variable Chla, which can be adequately forecasted by the input parameters and may be understood as an indication of the accuracy achieved by the model [142]. This behavior can be better visualized in Figure 7.

Figure 7 shows how the predicted data compares with the actual measured Chla values. We used log values to facilitate the comparison due to differences in the variables’ scales. The histograms showcase the normal distribution of the data after the logarithm transformation. In Figure 7, given the log scale, minor fluctuations in Chla values correspond to negligible alterations in the satellite-captured data. This is not related to the used GMDH model itself but is a consequence of the chosen scale. In fact, the remote sensed data lacks the sensibility to convey sufficient information for low Chla values, since the main characteristics of high Chla concentrations are their tendency to converge to a specific spectral pattern of green, while for low values they tend to converge to different spectral patterns that depend on the characteristics of each reservoir.

Still in Figure 7, it is possible to observe that both predicted and real data present similar distributions, which indicates the GMDH model's good performance in predicting Chla levels. The regression line displayed in the plotting area shows a positive correlation between the predicted and real Chla data, further attesting to the robust GMDH performance. The points clustered around the regression line also indicate the superior performance of the GMDH algorithm, especially for extreme values as depicted by the top-rightmost points in the graph area. Comparing the results in Table 2, we can observe that the highest R² value reached by the GMDH model was 0.91, representing a significant improvement of 57% over the SVR result, the second best-performing approach regarding the same parameter. The third- and fourth-best performing models were the tree-based ones XGBoost and RF, which achieved R² equal to 0.50 and 0.48, respectively. The k-NN performance, R² equal to 0.38, and the LASSO model, R^2, equal to 0.38, achieving a coefficient of determination equal to 0.41, was the worst-performing modeling paradigms. Analysis on the values of MAE and MAPE indicates that the GMDH is able to output more accurate and precise results when compared to the others assessed models. The negative values for MBE is a trend presented in all the investigated paradigms. It indicates that the ML models share the tendency to underestimate the outputs for Chla value.

4. Discussion of the results

Mainly in tropical regions, the mechanisms that regulate phytoplankton growth require an advanced analysis due to the complexity and nonlinearity of the relationship between chlorophyll and physicochemical/environmental factors [13,143]. The approach presented in this study achieves this analysis using a highly heterogeneous collection of in-situ observations, investigating the performance of different ML models to determine Chla levels in 149 dams in the Brazilian state of Ceará. Although most dams were predominantly eutrophic throughout the years, different human activities and pollution sources have contributed to eutrophication processes lead to Chla spatiotemporal fluctuations and algal blooms [17,20,144,145].

We applied the forward stepwise approach to select the parameters, including bands and spectral indexes. To our knowledge, this method has not yet been applied in previous remote sensing studies. The forward stepwise approach s, significantly improving the Chla determination for the ML models. Existing literature regarding the influence of spectral bands on Chla determination can be found. In work [146], the authors applied the SHAP analysis to explain the influence of the Sentinel-2 spectral bands for estimating Chla values. Their results elucidated that band 3, band 2, and band 8 were the top three most influential parameters for Chla determination. Bands 3 and 2 exhibited a positive correlation with the Chla, while band 8 showed a relevant negative correlation with the same parameter. A similar approach was conducted by Kim et al. [147]. In their work, the SHAP analysis showed relevant participation in Chla prediction of red bands, i.e., bands 4, 5, and 6, and blue and green bands (Table 1), a similar conclusion found by Ha et al. [148].

The influence of different spectral indexes for Chla modeling was investigated in previous works. Castro et al. [149] showed that indexes merging red and NIR bands yielded the best outcomes for determining Chla concentrations in small reservoirs. Similar conclusions were drawn by Buma and Lee [150] and by Aubriot et al. [151], which attested to the importance of bands within the red spectral range for Chla characterization in a lake in Chad, and for the Rio de La Plata, respectively. On the other hand, Viso-Vazquez et al. [152], showed that the green band, i.e., band 3, promoted the highest correlation between the remote sensing data and Chla levels.

Our results show the GMDH methodology could efficiently identify the latent non-linear ties ruling input and output attributes. Moreover, the GMDH approach proved to be a more robust model for analyzing satellite imagery, being more resilient to noises and satisfactory results when handling different band information, which also can be attributed to its improved performance.

To better understand where the GMDH results lie within the literature, we compared our results with those found in previously published works. Such evaluation, however, may not be representative given the different methodologies used in different studies, models, input attributes, and the study area. Yet, comparing their evaluation metrics is still a viable approach to assess the performance of different models [153]. Table 4 compiles the results for the GMDH model, while Table 5 gathers the results found in the literature.

Table 5 shows that the models based on the deep learning methodology all performed remarkably well in determining Chla levels. Compared to the results found in our study, it is possible to notice that the R-squared values for references [154,155] are in the same range, over 90%. Nevertheless, regarding the work by Guo et al. [154], the combined utilization of Sentinel-2 and Landsat-8 data significantly enhances the machine learning model's performance. This improvement is primarily due to the reduction in the revisit time, which subsequently minimizes the variance in the dataset. Consequently, this leads to a more robust and accurate machine-learning model. It is important to note that, although their results are slightly superior, the study location is limited to only one place.

Adding more water bodies is expected to add variance to the dataset. The model in reference [156] was built using modeled data, including 11 times more water bodies than the previous studies. Their results proved to be slightly inferior to the other two previous works and more relatable to the ones found in this present assessment, with RMSE in the same magnitude order.

Regarding the R², our GMDH model showed a nearly 10% superior value. However, it must be stated that since the data used for training the ML model was simulated, it may not consider several naturally occurring situations. This would lead to a more homogeneous dataset with less variance, improving the ML model performance compared to the model proposed in our study. Another significant difference is the time window used for testing the developed model in reference [156], which is significantly smaller than the one used in our model [160], which reduces the dataset variance, improving the ML performance.

The DL approaches presented in Table 5 still require a considerable amount of data to yield reliable outputs. Depending on the data availability for the study location, this characteristic may pose a major drawback. On the other hand, the GMDH paradigm can be promptly implemented with less information available and does not require extensive training datasets [40], deeming them more flexible for different situations.

References [157,158] applied SVR to determine Chla. It is noticeable that although the models are the same, their methodology was disparate. Regarding reference [157], the RMSE values were within the same order of magnitude. Besides that, the R² values were relatively close to each other. Contrarywise, reference [158] performed a much broader study regarding several lakes spread across the Chinese territory. As previously mentioned, the inclusion of different lakes allows the forecasting model to generalize its results better, and thus provide more robust outcomes. Another critical difference between these two works is that the former implemented as input information just the spectral bands of the MSI onboard of Sentinel-2, while the latter used both bands and indexes. In this aspect, the present study achieved superior performance, attesting the GMDH model benefited from the inclusion of spectral indexes, which improved their Chla outcomes.

The work conducted by Aranha et al. [13] shares the same location as the one used in this study. However, they used only a subset of 5 reservoirs of 149. In their approach, the authors fitted a regression line to their dataset using the three-band spectral index (3BSI), showing a good agreement between the index and the Chla values. A similar methodology was implemented in reference [152], where the Toming’s index was used to fit a regression line for Chla values. These two studies implemented spectral bands to estimate Chla concentrations. By evaluating the R² metric values for [13,152], the proposed GMDH was superior over both studies, with significant improvements of 12% and 5%, respectively. Furthermore, a major difference between these two studies and the methodology presented in this work, is the data handling. The other studies used methodologies of fitting a regression line using the proposed indexes over the entire dataset, making no distinction between training and testing datasets. This is analogous to assessing our ML model's performance considering the training dataset only. Therefore, their methodologies lack generalization, being bound to a particular time and geographic location. In addition, while these approaches are considerably less complex than the proposed ML model, they provide valuable insights into Chla's behavior when analyzed using the evaluated indexes.

In the work conducted by Pompêo et al. [159], the Sentinel Application Platform (SNAP) algorithm was evaluated for modeling the Secchi disk depth, the Chla concentration and the number of Cyanobacteria cells. The study location was the Cantareira System (CS) in the Brazilian state of São Paulo, Brazil. In their study, authors used in situ collected water samples from 6 reservoirs in the CS as ground truth for the evaluated water quality parameters. These data were latter compared with the Case 2 Regional Coast Color (C2RCC) atmospheric correction algorithm. It is a machine learning paradigm based on neural networks trained to reproduce top-of-atmosphere reflectance [161]. The results of their study showed a good correlation between the modeled data and the real collected samples for the Chla, achieving an RMSE value of 2.3 μg/L and R² of 0.75.

A direct comparison of these results with the ones found in our study shows an R-squared value 16% superior for the GMDH approach, while the C2RCC had better RMSE values. The reasoning behind this discrepancy for the error metric is twofold. Firstly, even though both studies were conducted in Brazilian territory, the state of São Paulo is characterized with having a subtropical and tropical climate type [162]. Such climatic configuration is much less prone to dry seasons compared to the studied semi-arid region of the Ceará state. This consistency leads to less fluctuation in Chla levels. Consequently, this could decrease the model’s variance, thus enhancing the model’s performance. Secondly, the work presented by Pompêo et al. used significantly fewer reservoirs compared to the present work. As previously mentioned, the reduced number of reservoirs hinders the model capacity to generalize unknown data while diminishing the dataset's variance, leading to improved error outcomes. However, despite improved error, the C2RCC model is a less robust approach when compared to the GMDH algorithm, as evidenced by comparing the R² metric of both models.

It is important to disclose that our study was conducted in a tropical semi-arid region, considering 149 dams, considerably larger than any of the presented literature. In addition, a broader investigation of the spectral indexes for the Sentinel-2 constellation was assessed compared to the references in Table 5. This allows us to gather further insight into the impact of the spatiotemporal influence of both bands and indexes on the Chla prediction. From the results for parameter selection in Figure 4 and Figure 5, we observe that bands 8 and 11 do not compose the selected bands set. However, they are still present in the form of indexes. Moreover, a careful assessment of the GMDH results (available at Supplementary Spreadsheet 1) showcases the relevant contribution of bands 3, 4, 5, 7, 8, and 11 to the Chla prediction. This indicates that the predictive models benefit from a higher spatial resolution, as well as from green, red and infra-red bands composing the indexes, according to the literature [13,29,64].

5. Conclusions

This study evaluated several input parameters for Chla modeling in a vast spatial coverage of 149 freshwater reservoirs spanning a semi-arid tropical region across the Brazilian state of Ceará. This evaluation was conducted based on satellite remotely sensed data and ground-truth Chla concentration measurements that reflected the temporal and spatial distribution, notably impacted by interannual rainfall variability. To this end, we investigated the performance of several ML approaches, namely the k-nearest neighbours, random forests, extreme gradient boosting, the least absolute shrinkage and selection operator, support vector machine, and GMDH.

Forward stepwise parameter selection was implemented to determine the best input attribute selection among the 13 bands from the Sentinel-2 constellation and 16 spectral indexes derived from such bands. The usage of this methodology applied to the remote sensing field is, to the best knowledge of the authors, new and proven to improve the outcomes of the investigated Chla forecasting models. Figure 7 elucidates the performance of the GMDH model, showing that the model under investigation yielded excellent values for the Chla concentrations, achieving R² value greater than 90%.

It is important to emphasize that our methodology ensures a clear separation between training and testing data, unlike some studies found in the literature. The approach is substantial for applying ML models, as a way to promote more robust results and capacity to generalize unseen data. This approach is especially crucial for a large dataset such as the one used in the present work, comprising 149 reservoirs in a semi-arid region similar in size to England, where the spatiotemporal variation of the data is substantial. Regarding the spectral bands and indexes, we showed that the Chla modeling benefited primarily from red, NIR, and green bands.

Extensive comparison with the literature found studies showed that the used models in this work offer competitive results. It is important to note that such a comparison may be misleading due to the disparate methodologies and study locations, directly influencing the ML outcomes.

In future works, the addition of more indexes, as well as the merge of Sentinel-2 data with Landsat-8 could be investigated. Including more spectral indexes and the Landsat-8 MODIS data would provide more spatiotemporal information and reduce data variance due to the finer temporal resolution. Furthermore, implementing atmosphere correction preprocessing could also benefit the predictive paradigms, as it would reduce data noise, diminish the error variance, and improve the forecasting of Chla concentration.