Streamflow forecasting in Tocantins river basins using machine learning Open Access

Brazil is the country with the highest amount of fresh water in the world. This resource has been used mainly in agriculture, energy production, economic growth, and sanitation (Cantelle et al. 2018). Due to population growth and economic development, it is estimated that in the last two decades there has been an increase of 80% in total water withdrawal from water bodies, and that there will be an increase of 23% by the year 2030 (ANA 2020). Based on the foregoing, understanding the river system behavior in basins is fundamental to hydroelectricity planning, characterization of grant flows, flood forecasting, and assessment of impacts of climate change and soil use, among others (Bourdin et al. 2012; Tongal & Booij 2018; Yaseen et al. 2018).

The streamflow forecasting and prognoses of hydrological behavior are often performed using hydrological models, which can be classified as conceptual or empirical (Debastiani et al. 2019). Conceptual hydrological models are based on physical characteristics of basins and require large amounts of data, which in many contexts are difficult to be acquired, are unavailable or insufficient for covering all the spatial and time variability (Yang et al. 2019). Empirical models, on the other hand, use a system's data series for mathematical functions to establish connections between the target variable of the estimate and the input variables in the system, disregarding the intervening physical processes (Uliana et al. 2019).

Among the empirical models are the machine learning algorithms, which have been used in basins around the world due to their relatively easy adjustment and practicality, and they are robust tools for analyzing complex systems (Tongal & Booij 2018). Jimeno-Sáez et al. (2018) compared the performance of the Soil and Water Assessment Tool (SWAT) conceptual model to the artificial neural networks for daily streamflow forecasting in the Miño-Sil and Segura watersheds, Spain, and concluded that both are efficient tools for forecasting. Further, they concluded that the conceptual model was superior in the lower streamflows forecasting, whereas the empirical model showed better performance in the highest streamflows forecasting. Adnan et al. (2018) analyzed the least square support vector machine (LSSVM), fuzzy genetic algorithm, and M5 model tree models in the daily and monthly streamflow forecasting in the Hunza River basin, Pakistan, and concluded that these models are efficient for forecasting in the studied river basin, especially the LSSVM model. Kabir et al. (2020) evaluated the wavelet-based artificial neural networks, support vector regression, and deep belief network models for hourly streamflow prediction in three locations in the United Kingdom and observed good performance of the models for a 2-hour forecast horizon.

The use of empirical hydrological models requires the appropriate input variables selection to make the learning process less complex, the interpretation of the results simpler, and reduce the computational cost (Dariane et al. 2019). Because of a large amount of hydrometeorological data available for some basins, the variables selection process is essential to eliminate correlated variables, noisy or non-significant relationships with the dependent variable (Jain & Zongker 1997; Bowden et al. 2005; Prasad et al. 2017; Hadi et al. 2019; Zhu et al. 2019; Afan et al. 2020; Ren et al. 2020).

According to the National Water Agency (2009), one of the threats to the Tocantins-Araguaia hydrographic region is the entry of large enterprises, which can pressure the environment if implemented without proper planning. In view of the need to support the planning and management of water resources in the Sono, Manuel Alves da Natividade and Palma River basins, the objectives of this study were to: (i) select the best sets of hydrometeorological variables using selection algorithms, namely genetic algorithm, recursive feature elimination and vsurf, and (ii) evaluate the performance of the random forest, artificial neural networks, support vector regression, and M5 model tree models in forecasting daily streamflow of the basins.

The study hydrographic basins are the Sono (SRB), Manuel Alves da Natividade (MRB), and Palma (PRB), located in the Tocantins-Araguaia hydrographic region, on the right side of the Tocantins River (Figure 1). The Tocantins-Araguaia hydrographic region has an area of approximately 920,087 km², which starts in the Center-West region and downstream to the bay of Ilha de Marajó, in the north. It is the largest hydrographic region fully inserted in the Brazilian territory (ANA 2015). Due to its great water availability, which is equivalent to approximately 6% of the national total, the region has potential for hydroelectricity, mining, livestock, irrigation, fishing, agriculture, transport, and tourism. The main consumptive water uses are irrigation and human, animal, and industrial consumption (ANA 2009). The SRB (45,042 km²), MRB (14,344 km²), and PRB (17,468 km²) are important hydrological and environmental units of the upper and middle courses of the Tocantins River. Entirely inserted in the Cerrado biome, they also stand out in non-consumptive water uses for the ecosystem's conservation, supply of four hydroelectric reservoirs along the Tocantins River (Peixe Angical, Luís Eduardo Magalhães, Estreito, and Tucuruí), and tourism (Mauriz 2008; ANA 2009; Rodrigues et al. 2020).

Hydrometeorological data were obtained from the Hydrological Information System (HIDROWEB) of the National Water Agency (ANA) and the Meteorological Database for Teaching and Research (BDMEP) of the National Institute of Meteorology (INMET). The codes and location of the stations used for obtaining the daily streamflow series (m³ s⁻¹), average rainfall (mm), and potential evapotranspiration (mm) are presented in Figure 1. In the gauged stations 1, 2, and 3 streamflow series used in the study of SRB, MRB, and PRB, respectively, were monitored. The spatial average rainfall in the basins was obtained by the Thiessen Polygon method (Macêdo et al. 2013), considering the rainfall stations n.6 to n.11, n.13 and n.15 for SRB; n.11 to n.15 for MRB; and n.15 to n.20 for PRB. The weather station data n.4 was used for SRB and n.5 for MRB and PRB for calculating the evapotranspiration by the Penman-Monteith equation (Allen et al. 1998).

As input variables in the selection and modeling process the Julian day (JD), ordered from January 1, and the observed values in the previous 3 days of the respective streamflow (Q_t−1, Q_t−2, Q_t−3), average rainfall (R_t−1, R_t−2, R_t−3) and evapotranspiration (ET_t−1, ET_t−2, ET_t−3), totaling 10 variables were considered. The use of lagged variables in time, mainly from the rainfall series (R_t−1, R_t−2, R_t−3), is important in empirical models, because, among other reasons, it covers the time of concentration of the basins in the modeling, which for this study are 63.35 (SRB), 37.82 (MRB) and 61.72 hours (PRB) according to the Giandotti equation (Giandotti 1940).

The periods available and used in the evaluations were from 7/1/1977 to 5/28/2017, 12/30/2017, and 9/25/2006 for SRB, MRB, and PRB, respectively. Because of the existence of gaps in some series and the large amount of observed data, dates with missing data were removed from the database, leaving a total of 10,194 (SRB), 12,019 (MRB), and 5830 observations (PRB). The dataset was divided into training (70%) and validation (30%) sets to assess the efficiency of the models. All statistical procedures were performed in the R Statistical Environment (R Development Core Team 2019).

The feature selection algorithms used (genetic algorithm, recursive feature elimination, and vsurf) belong to the wrapper methods, whose approach consists of the search for subsets of variables that minimize the estimation error of a machine learning algorithm (Kohavi & John 1997). The random forest (RF) machine learning model was used in all feature selection approaches due to its accurate forecasting and resistance to overfitting (Breiman 2001; Cutler et al. 2012). The selected algorithms are presented below.

Inspired by evolutionary biology, the GA for feature selection considers a set of variables to be a chromosome/individual, and a set of individuals to be a population. By a fitness function, the algorithm evaluates and ranks the individuals according to a specific metric. The ordered individuals go through genetic operators (elitism, selection, crossing, and mutation) for creating new individuals that will integrate to a new generation, which will be evaluated and ordered by the fitness function recursively until the stopping criterion is met and the best individual is selected (Holland 1975; Bento & Kagan 2008; Xue et al. 2018). The algorithm was fitted using the gafs function of the ‘caret’ package (Kuhn et al. 2019), with the RF model as fitness function and the root mean square error (RMSE) as the objective function to be minimized. The internal parameters used in the GA configuration were defined through a trial-and-error process, guided by related studies (Dariane & Azimi 2016, 2018; Afan et al. 2020), and the parameters are the number of generations (20), population size (30), and probability of crossing (0.8), mutation (0.1), and elitism (0.0).

The RFE algorithm consists of the retroactive process of feature selection, which occurs in the fitting of a machine learning model, allocation importance to the available features, and removing the least relevant feature in the predictive process. The recursive elimination of features is gradual, starting from the total set of features for a predetermined maximum number (Guyon et al. 2002; Granitto et al. 2006). The execution of the algorithm was performed using the rfe function of the ‘caret’ package (Kuhn et al. 2019), with the maximum number of features regressing from 10 to 1.

The VS selection algorithm operates in three steps: the first one is the classification and ordering of features as to their importance, according to the results of the RF model, eliminating those considered irrelevant to the process; the second one selects the features related to the dependent features that lead to lower forecasting error; and in the third step there is the fitting of models with the sequential introduction of the features selected in the previous step, evaluating the reduction of the error for the decision regarding the permanence of each one (Genuer et al. 2015). The algorithm was implemented using the ‘VSURF’ package (Genuer et al. 2019), with 500 trees developed in the RF (ntree) and the standard number of features randomly sampled for the growth of each tree (mtry). The selected features in the third step were adopted.

Machine learning algorithms must be able to generalize their estimates from the input data, and an important step is the appropriate choice of parameters of the model. One of the main methods for selecting the best model with optimization of parameters is the cross-validation, applied in the models training step (Elshorbagy et al. 2010; Shortridge et al. 2016; Yang et al. 2017). The cross-validation k-fold was implemented in this study with the ‘caret’ package (Kuhn et al. 2019), with 10 resamples. The evaluated models were: random forest, artificial neural networks, support vector regression, and M5 model tree.

Introduced by Breiman (2001), the RF is a classification or regression algorithm that aggregates the results of a set of decision trees, grown with random variables/features, independently sampled, and equally distributed. For regression problems, the result of the model is given by the average of the predictions made by all the trees in the forest (Breiman 2001). The RF models were fitted using the rf function of the ‘randomForest’ package (Liaw & Wiener 2002), in which the optimized parameter was mtry, which varied from 2 to 10, with ntree equal to 500.

ANN are algorithms inspired by the human nervous system. The basic structure has an input layer containing the independent variables, one or more hidden layers for processing the data, and an output layer with the results of the iterations. The neurons are interconnected by weights so that the receiving neuron aggregates the weights of the previous layer, adds a bias, and forwards the result through a transfer function (Dastorani et al. 2018; Dariane et al. 2019). To fit the ANN, the nnet function of the ‘nnet’ package were used (Venables & Ripley 2002) with feedforward multilayer perceptron architecture, one intermediate layer, and linear output units. The input data were normalized by the min-max method to be between 0 and 1 (Riad et al. 2004; Selvi & Huseyinov 2020). The optimized parameters were the number of neurons in the hidden layer (size), from 1 to 10, and the decay of the weights (decay), used mainly to avoid overfitting the model to the data, considering 0.00, 0.01, 0.05, 0.10, 0.50, 1.00, and 2.00.

The SVR, introduced as support vector machine for classification by Vapnik (1995), is the projection of the features in high-dimensional space to map data sets by fitting a curve (kernel function) between two marginal hyperplanes, to minimize the regression error. Among the available kernel functions are linear, polynomial, radial, and sigmoid (Vapnik 1995; Yang et al. 2017; Dastorani et al. 2018; Shamshirband et al. 2020). The function used to fit the SVR was the svmRadialSigma from the ‘kernlab’ package (Karatzoglou et al. 2004), in which the kernel function used is radial and the parameters regularization (C) and kernel function (sigma) are required. The values for C optimization ranged from 1 to 20 and were tested for sigma 1 10⁻⁵, 1 10⁻⁴, 1 10⁻³, 1 10⁻², and 1 10⁻¹.

The M5T algorithm, proposed by Quinlan (1992), presents a tree structure developed in two stages. The first step consists of growing a decision tree using independent variables to homogenize the responses to the variable of interest, using the reduction of the standard deviation as a criterion for dividing the nodes. Successive divisions generate excessive branches, which may cause overfitting of the data. Therefore, the second step is pruning the tree grown by replacing the branches with linear regression functions according to the independent variables used in the divisions of the pruned branch. Thus, for each homogenized subspace of the target variable, a linear model is fitted (Quinlan 1992; Pal & Deswal 2009; Shamshirband et al. 2020). The ‘Cubist’ package (Kuhn et al. 2020) was used to fit the models, with the number of interactive model trees (committees) ranging from 10 to 100 adding 10.

The NSE proposed by Nash & Sutcliffe (1970) is one of the main metrics for evaluating the performance of hydrological models. The classification of the performance of the fitted models with data on a daily scale proposed by Moriasi et al. (2015), suggests: NSE>0.80 as ‘very good’; 0.70<NSE≤0.80 ‘good’; 0.50<NSE≤0.70 ‘satisfactory’; and NSE≤0.50 as ‘not satisfactory’.

Table 1 presents the sets of total independent variables/features and the selected ones by the selection algorithms for each basin, with the latter being ordered in descending order according to the importance attributed to them during the selection processes. The VS algorithm selected a smaller number of variables, followed by the GA, which removed one variable from the total set of the SRB and MRB and four from the PRB, and RFE which selected nine in each basin. The variables Q_t−1 and P_t−1 were selected by all algorithms, followed by Q_t−2 and ET_t−1, selected in eight of the nine selections (3 basins×3 selection), and JD selected in seven. The variables that were more frequently removed were ET_t−2 and ET_t−3.

The most important variables by most algorithms were the Q_t−1 and Q_t−2 streamflows, with only the selection made by the GA in the PRB presenting the Q_t−3 as the second most important variable. Previous streamflows are widely used in empirical hydrological models due to their strong correlations with the target streamflow (Yang et al. 2017; Liu et al. 2019; Yaseen et al. 2019; Zhu et al. 2019). Following the streamflows, overall, the JD and the lagged average rainfall were the most important, since the JD temporarily marks the changes in the streamflow, and rainfall provides the main water entry into the basins (Dinpashoh et al. 2019).

The model parameters, optimized during the cross-validation process, are shown in Table 2, and the results of the performance of the models using different sets of variables, in the training and validation periods, for the SRB, MRB, and PRB are present in Tables 3–5, respectively. The models' performances were classified as very good, according to the classification proposed by Moriasi et al. (2015), with NSE ranging from 0.849 to 0.926 in the validation. Only the ANN model, fitted with the set of variables selected by VS for PRB presented a performance classified as good (0.769). The results of the NSE elucidate the precision of the models in the peak flows forecasting, which has been a challenge for hybrid hydrological models, such as SWAT, LASH, MGB-IPH among others. Beyond this, the LNSE values ranged from 0.782 to 0.969, revealing a good accuracy in the forecasting of the baseflow in the basins.

The models presented better performance when fitted with all available variables or with those selected by the GA and RFE algorithms, corroborating the increasing use of selection algorithms in hydrological studies. Dariane & Azimi (2016) evaluated the GA to select input variables in the ANN model, which was applied for forecasting the monthly streamflows in the Ajichai sub-basin, Iran, and obtained superior performance using the selected variables (NSE=0.870) instead of all the variables (NSE=0.780). Khan et al. (2020), for forecasting droughts in Pakistan, used the RFE to select variables for the ANN, support vector machine, and k-nearest neighbor models and obtained results considered satisfactory. Dariane et al. (2019) using ANN for forecasting rainfall at six stations in Iran, evaluating the selection methods GA, wingamma, and self-organized map, concluded that the GA is the most reliable for input feature selection.

According to the results presented in Tables 3–5, the performance of the M5T model using the set of variables selected by the RFE was superior to the other models for SRB and PRB, with lower RMSE (99.254 and 59.671 m³ s⁻¹) and higher NSE (0.926 and 0.915) and LNSE (0.959 and 0.964) in the validation periods, respectively. For MRB, the total set of input variables provided the best fit of the M5T model, considered the best for the basin (RMSE=58.039 m³ s⁻¹, MAE=22.857 m³ s⁻¹, NSE=0.899 and LNSE=0.969).

Yin et al. (2018) fitted the SVR, multivariate adaptive regression splines, and M5T models to streamflow forecasting in the Pailugou River basin, China. They concluded that the M5T model was superior to the others, with an average NSE of 0.890 in the validation. Diverging from this study, Adnan et al. (2019) observed, when evaluating six models in the daily forecasting of the monitored streamflows in two stations on the Fujiang River, China, lower performance of the M5T model, with average NSE values of 0.652 and 0.670, evidencing the need to study the models for different basins. Tongal & Booij (2018) verified the accuracy of the SVR, ANN, and RF models for four river basins in the United States, and they observed that the NSE varied between 0.880 and 0.980 in the validation, which was evaluated as very good by the researchers and similar to those obtained in this study. Liu et al. (2019) fitted the Gaussian mixture regression, SVM, and ANN models for forecasting streamflows in two sections of the Jinsha River basin, China, and obtained NSE indexes between 0.850 and 0.860, corroborating the results found for the river basins in this study. Regarding the above results, the good performance of statistical models fitted in this study was confirmed, especially the M5T model.

Figure 2 shows the hydrographs observed and simulated by the models using the variables that provided the least statistical errors for the SRB, MRB, and PRB basins in the validation. Good adherence of the forecasted hydrographs to those observed can be seen, mainly in the periods of streamflow recession. According to Rodrigues et al. (2021), the best models for these periods should be based on the relationship between the groundwater flow, which is predominant in the recession periods, and the discharge from the aquifer. Peak flows were underestimated by the models, mainly due to the difficulty in capturing the many variables involved in the streamflow process, the difficult representation of the spatial variability of rainfall over the basin due to the low density of rain-gauge stations in the basins, and even errors from the extrapolation of the stage-curve of the river (Viola et al. 2009; Silva Neto et al. 2020).

Figure 3 presents the flow duration curves (FDC) of the observed streamflows and forecasted by the models using the variables that provided the best fitting for each river basin, considering the validation period. This analysis allows inferring about the percentage of time in which a given streamflow is exceeded or equalized and to establish reference values for water management and projects. The FDC forecasted by the M5T and ANN models presented a great agreement with the observed FDC for the SRB and MRB, as well as by the M5T and RF models for the PRB. The inferior performances were observed for the RF and SVR models for MRB, and ANN for PRB, mainly with overestimates of the lowest streamflows.

Table 6 shows the observed streamflows with 95% (Q₉₅), 90% (Q₉₀), 50% (Q₅₀), 10% (Q₁₀), and 5% (Q₅) of exceedance and forecasted by the models in the validation periods. The Q₉₅ and Q₉₀ streamflows indicate the values that are equaled or exceeded for 95 and 90% of the time, respectively, and have been used as a reference for water rights in several Brazilian states and for planning irrigation systems, dams, and hydroelectric plants since they represent the minimum flows (recession period). The estimated Q₉₅ were close to those observed for the SRB and PRB, with errors ranging from 0.07 to 28.53%, while for MRB the best forecasts were obtained by the ANN and M5T models, with errors of 10.69 and 7.64%, respectively.

According to Decree No. 2432 of June 6, 2005, the Q₉₀ streamflow is the reference used for granting in the Tocantins state. The streamflow forecast with 90% of exceedance was performed satisfactorily by the models, highlighting M5T. Such value was forecasted for SRB as equal to 260.15 m³ s⁻¹ (observed=258.27 m³ s⁻¹) and for MRB, 23.16 m³ s⁻¹ (22.78 m³ s⁻¹). The RF model forecasted this streamflow as 119.98 m³ s⁻¹ for PRB (observed=119.85 m³ s⁻¹). These results, together with those found for Q₉₅, showed the precision of the streamflow forecasting by the models for the recession periods and its potential use as a tool in the management of water resources.

The other reference streamflows (Q₅₀, Q₁₀, and Q₅) represent the values with medium and highest values (low durations) in the basins. These streamflows (Q₁₀ and Q₅) were underestimated by most models; however, the forecasts for this complex hydrological process were satisfactory, with errors between 0.14 and 8.57% for Q₁₀ and −0.34 and −7.95% for Q₅. The forecast for the Q₅₀ streamflow resulted in average errors of 4.18% (SRB), 7.25% (MRB), and 9.88% (PRB), with better performances of the M5T model for the SRB and PRB, and ANN for the MRB. Accurate forecasting of peak flows is necessary to plan actions to mitigate floods and support the drainage and dam projects (Tucci 2004).

Rodrigues et al. (2021) evaluated the performances of the SWAT and ANN models for forecasting daily streamflow of the MRB and obtained NSE indexes in the validation step of 0.610 and 0.910, respectively, with errors of −14.10% in estimating Q₉₀ using SWAT and −10.60% by ANN. For the same basins, Rodrigues et al. (2020) found, by the SWAT model, NSE values of 0.730 (SRB), 0.810 (MRB), and 0.700 (PRB) in the validation period. In view of these results for the same basins in this study, the good performance of the evaluated models was confirmed, since they presented results superior to those of the conceptual SWAT model and values close to that of the empirical ANN, mainly the M5T model using the sets of variables selected by the RFE for SRB and PRB and total set for MRB. However, although machine learning algorithms showed good streamflow simulation capabilities, the models are based on pattern recognition and do not consider the hydrological processes involved. Therefore, although new data is incorporated, it only operates according to previously learned patterns and is not recommended for situations in which there are alterations of the physical characteristics of the basin (Rodrigues et al. 2021).

The analyzed models proved to be suitable for application in the management and planning of water resources. Considering the short forecast horizon (1 day) the models can be useful tools in the operation of hydraulic works, river navigation, irrigation, hydroelectric power generation, water supply, flood control, among other purposes requiring short-time streamflow predictions in the SRB, MRB, and PRB basins.

The random forest, artificial neural networks, support vector regression, and M5 model tree models, fitted with the total set of variables and sets selected by the genetic algorithm, recursive feature elimination, and vsurf, showed good performances in the daily streamflow forecasting of the Sono, Manuel Alves da Natividade and Palma River basins.

The performance of the M5 model tree model stood out from that of the other models in the calibration and validation periods using the variables selected by the recursive feature elimination algorithm for the Sono (Q_t−1, Q_t−2, JD, Q_t−3, R_t−1, R_t−2, ET_t−1, ET_t−2, R_t−3) and Palma (Q_t−1, Q_t−2, Q_t−3, JD, P_t−1, ET_t−3, R_t−2, ET_t−1, R_t−3) River basins, and total set of variables for the Manuel Alves da Natividade River basin (JD, Q_t−1, Q_t−2, Q_t−3, R_t−1, R_t−2, R_t−3, ET_t−1, ET_t−2, ET_t−3), considered the best sets of variables for fitting the M5 model tree algorithm in the studied basins.

The hydrographs, FDC, and forecasted reference streamflows confirmed the satisfactory performance of the models, considering the agreement between the hydrographs (observed and forecasted) and the forecasted and observed FDC, and the small estimation errors in the peak and recession streamflows periods. In view of these results, it can be concluded that the evaluated models are potential tools to assist in the different demands faced by the water resources management in the river basins, especially the M5 model tree model using streamflow lags, average rainfall, and evapotranspiration as inputs.

The main limitation of this study was the difficult representation of the spatial variability of rainfall over the basin due to the low density of rain-gauge stations in the basins. Due to the successful application of the machine learning algorithms in this study, it is suggested for further studies cover a longer forecast horizon.

This study was financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior – Brasil (CAPES) - Finance Code 001; The second author would like to thank the Conselho Nacional de Desenvolvimento Científico e Tecnológico – CNPq for granting fellowship of research productivity (PQ) - number 311191/2021-5. The authors would like to extend thanks to the National Water Agency (ANA) and the National Institute of Meteorology (INMET) for providing the data for implementing this study.

All relevant data are included in the paper or its Supplementary Information.