Machine learning models coupled with empirical mode decomposition for simulating monthly and yearly streamflows: a case study of three watersheds in Ontario, Canada

Abstract

This paper presents a novel approach for enhancing long-term runoff simulations through the integration of empirical mode decomposition (EMD) with four machine learning (ML) models: ensemble, support vector machine (SVM), convolutional neural networks (CNN), and artificial neural networks with backpropagation (ANN-BP). The proposed methodology uses EMD to decompose precipitation and temperature time-series into intrinsic mode functions, thereby revealing underlying data patterns. Subsequently, these components are incorporated into the ML models to simulate the runoff time-series. The effectiveness of the hybrid models is evaluated using streamflow runoff data obtained from the Grand, Winnipeg, and Moosonee Rivers in Ontario, Canada. Four widely used performance indices, namely, correlation coefficient, root mean square error (RMSE), mean absolute relative error, and Nash–Sutcliffe efficiency, are employed to assess the models’ performance. The results demonstrate that the hybrid EMD-ML models exhibit significantly superior performance compared with the standalone ML methods. During the validation phase, the EMD-Ensemble, EMD-SVM, EMD-CNN, and EMD-ANN-BP models exhibit notable reductions in the RMSEs of monthly streamflow estimates for the Grand River, amounting to 11%, 22%, 8%, and 33%, respectively, compared with their non-EMD counterparts. Additionally, these hybrid models exhibit improved RMSEs for yearly simulations in the Winnipeg River, with reductions of 54%, 0.08%, 6%, and 4.5% respectively. To further enhance the accuracy of monthly and yearly streamflow estimates, an SVM-recursive feature elimination technique is employed to select a more appropriate EMD dataset in all study cases. This research underscores the potential of integrating EMD with ML models to enhance long-term runoff simulations. The outcomes highlight the superior performance of the hybrid EMD-ML models, demonstrating their ability in generating lower biases than the standalone ML methods. These findings hold significant implications for the field of computational fluid mechanics and can contribute to the understanding of hydrological processes.

KEYWORDS:

• Machine learning
• runoff estimation
• EMD method
• ensemble model

Introduction

A key task in hydrological modelling studies is to accurately simulate the runoff in a basin over a long-term period, for example, monthly and yearly time-series, which affects the water supply efficiency, flood management, and water resource management (Mohammadi, 2021; Soltani et al., 2021). Several methods have been applied to simulate the runoff, including conceptual models, machine learning (ML) models, and coupled models (Kratzert et al., 2018). In particular, several researchers have developed time-series statistical models such as moving average (MA), autoregressive (AR), and/or autoregressive moving average for runoff time-series simulation (He et al., 2019; Weeks & Boughton, 1987). Recently, ML and deep learning (DL) models have been applied for runoff modelling (Hu et al., 2018; Mallick et al., 2022; Parisouj et al., 2020, 2022).

ML is a broad field that covers artificial intelligence, probability, psychology, and statistics, among other domains. ML methods can be used to simplify and solve problems (Nasteski, 2017). Recently, DL methods, as a subset of ML methods, have received considerable attention in hydrology and community studies (Barzegar et al., 2021; Goliatt et al., 2021; Nasteski, 2017). Artificial neural networks (ANNs) are widely used in hydrology for runoff simulation. However, their learning process is slow, and they are prone to being trapped in local minima (Parisouj et al., 2020). To address this limitation, Cortes and Vapnik (1995) proposed the support vector machine (SVM) technique to efficiently find global optimum solutions (Meng et al., 2019). Several researchers have highlighted the superiority of support vector regression (SVR) over artificial neural network with backpropagation (ANN-BP) for simulating runoff (Bafitlhile & Li, 2019; Kalteh, 2013). DL architectures, such as long short-term memory (LSTM), deep neural networks, and convolutional neural networks (CNNs) are being widely used to predict the time-series of wind, solar, and streamflow (Ghimire et al., 2021). Among these, the CNN algorithm lacks post-processing capability (Liu et al., 2022).

Other promising ML models include ensemble models such as extra tree regressor, bagging regressor, adaptive boosting regression (AdaBoost), and stack generalisation. Ensemble models combine multiple models, and the output of one model is input to another model to increase the prediction accuracy (Sagi & Rokach, 2018). Recently, ensemble models have been applied for predicting rainfall runoff (Barrera-Animas et al., 2022; Jose et al., 2022; Zhao et al., 2022). Tarfaya et al. (2022) used ensemble models for predicting the index rainfall and showed that the extra tree model yielded reasonable results. Liu et al. (2014) applied AdaBoost to enhance the accuracy of runoff prediction. Elbeltagi et al. (2022) predicted the river flow rate in the Moines watershed by applying ML models, and the bagging model was noted to exhibit acceptable performance.

The prediction accuracy of ML models can be enhanced by extracting trends and harmonics from hydrological time-series and removing noise through appropriate data preprocessing techniques, such as genetic algorithm optimisation, MA, principal component analysis, singular spectrum analysis, wavelet analysis, and gamma testing (Band et al., 2021; Bartoletti et al., 2018; Cui et al., 2021; Golshan et al., 2020).

Huang et al. (1998) proposed the empirical mode decomposition (EMD) technique for noise-assisted data analysis. Several researchers have applied EMD to extract signals from noisy nonstationary data in the analysis of several aspects, such as hydroclimatic processes, solar radiation, wind speed, speaker recognition, and ice-snow coverage (Lee & Ouarda, 2012; Metzger et al., 2020; Prasad et al., 2019; Sánchez-Martínez et al., 2022; Zhang et al., 2018). The EMD technique is self-adaptive, which makes it preferable to other traditional approaches (Tayyab et al., 2018). Temperature, precipitation, and runoff represent nonlinear and nonstationary time-series (Chen et al., 2018). Therefore, the EMD approach can be used to analyse these hydrological time-series data.

The objective of this study was to use the novel hybrid EMD-Ensemble method to estimate runoff at monthly and yearly time scales. To the best of our knowledge, none of the existing studies have used the combined EMD-Ensemble approach to estimate monthly and yearly streamflows. The hybrid EMD-Ensemble was compared with hybrid EMD-SVR, EMD-CNN, and EMD-ANN-BP. In addition, the performance of hybrid models was compared with the standalone ensemble, SVR, CNN, and ANN-BP models. Three monthly and annual runoff time-series from the Grand, Winnipeg, and Moosonee Rivers in Canada were investigated to ensure the applicability of the proposed framework. These three rivers were chosen because of differences in their drainage area (watershed), landcover, discharge, watershed topography, and climate.

Method and materials

Study area and data description

The study areas included the Grand, Winnipeg, and Moosonee River basins in Ontario, Canada, which have different sizes, geology, and hydro climatology. Figure 1 shows the locations of the three basins in Ontario.

Figure 1. Study area.

The Grand River watershed, with a drainage area of 6965 km² and length of 280 km, is located in southern Ontario. Bahamonde et al. (2015) reported that 76% and 17% of the watershed pertain to agricultural land and forest areas, respectively. Thirty treated effluents from municipal wastewater are dumped into the Grand River basin. The elevation difference in the regions upstream and downstream of the Grand River basin is approximately 350 m. The northwest region receives more rain than the southeast, with annual precipitation levels averaging 850 mm with a peak of 1,000 mm. January and February are the driest months, and July and August are the wettest months. Additionally, the mean annual temperature varies from 5°C in the higher elevations in the north to 8°C along the lakeside (Krause et al., 2001). Great Lakes, the Arctic region, and the Gulf of Mexico affect the climate of the Grand River basin.

Most of the Winnipeg River watershed, with a drainage area of approximately 150,000 km², is located in the northwestern Ontario province, with part of it lying it the southeastern Manitoba province. The basin’s main branch passes predominantly through forest areas between the northern United States and northern Canada. The Winnipeg River flows to the west, through Manitoba province, into lake Winnipeg. The elevation difference in the regions upstream and downstream of the Winnipeg River basin is approximately 217 m. Approximately 30% of the precipitation is snowfall, with annual precipitation averaging 780 mm. Similar to the Grand River basin, January and February are the driest months, and June and July are the wettest months. The annual mean temperature varies between −18°C and 25°C at Slave Falls (St. George, 2007).

The Moosonee River watershed, with a drainage area of approximately 109,000 km², is located in northeast Ontario, southwest of the James Bay region. The watershed consists of three main tributaries, including the Missinaibi and Mattagami Rivers that constitute the Moose and Abitibi Rivers. The basin is divided into northern and southern portions. This study focused on the southern portion, which is more topographically diverse. The elevation difference in the regions upstream and downstream of the Moosonee River basin is approximately 580 m. The mean annual precipitation varies from 650 to 1,000 mm, and the mean annual temperature is between −9°C and 8°C. The correlation between precipitation and temperature is typically positive, but the wettest period is July–December. Approximately 35% of the mean annual precipitation is snowfall, and more than half of the annual Moosonee River watershed discharge occurs in spring (Ho et al., 2005; Story & Buttle, 2001).

The daily discharge data from January 1, 1973, to December 31, 2020, was selected from the Water Survey of Canada for three measurement stations including the watersheds (https://wateroffice.ec.gc.ca/). A daily dataset including precipitation (P), maximum temperature (T_max), and minimum temperature (T_min) was selected, consistent with the discharge period from the nearest weather stations. Records of weather observations are available at https://climate.weather.gc.ca/. Figure 1 shows the locations and information of the streamflow and weather stations.

Methodology

EMD model

Huang et al. (1998) developed the EMD algorithm for nonlinear and nonstationary datasets. The core of this algorithm is that most raw time-series have multiple frequencies of different scales (Karthikeyan & Nagesh Kumar, 2013). The algorithm defines the decomposition dataset to a group of bands by several intrinsic mode functions (IMFs) and a residual function that is calculated for all IMFs and remains constant. Owing to its simplicity, EMD has been widely applied in hydrology (Karthikeyan & Nagesh Kumar, 2013). Kamath and Senapati (2021) coupled the EMD model with ANN to predict 24-h wind speed. The results revealed that EMD-ANN outperformed ANN. Sibtain et al. (2021) predicted the runoff and explored the effect of EMD on the ANN-BP model. The EMD-ANN-BP was noted to outperform the standalone ANN-BP. Yuan et al. (2021) incorporated ensemble EMD (EEMD) in an LSTM to predict daily runoff and proved that the use of the EEMD model output as input data to the LSTM model enhanced the model performance. In general, EMD intermittently extracts the scales of the IFM time-series. The IMFs must satisfy two requirements:

1. The absolute value of the total local minima and local maxima minus zero-crossing equals zero or one.

2. The mean between the upper and lower envelopes in each time-series point equals zero.

To satisfy these requirements, the IMF must be generated by making the function smooth with respect to zero. The EMD breaks down time-series into IMFs through ‘sifting’. Chu and Huang (2020)_ENREF_11 demonstrated that EMD can satisfactorily treat nonstationary time-series in hydrological analyses. The process flow of the EMD algorithm is illustrated in Figure 2.

Figure 2. Process flow of empirical mode decomposition (EMD).

SVR model

SVR, a ML model introduced by Cortes and Vapnik (1995), has been widely applied in hydrology (Achite et al., 2022; Kolachian & Saghafian, 2021; Mozaffari et al., 2022; Sun et al., 2021). Unlike methods such as ANN, which implement empirical risk minimisation, SVR implements the concept of structural risk minimisation. The SVR algorithm for regression computation can be expressed as follows: (1) $f\left(x\right)=a.k\left(x\right)+b$(1) where $a$ represents the weight vectors or coefficients, $k\left(x\right)$ is the kernel function, and $b$ is a bias term. In this study, the radial basis function (RBF) kernel is used to solve Equation (1). (2) $\begin{array}{rl}& \text{maximize}\{\begin{array}{rl}& -\frac{1}{2}\sum _{i,j=1}^l\left({\alpha }_i-{\alpha }_i^{\ast }\right)\left({\alpha }_j-{\alpha }_j^{\ast }\right)k(x_i-x_j)\\ & +\sum _{i,j=1}^l{y}_i\left({\alpha }_i-{\alpha }_i^{\ast }\right)\end{array}\end{array}$(2) (3) $\begin{array}{rl}& \text{subject }\text{to}\{\begin{array}{rl}& \sum _{i=1}^l\left({\alpha }_i-{\alpha }_i^{\ast }\right)=0\\ & \sum _{i=1}^l\left({\alpha }_i-{\alpha }_i^{\ast }\right)\le C\nu l\\ & \begin{array}{cc}{\alpha }_i,{\alpha }_j& \in \left[0,C\right]\end{array}\end{array}\end{array}$(3) where $l$ is the sample size, $\alpha$ and ${\alpha }^{\ast }$ are Lagrange multipliers, $C$ is the cost of the kernel function, $y_i$ is the output, and $k(x_i-x_j)$ is the kernel function. Equation (2) must satisfy the Karush–Kuhn–Tucker requirements to map the dataset, which can be defined as follows: (4) $\{\begin{array}{rl}& {\alpha }_i^{\ast }\left(f\left(x_i\right)-y_i-\eta -{\omega }_i^{\ast }\right)=0\\ & {\alpha }_i\left(y_i-f\left(x_i\right)-\eta -{\omega }_i\right)=0\\ & {\alpha }_i^{\ast }{\alpha }_i=0;{\omega }_i^{\ast }{\omega }_i=0\\ & \left(C-{\alpha }_i^{\ast }\right){\omega }_i^{\ast }=0;\left(C-{\alpha }_i\right){\omega }_i=0\end{array}$(4) where $\eta$, ${\omega }_i^{\ast }$, and ${\omega }_i$ are slack variables. Finally, the SVR can be solved using the following set of equations: (5) $\begin{array}{rl}& \sum _{i,j=1}^l\left({\alpha }_i-{\alpha }_i^{\ast }\right)\cdot k(x_i-x_j)+b,\end{array}$(5) (6) $\begin{array}{rl}& b=y_i+\eta -\sum _{i,j=1}^l\left({\alpha }_i-{\alpha }_i^{\ast }\right)\cdot k(x_i-x_j)+b\end{array}$(6)

CNN model

ML research focuses on DL modelled on the human brain. ANNs are representative computational systems in this context. Neural networks in DL must train computers to have the same functionality as that of the neural system in the human brain. LeCun et al. (1998) designed CNNs, which are DL models that have been widely used in classification and regression tasks in several fields, especially hydrology (Hussain et al., 2020; Sadeghi et al., 2019; Tu et al., 2021). Unlike traditional neural networks, CNNs incorporate multiple architectures such as pooling, local connections, and shared weights. CNNs work on the principle that the input dataset consists of images or data that can be represented as images. Consequently, the processing time and number of parameters are reduced. A CNN typically includes convolutional layers, pooling layers, and fully connected layers. Convolutional layers, as key components, include filters known as kernels that apply convolutional functions to the input dataset and prepare pixels for the next process. Pooling layers help the CNN model control overfitting, thereby limiting the required computation and parameters by reducing the representation size in convolutions (Tu et al., 2021). The LeakyRelu function can accelerate the convergence of the CNN model and facilitate the learning of the neuron weights, even if the input includes zero values. Several structures have been developed based on the type of input data and research objectives, such as InceptionV3 (Szegedy et al., 2016), VGG16 (Simonyan & Zisserman, 2014), ResNet50 (He et al., 2016), Xception (Chollet, 2017), and InceptionResNetV2 (Szegedy et al., 2017), and the corresponding layers, learning parameters, and training process have been elucidated.

ANN-BP model

ANN-BP is a three-layer traditional ANN with feedforward that uses BP in the training dataset (Parisouj et al., 2020, 2022). ML algorithms can be divided into supervised and unsupervised methods. The ANN-BP is a supervised learning method in which data passes through the input layers. The weights minimise the error through hidden layers and add biases to the calculation. The output is generated by the output layer (Sudheer et al., 2002). BP is implemented by assigning a gradient of the loss function to each weight in the chain rule. The gradient is computed one layer at a time, iterating backward from the last layer to avoid superfluous intermediate term calculations. Stochastic gradient descent is a representative learning algorithm that uses BP to compute a gradient.

Ensemble model

Ensemble learning is an ML algorithm that uses a group of base learners to assess and solve real-world issues. Meta-learning pertains to learning from base learners, and ensemble ML (EML) techniques are meta-learning methods (Sagi & Rokach, 2018; Tyralis et al., 2021; Zhang & Ma, 2012) that merge two or more models to enhance the generalisation and performance.

AdaBoost: AdaBoost is an effective ensemble technique that adapts a sequence of base learners to a larger dataset using a more recent dataset. To reproduce the final prediction, base learner predictions are merged using a weighted summation (Idris et al., 2012). The training set is enhanced by the addition of weights ${\omega }_1$, ${\omega }_2$, … , ${\omega }_N$ in each boosting iteration. The initial boosting iteration uses the same weights and data. Subsequently, the learner algorithm is applied to the newly weighted data. In subsequent rounds, the weights of the incorrectly (correctly) predicted training data are increased (decreased). Eventually, each poor learner is forced to focus on the samples missed by the previous learners (Liu et al., 2014).

Bagging regressor: The bagging regressor is a bootstrap aggregation-based ensemble meta-estimator. There exist m bootstrap copies of a sample data point drawn with replacement, and the base learner is used for each bootstrap sample. Finally, the results of each base learner are averaged or voted on. In certain cases, the base learner is a regression or a classification algorithm. Aggregation helps reduce the variance of an individual base learner (Breiman, 1996; Meddage et al., 2021; Singh et al., 2022). Notably, research on runoff prediction using the bagging regressor is limited at present.

Extra tree regressor: Geurts et al. (2006) presented the extra tree model as a decision tree based on an ML algorithm. The model builds a group of decision trees that do not have to be pruned (the trees grow in a top-to-down configuration). The advantages of the extra tree model can be summarised as follows: (1) can easily avoid overfitting, (2) is robust to noise, and (3) can efficiently handle high-dimensional data without feature selection. Similar to other tree-based ensemble approaches, the extra tree regressor produces a collection of decision trees but emphasises randomisation to reduce variance without increasing bias (Eslami et al., 2020; Geurts et al., 2006). The extra tree regressor technique generates random split nodes, allowing it to be implemented faster than other decision-tree-based approaches. To prevent any subsequent increase in bias, three parameters are imported: (1) number of randomly selected attributes at each node (random state), (2) minimum number of samples required to split an internal node (min samples split), and (3) number of trees in the forest (n estimators).

Stack generalisation: Wolpert (1992) introduced the stacking generalisation method, in which several models are assembled to develop an efficient meta-learner. This model takes advantage of singular models to enhance the generalisation. Several models are used as estimators, and one model is used as the final estimator. The results of estimators are used as the input of the final estimator. In this manner, the stacking generalisation method can enhance the performance of the final estimator model. To avoid overfitting, the meta-model does not directly learn the outputs of the base models. The model can be mathematically expressed as follows: (7) $\hat{y}\left(x\right)=\sum _{i=1}^m{\omega }_i{h}_i\left(x\right)$(7) where ${\omega }_i$ is the weight determined for each base learner, and $h_i\left(x\right)$ is the model prediction.

The optimal final prediction is defined by minimising the set of stacking weights using mean square linear regression. Equation (8) mathematically represents the least-squares regression: (8) $\{\begin{array}{rl}& {\omega }^{\ast }=argnmi{n}_{\omega }\sum _{j=1}^n{\left(y\right(x_j)-\sum {\omega }_i{h}_i^{\left(-i\right)}(x_j\left)\right)}^2\\ & {\omega }_i\ge 0\\ & \sum _{i=1}^M{\omega }_i=1\end{array}$(8) where n is number of samples, $y\left(x_j\right)$ denotes the observed values, $h_i^{\left(-i\right)}\left(x_j\right)$ represents the output or prediction of the ith base learner for the jth data point, and ${\omega }^{\ast }$ = (${\omega }_1,{\omega }_2,\dots ,{\omega }_m$) is the set of weights assigned to the base learners.

Model development and input

The ability of four methods (ANN-BP, SVR, CNN, and ensemble stacking model) in simulating the monthly and yearly runoff in the three selected basins was evaluated. To enhance the performance, the IMFs of P, T_min, and T_max were applied as input variables. Different numbers of IMFs were set for the monthly and yearly simulations. The optimal group of variables of the input dataset was selected by applying two feature algorithms: SVM-recursive feature elimination (SVM-RFE) and random forest-Boruta (RF-Boruta) feature selection algorithm for all three basins and models. For a more comprehensive understanding of the SVM-RFE and RF-Boruta feature selection algorithms, please refer to the works of Ahmadpour et al. (2021), Parisouj et al. (2022), Jamei et al. (2023), Farhana et al. (2023), Kursa and Rudnicki (2010), and Maguire et al. (2022). Five groups of input variables were considered for each model: (1) P, T_min, and T_max, and their IMFs selected by SVM-RFE; (2) P, T_min, and T_max, and their IMFs selected by RF-Boruta feature selection; (3) IMFs of the main values obtained using SVM-RFE; (4) IMFs of the main values obtained using RF-Boruta feature selection; and (5) only the three main variables (P, T_min, and T_max). Table 1(a–c) presents the input variables of each method for the Grand, Winnipeg, and Moosonee River basins at monthly and yearly time scales, respectively.

Table 1. (a) Variables selected for each method for the Grand River basin. (b) Variables selected for each method for the Winnipeg River basin. (c) Variables selected for each method for the Moosonee River basin.

(a)
Time scale	Method	Feature selection	Variables
Monthly	SVR-RFE	Main variables and IMFs	P, P′4, P′5, Tmax, Tmax′1, Tmax′2, Tmax′3, Tmax′7, Tmin, Tmin′2, Tmin′3, Tmin′5
		IMFs	P′1, P′2, P′3, P′4, P′5, Tmax′1, Tmax′2, Tmax′3, Tmin′2, Tmin′3
	RF-Boruta	Main variables and IMFs	P, Tmax, Tmax′5, Tmin
		IMFs	P′1, P′2, P′4, P′5, Tmin′2
Yearly	SVR-RFE	Main variables and IMFs	P, P′1, Tmax, Tmax′1, Tmax′2, Tmin′1, Tmin′4
		IMFs	P′1, P′2, P′4, Tmax′1, Tmin′3, Tmin′4
	RF-Boruta	Main variables and IMFs	P, P′1, Tmax, Tmax′1
		IMFs	P′1, Tmax′1
(b)
Monthly	SVR-RFE	Main variables and IMFs	P′2, P′4, P′7, P′8, Tmax, Tmax′2, Tmax′3, Tmax′4, Tmin, Tmin′1, Tmin′6, Tmin′7
		IMFs	P′3, P′4, P′7, P′8, Tmax′2, Tmax′3, Tmax′4, Tmin′3, Tmin′6, Tmin′7
	RF-Boruta	Main variables and IMFs	P′3, P′4, P′5, P′6, P′7, P′8, Tmax′1, Tmax′2, Tmax′4, Tmax′5, Tmax′6, Tmax′7, Tmin, Tmin′1, Tmin′2, Tmin′3, Tmin′4, Tmin′5, Tmin′7
		IMFs	P′3, P′4, P′5, P′6, P′7, P′8, Tmax′1, Tmax′2, Tmax′4, Tmax′5, Tmax′6, Tmax′7, Tmin′1, Tmin′2, Tmin′3, Tmin′4, Tmin′5, Tmin′7
Yearly	SVR-RFE	Main variables and IMFs	P, P′2, P′3, Tmax, Tmax′1, Tmin, Tmin′1
		IMFs	P′2, P′3, Tmax′1, Tmax′2, Tmax′3, Tmin′1
	RF-Boruta	Main variables and IMFs	P, P′2, Tmax, Tmin′1
		IMFs	P′3, Tmax′1
(c)
Monthly	SVR-RFE	Main variables and IMFs	P, P′2, P′3, P′4, P′6, Tmax′2, Tmax′5, Tmin, Tmin′1, Tmin′2, Tmin′5, Tmin′6
		IMFs	P′1, P′2, P′3, P′4, P′6, Tmax′1, Tmax′5, Tmax′6, Tmin′1, Tmin′2, Tmin′5
	RF-Boruta	Main variables and IMFs	P, P′1, Tmax, Tmax′2, Tmin
		IMFs	Tmax′1, Tmax′2, Tmin′1, Tmin′2
Yearly	SVR-RFE	Main variables and IMFs	P, P′1, P′2, Tmax, Tmax′2, Tmax′4, Tmin
		IMFs	P′1, P′2, P′3, P′4, Tmax′2, Tmax′4
	RF-Boruta	Main variables and IMFs	P, P′1
		IMFs	P′1, P′2

Note: P.: precipitation; T.: temperature; min.: minimum; max.: maximum; the symbol ′ represents the IMF.

The average monthly and yearly P, T_min, and T_max and their IMFs were applied as the model input, and the observed streamflow variable was used to evaluate the prediction accuracy of each model. Data from January 1973 to August 2006 were used for model training and those from September 2006 to December 2020 were used for model testing. The normalisation method was used to normalise input variables to improve the model’s learning ability. The mean and standard deviation (μ and σ, respectively) were used in the training phase, defined as in Equation (9). (9) $\text{normalized }x=\frac{x-\mu }{\sigma },$(9) where $x$ represent the original value.

Model parameterisation

To avoid overfitting, the 10-fold cross-validation method was applied to minimise the root mean square error (RMSE) to determine the optimal parameters of the ANN-BP, SVR, CNN, and ensemble models for training each model. The optimal hyperparameters were applied to build the optimal model for simulating the runoff for the training and testing phases. Python 3.8 was used for preparing scripts.

The ANN-BP model was built based on a three feedforward layer configuration with the limited-memory Broyden–Fletcher–Goldfarb­–Shanno solver. The hidden layer consisted of the logistic sigmoid function that delineates the input weights that are transformed into an output from the nodes. The learning rate and max-iteration of the ANN-BP were 0.07 and 10,000, respectively. The optimal hyperparameters of ANN-BP were defined by random generalisation and 10-fold cross-validation for the training model. Three hidden layers were applied to generate the highest runoff accuracy by minimising the RMSE function. Therefore, the network structure was (number of inputs for each state):5:1 for the input, hidden, and output layers.

The SVR model was built based on the RBF kernel, and random generalisation was applied for each parameter to determine the optimal hyperparameters. The three hyperparameters of the RBF kernel were set to have the following ranges to tune the model: γ (30 values from 0.0001 to 1) C (100 values from 2000 to 5000), and ε (21 values from 0.5 to 3).

The CNN architecture included two one-dimensional CNN layers, a dropout layer for regulation, and a pooling layer. CNN layers are frequently created in pairs to aid the model in learning properties from the input data. The dropout layer aims to lower the learning rate of the CNNs to produce a more accurate final model. The pooling layer decreases the size of the learned features to a quarter of their original size, allowing them to be concentrated on the most relevant aspects. Next, the flatten function is used to flatten learned features as a vector and pass them into a fully connected layer. The fully connected layer prevents the learned features and target value to define the learned features before making a prediction. In this study, 512 and 256 filter features with kernel sizes 5 and 2 were used for the first and second CNN layers, respectively. RMSprop was used as the network optimiser. LeakyRelu was used as an activation function. A random generalisation process and 10-fold cross-validation were applied to determine the parameters.

The AdaBoost model involves three key factors: estimators, learning rate, and loss. The estimator and learning rate were set to range as follows: 30 values from 1 to 200 and from 0.0001 to 1, respectively. In the training phase, linear, square, and exponential loss functions were evaluated. In the bagging regressor and extra tree regressor models, the optimal estimator values were identified from a candidate set of thirty variables ranging from 1 to 200 and from 300 to 500, respectively, during training through 10-fold cross-validation. Stacking generation was incorporated into the models to enhance the accuracy of target data. The bagging regressor and extra tree regressor models were selected as estimators, and the AdaBoost model was chosen as the final estimator to simulate the runoff values. To enhance the model performance, 10,000 runs of each model were performed, and the RMSE was used to evaluate the highest accuracy among training periods.

Model evaluation

The coefficient of correlation (R), RMSE, mean absolute relative error (MARE), and Nash–Sutcliffe efficiency (NSE) were applied to assess the runoff accuracy for ANN-BP, SVR, CNN, and ensemble model for the training and testing periods. The definitions of these statistical indices are presented in the following equations: (10) $\begin{array}{rl}\text{RMSE}& =\sqrt{\frac{1}{n}\sum _{i=1}^n{\left(s_i-o_i\right)}^2}\end{array}$(10) (11) $\begin{array}{rl}R& =\frac{\left(\frac{1}{n}\right)\sum _{i=1}^n{\left(o_i-\overline{o}\right)}^2\left(s_i-\overline{s}\right)}{\sqrt{\left(\frac{1}{n}\right)\sum _{i=1}^n{\left(o_i-\overline{o}\right)}^2}\times \sqrt{\left(\frac{1}{n}\right)\sum _{i=1}^n{\left(s_i-\overline{s}\right)}^2}}\end{array}$(11) (12) $\begin{array}{rl}\text{NSE}& =1-\frac{\sum _{i=1}^n{\left(s_i-o_i\right)}^2}{\sum _{i=1}^n{\left(s_i-\overline{o}\right)}^2}\end{array}$(12) (13) $\begin{array}{rl}\text{MARE}& =\frac{1}{n}\sum _{i=1}^n\left|\frac{s_i-o_i}{o_i}\right|\times 100\end{array}$(13) where $o_i$ and $s_i$ refer to the observed and estimated values, respectively; and $\overline{o}$ and $\overline{s}$ are the average observed and estimated values, respectively.

Results and discussion

EMD was performed to simulate the monthly and yearly streamflow based on five scenarios: (1) using precipitation (P), maximum temperature (T_max), and minimum temperature (T_min); (2) using the RF-Boruta algorithm to select appropriate variables among decomposed variables; (3) using the RF-Boruta algorithm to select appropriate variables among decomposed variables including P, T_max, and T_min; (4) using the SVM-RFE algorithm to select appropriate variables among decomposed variables; (5) using the SVM-RFE algorithm to select appropriate variables among decomposed variables including P, T_max, and T_min. All the models were applied for the Grand, Winnipeg, and Moosonee River basins. To clarify the influence of EMD on the simulation performance, simulation models both with and without the EMD dataset were prepared. Additionally, the results for the different feature selection methods were compared to identify the most accurate model for simulating the runoff in Canadian basins.

Decomposing monthly and yearly runoff time-series using EMD

The EMD approach was used to decompose the runoff time-series for the three basins into IMFs at the monthly and yearly scales for the precipitation and maximum and minimum temperatures. Notably, several researchers have applied decomposition methods for analysing the runoff fluctuations of rivers and investigated the causes of cyclical changes in hydrological data and corresponding occurrence mechanisms (Pekárová et al., 2003; Wang et al., 2015; Williams, 1961). Decomposition can help enhance the prediction ability by transforming nonlinear and nonstationary time-series into stationary time-series.

Simulation results using a combination of P, T, and EMD dataset

Monthly time-series

Table 2(a–c) summarises the calibration and validation phase results for Grand, Winnipeg, and Moosonee Rivers. When combined with EMD and SVM-RFE or RF-Boruta for input selection, the ANN-BP and ensemble models are more efficient during calibration. For example, the ANN-BP model with SVM-RFE outperforms the RF-Boruta in the Grand River calibration phase by 48% and 22% in terms of the RMSE and NSE, respectively. In the calibration phase for Winnipeg River, the ensemble model using SVM-RFE outperforms RF-Boruta with a 27% lower RMSE and 26% higher NSE. In most cases, SVM-RFE performs better than RF-Boruta during calibration. However, in the Moosonee River watershed, the situation is reversed. In this case, the SVR-EMD model optimised with RF-Boruta exhibits a 24% higher RMSE and 19% lower NSE compared with SVM-RFE.

Table 2. (a) Monthly model performance indicators during the calibration and validation phases for the Grand River basin (combination of P, T, and EMD dataset). (b) Monthly model performance indicators during the calibration and validation phases for the Winnipeg River basin (combination of P, T, and EMD dataset). (c) Monthly model performance indicators during the calibration and validation phases for the Moosonee River basin (combination of P, T, and EMD dataset).

(a)
		Calibration				Validation
Model	Input Selection	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR-EMD	SVM-RFE	30.95	0.79	0.57	25.96	33.57	0.71	0.49	37.87
	RF-Boruta	34.36	0.7	0.47	38.9	40.09	0.57	0.27	50.14
ANN-BP-EMD	SVM-RFE	11.55	0.97	0.94	12.96	40.22	0.63	0.27	50.87
	RF-Boruta	22.52	0.88	0.77	25.99	44.86	0.52	0.09	56.75
CNN-EMD	SVM-RFE	37.16	0.67	0.38	53.36	39.4	0.67	0.3	39.64
	RF-Boruta	36.2	0.64	0.41	54.44	38.86	0.58	0.31	52.65
Ensemble-EMD	SVM-RFE	21.47	0.9	0.79	26.09	34.62	0.75	0.46	47.22
	RF-Boruta	15.19	0.97	0.9	21.42	37.28	0.62	0.37	50.1
(b)
SVR-EMD	SVM-RFE	285.17	0.77	0.58	29.14	329.62	0.67	0.44	30.53
	RF-Boruta	286.5	0.77	0.57	28.34	348.17	0.62	0.37	37.36
ANN-BP-EMD	SVM-RFE	205.01	0.88	0.78	18.4	356.43	0.63	0.34	35.92
	RF-Boruta	101.18	0.97	0.95	5.41	374.72	0.55	0.27	40.50
CNN-EMD	SVM-RFE	316.15	0.71	0.48	29.78	493.46	0.64	−0.26	63.27
	RF-Boruta	333.89	0.68	0.42	32	431.72	0.6	0.03	52.38
Ensemble-EMD	SVM-RFE	188.36	0.92	0.82	22.8	345.85	0.7	0.38	39.81
	RF-Boruta	258.93	0.81	0.65	25.67	340.3	0.66	0.40	39.84
(c)
SVR-EMD	SVM-RFE	62.78	0.79	0.60	64.39	71.48	0.74	0.51	68.87
	RF-Boruta	50.5	0.87	0.74	47.04	75.41	0.68	0.46	82.43
ANN-BP-EMD	SVM-RFE	60.56	0.79	0.63	78.62	70.56	0.73	0.52	79.04
	RF-Boruta	59.39	0.80	0.64	83.48	67.43	0.76	0.57	88.47
CNN-EMD	SVM-RFE	70.88	0.71	0.49	86.15	79.06	0.65	0.40	91.24
	RF-Boruta	71.88	0.70	0.48	106.92	81.30	0.61	0.37	116.04
Ensemble-EMD	SVM-RFE	49.86	0.87	0.75	47.43	66.35	0.77	0.58	56.83
	RF-Boruta	31.03	0.95	0.90	35.76	72.26	0.71	0.50	88.26

During the validation phase, the SVR model outperforms the other models in terms of the RMSE, NSE, and MARE for the Grand River, although its R value is inferior. The SVR model using SVM-RFE has an 81% higher NSE and a 16% lower RMSE, which confirms (Luo et al., 2022) that SVM-RFE typically selects better variables than RF-Boruta. The ensemble model outperforms other models for Winnipeg and Moosonee Rivers, regardless of its integration with SVM-RFE or RF-Boruta. Compared with RF-Boruta, SVM-RFE and EMD exhibit 1.3% and 8.1% lower RMSEs and 5% and 16% higher NSEs, respectively, for the Winnipeg and Moosonee Rivers. The ANN-BP model outperforms the CNN model in the Moosonee River, with a 14% lower RMSE and 20% higher NSE.

The R values of the ensemble model match those of the SVR model, suggesting a similar correlation between the observed and modelled values. SVR and ensemble models exhibit reduced biases, indicating improved generalizability. The SVR model utilises a structural risk minimisation approach, resulting in a superior solution, whereas the ensemble model adopts multiple weak learners, yielding accurate predictions (Gizaw & Gan, 2016; Htike, 2017; Kumar et al., 2019; Shrestha & Shukla, 2015). Despite the overall superior performance of the ensemble model against SVR, ANN-BP, and CNN models, the SVR model is significantly superior in the Grand River basin.

Figures 3 and 4 show that the simulated and observed streamflows are consistent for both the calibration and validation phases, especially for the Moosonee River. Figure 4 demonstrates that the values predicted by the CNN and ANN-BP models lie in ranges of 50–250 m³s⁻¹ for the Grand River, 550–2,500 m³s⁻¹ for the Winnipeg River, and 50–300 m³s⁻¹ for the Moosonee River. These variances indicate that the ensemble and SVR models are more flexible and accurate in their predictions. In summary, the combination of ensemble and SVR models with SVM-RFE yields superior results, rendering them promising alternatives for streamflow prediction and similar tasks. This observation is consistent with previously reported findings and highlights exciting possibilities for future work.

Figure 3. Best monthly line-graph using a combination of P, T, and EMD dataset.

Figure 4. Scatter plots of monthly streamflow simulation using a combination of P, T, and EMD dataset in the validation period. The different rows present the results of different models.

Yearly time-series

Table 3(a–c) summarises the performance metrics of various models for the Grand, Winnipeg, and Moosonee River basins at the yearly scale. The findings, especially those of the NSE and MARE, emphasise the advantage of SVM-RFE over RF-Boruta in feature selection across these basins. For the Grand River basin, the SVR-EMD model using SVM-RFE exhibits a 30% higher NSE and lower MARE during the calibration phase compared with RF-Boruta. This enhancement is also observed in the validation phase, with an 18% higher NSE. Similar trends can be observed in the cases of the Winnipeg and Moosonee River basins.

Table 3. (a) Yearly model performance indicators during the calibration and validation phases for the Grand River basin (combination of P, T, and EMD dataset). (b) Yearly model performance indicators during the calibration and validation phases for the Winnipeg River basin (combination of P, T, and EMD dataset). (c) Yearly model performance indicators during the calibration and validation phases for the Moosonee River basin (combination of P, T, and EMD dataset).

(a)
	Calibration					Validation
Model	Input Selection	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR-EMD	SVM-RFE	3.4	0.96	0.92	2	10.81	0.91	0.59	12.19
	RF-Boruta	6.61	0.84	0.71	7.93	11.88	0.76	0.5	15.91
ANN-BP-EMD	SVM-RFE	0	1	1	0	12.44	0.76	0.45	15.68
	RF-Boruta	0	1	1	0	24.16	0.26	−1.07	28.3
CNN-EMD	SVM-RFE	10.01	0.81	0.33	12.52	11.29	0.78	0.55	15.96
	RF-Boruta	8.4	0.8	0.53	10.51	18.26	0.47	−0.18	21.4
Ensemble-EMD	SVM-RFE	5.17	0.91	0.82	6.04	12.35	0.88	0.46	13.73
	RF-Boruta	8.64	0.73	0.5	11.4	12.88	0.79	0.41	16.51
(b)
SVR-EMD	SVM-RFE	148.00	0.93	0.75	13.67	191.23	0.87	0.52	16.49
	RF-Boruta	156.90	0.85	0.72	12.29	244.92	0.55	0.21	21.62
ANN-BP-EMD	SVM-RFE	27.14	1.00	0.99	1.26	261.63	0.71	0.10	20.38
	RF-Boruta	92.07	0.95	0.90	6.04	432.03	0.19	−1.45	37.14
CNN-EMD	SVM-RFE	183.33	0.84	0.61	18.76	221.49	0.80	0.36	17.90
	RF-Boruta	221.51	0.75	0.43	21.66	515.22	0.42	−2.49	49.29
Ensemble-EMD	SVM-RFE	36.56	0.99	0.98	2.70	96.22	0.94	0.88	7.73
	RF-Boruta	117.19	0.93	0.84	11.97	230.10	0.62	0.30	18.55
(c)
SVR-EMD	SVM-RFE	6.38	0.95	0.88	6.06	16.71	0.82	0.36	14.54
	RF-Boruta	11.08	0.81	0.65	9.95	14.33	0.73	0.53	11.41
ANN-BP-EMD	SVM-RFE	0.00	1.00	1.00	0.00	19.15	0.53	0.16	15.80
	RF-Boruta	0.00	1.00	1.00	0.00	23.10	0.29	−0.23	18.60
CNN-EMD	SVM-RFE	9.80	0.86	0.72	8.36	13.37	0.79	0.59	11.83
	RF-Boruta	11.82	0.78	0.60	9.06	15.06	0.70	0.48	12.67
Ensemble-EMD	SVM-RFE	0.00	1.00	1.00	0.00	12.40	0.83	0.65	9.96
	RF-Boruta	1.21	1.00	1.00	0.63	13.14	0.79	0.60	11.35

Figures 5 and 6 present line graphs for the calibration and testing periods, and scatter plots for the testing period in the three Canadian basins. The ANN-BP-EMD model paired with SVM-RFE achieves perfect calibration across all studied basins. However, its inconsistent outcome in the validation phase may suggest potential overfitting. In contrast, the SVR-EMD and ensemble-EMD models combined with SVM-RFE yield more consistent results in both calibration and validation phases. Specifically, the ensemble-EMD model paired with SVM-RFE exhibits robust performance across all basins. For example, in the case of the Winnipeg River basin, it achieves an NSE of 0.98 and a MARE of 2.70 during calibration, significantly surpassing both CNN-EMD and SVR-EMD. This high performance is retained in the validation phase, with an NSE of 0.88. Similar trends are observed for the Moosonee River basin. Overall, in the annual simulations involving scenario 5, SVM-RFE consistently outperforms RF-Boruta in feature selection across all basins, and the ensemble-EMD model outperforms the other models.

Figure 5. Yearly line-graph using a combination of P, T, and EMD dataset.

Figure 6. Scatter plots of yearly streamflow simulation using a combination of P, T, and EMD dataset in the validation period.

To clarify the influence of the combined P, T, and EMD on the model efficiency, simulation models using only the EMD technique were evaluated, as described in the following section.

Simulation results using the EMD dataset

Monthly time-series

Valuable insights could be derived from the analysis of models across the three river basins. In the case of the Grand River basin (Table 4(a)), the ANN-BP-EMD model with SVM-RFE selection achieves the highest NSE of 0.93 during calibration, but this value sharply deteriorates to −0.63 in the validation phase, indicating potential overfitting. In contrast, the ensemble-EMD model with RF-Boruta consistently performs well in both stages with an NSE of 0.49 and 0.35 in the calibration and validation phases, respectively.

Table 4. (a) Monthly model performance indicators during the calibration and validation phases for the Grand River basin (using only the EMD dataset). (b) Monthly model performance indicators during the calibration and validation phases for the Winnipeg River basin (using only the EMD dataset). (c) Monthly model performance indicators during the calibration and validation phases for the Moosonee River basin (using only the EMD dataset).

(a)
		Calibration				Validation
Model	Input Selection	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR-EMD	SVM-RFE	16.13	0.95	0.88	3.82	42.97	0.43	0.16	56.08
	RF-Boruta	38.96	0.59	0.32	41.56	37.76	0.62	0.35	43.80
ANN-BP-EMD	SVM-RFE	12.05	0.97	0.93	13.59	60.02	0.22	−0.63	80.10
	RF-Boruta	22.46	0.88	0.77	31.74	52.57	0.38	−0.25	73.35
CNN-EMD	SVM-RFE	39.34	0.56	0.3	61.72	42.94	0.42	0.16	61.89
	RF-Boruta	42.43	0.45	0.19	64.98	42.34	0.45	0.19	60.49
Ensemble-EMD	SVM-RFE	23.5	0.92	0.75	18.04	39.07	0.56	0.31	54.81
	RF-Boruta	33.61	0.73	0.49	40.81	37.72	0.60	0.35	56.76
(b)
SVR-EMD	SVM-RFE	284.63	0.76	0.58	32.20	389.54	0.56	0.21	45.50
	RF-Boruta	290.59	0.77	0.56	28.15	377.18	0.56	0.26	42.59
ANN-BP-EMD	SVM-RFE	238.26	0.84	0.71	26.12	400.18	0.51	0.17	39.47
	RF-Boruta	66.97	0.99	0.98	4.49	635.00	0.40	−1.09	58.28
CNN-EMD	SVM-RFE	324.91	0.69	0.45	32.00	381.12	0.55	0.25	42.65
	RF-Boruta	325.32	0.69	0.45	30.95	479.92	0.47	−0.20	42.69
Ensemble-EMD	SVM-RFE	215.90	0.88	0.76	23.18	390.41	0.59	0.21	48.69
	RF-Boruta	230.15	0.86	0.73	25.09	347.61	0.65	0.37	40.82
(C)
SVR-EMD	SVM-RFE	75.58	0.71	0.42	52.35	92.05	0.51	0.19	94.85
	RF-Boruta	68.25	0.75	0.53	53.95	88.16	0.53	0.26	101.44
ANN-BP-EMD	SVM-RFE	89.53	0.46	0.19	92.30	100.95	0.46	0.03	97.01
	RF-Boruta	52.14	0.85	0.72	69.92	100.95	0.46	0.03	97.01
CNN-EMD	SVM-RFE	87.58	0.61	0.22	78.45	97.43	0.49	0.09	86.46
	RF-Boruta	78.72	0.62	0.37	121.53	87.54	0.53	0.27	122.22
Ensemble-EMD	SVM-RFE	57.50	0.86	0.66	41.53	89.01	0.51	0.24	120.60
	RF-Boruta	57.50	0.86	0.66	41.53	89.01	0.51	0.24	120.60

In the case of the Winnipeg River basin (Table 4(b)), the ANN-BP-EMD model with RF-Boruta demonstrates satisfactory calibration performance with an NSE of 0.98. However, its performance deteriorates during the validation phase. In comparison, the ensemble-EMD model with RF-Boruta consistently exhibits a strong performance, with an NSE of 0.73 and 0.37 during calibration and validation, respectively. Notably, the ensemble model outperforms all other models, and the worst performance corresponds to the ANN-BP model using RF-Boruta.

For the Moosonee River basin (Table 4(c)), the ANN-EMD model with RF-Boruta for feature selection demonstrates strong performance during calibration (NSE = 0.72), but its performance deteriorates in the validation phase (NSE = 0.03). In contrast, the ensemble-EMD model exhibits a high performance in both phases (NSE = 0.66 and 0.24 in calibration and validation, respectively), regardless of the feature selection method used (SVM-RFE or RF-Boruta).

In general, the inclusion of main variables such as precipitation and temperature in the EMD dataset improves model performance, except for the SVR model with RF-Boruta in the Grand River basin, where the improvement is minimal. This observation is consistent with that reported by Zhu and Pierskalla (2016), which emphasises the potential bias introduced by a large number of input features. Figures 7 and 8 present the performance metrics of models using the EMD dataset with RF-Boruta across the three Canadian basins. The model ranking in calibration and validation is ANN-BP > ensemble > SVR > CNN and ensemble > SVR > CNN > ANN-BP, respectively.

Figure 7. Best monthly line-graph using the EMD dataset.

Figure 8. Scatter plots of monthly streamflow simulation using the EMD dataset in the validation period.

Yearly time-series

Table 5(a–c) presents the RMSE, R, NSE, and MARE metrics for the model during the training and validation periods. Because of the inherent structure of neural network models, their performance during the validation phase is inferior to that during calibration. In the case of the Grand River basin, the SVR model using RF-Boruta and the ensemble model using SVM-RFE exhibit excellent performances, with the highest R values of 0.85 and 0.82, respectively, during validation. The corresponding RMSE values are 10.18 and 10.06 m³s⁻¹, and the NSE values are 0.63 and 0.64. Table 5(a) reveals that although the models combined with RF-Boruta generally outperform those combined with SVM-RFE, the ensemble model that uses SVM-RFE for feature selection demonstrates the best performance among all models.

Table 5. (a) Yearly model performance indicators during the calibration and validation phases for the Grand River basin (using only the EMD dataset). (b) Yearly model performance indicators during the calibration and validation phases for the Winnipeg River basin (using only the EMD dataset). (c) Yearly model performance indicators during the calibration and validation phases in the Moosonee River basin (using only the EMD dataset).

(a)
		Calibration				Validation
Model	Input Selection	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR-EMD	SVM-RFE	7.43	0.8	0.63	9.01	15.2	0.88	0.18	16.95
	RF-Boruta	8.8	0.71	0.49	10.95	10.18	0.85	0.63	11.66
ANN-BP-EMD	SVM-RFE	0	1	1	0	18.41	0.66	−0.2	18.57
	RF-Boruta	0	1	1	0	18.03	0.53	−0.15	20.51
CNN-EMD	SVM-RFE	19.69	0.8	−1.58	28.97	21.74	0.7	−0.67	29.51
	RF-Boruta	10.82	0.57	0.22	14.98	13	0.65	0.4	18.63
Ensemble-EMD	SVM-RFE	7.49	0.8	0.63	9.7	10.06	0.82	0.64	13.91
	RF-Boruta	8.66	0.78	0.50	11.78	12.81	0.66	0.42	18.12
(b)
SVR-EMD	SVM-RFE	139.57	0.93	0.78	12.23	191.39	0.85	0.52	16.90
	RF-Boruta	155.03	0.86	0.72	9.87	211.75	0.73	0.41	17.83
ANN-BP-EMD	SVM-RFE	0.01	1.00	1.00	0.00	274.01	0.68	0.01	22.05
	RF-Boruta	45.85	0.99	0.98	2.48	271.54	0.45	0.03	22.69
CNN-EMD	SVM-RFE	130.80	0.90	0.80	12.25	235.64	0.72	0.27	18.63
	RF-Boruta	170.93	0.82	0.66	17.33	231.12	0.60	0.30	22.68
Ensemble-EMD	SVM-RFE	42.08	0.99	0.98	3.48	213.66	0.76	0.40	15.47
	RF-Boruta	187.61	0.79	0.59	14.95	185.94	0.77	0.55	15.37
(c)
SVR-EMD	SVM-RFE	10.21	0.85	0.70	7.34	16.67	0.71	0.36	13.97
	RF-Boruta	13.76	0.71	0.45	11.65	16.71	0.67	0.36	14.85
ANN-BP-EMD	SVM-RFE	9.21	0.87	0.75	7.46	18.21	0.53	0.24	15.35
	RF-Boruta	9.49	0.86	0.74	8.37	22.52	0.47	−0.17	18.85
CNN-EMD	SVM-RFE	16.10	0.73	0.25	11.62	34.71	0.68	−1.77	31.36
	RF-Boruta	15.68	0.63	0.29	11.15	25.28	0.55	−0.47	21.03
Ensemble-EMD	SVM-RFE	12.45	0.76	0.55	11.83	14.71	0.72	0.50	12.50
	RF-Boruta	8.18	0.94	0.81	4.72	16.11	0.71	0.40	11.74

Table 5(b) shows that for the Winnipeg River basin, the SVR model based on SVM-RFE outperforms the other models, with the highest R values. Moreover, the NSE is 0.52 and RMSE is 191.39 m³s⁻¹. For the Moosonee River basin (Table 5(c)), the ensemble model exhibits excellent performance during the validation period, achieving the highest R values of 0.72 and 0.71 when coupled with the SVM-RFE and RF-Boruta selection methods, respectively. The corresponding NSE values are 0.50 and 0.40, and the RMSE values are 14.71 and 16.11 m³s⁻¹. These metrics are highly improved compared with the other models, and in particular, the RMSE is 5% lower.

Figures 9 and 10 show the plots of the observed and simulated yearly streamflows. The simulation results obtained using the SVR and ensemble models are closer to the observed streamflows than those of the CNN and ANN-BP, which tend to overestimate the yearly streamflow values. The overall performance ranking of the models can be summarised as follows: SVR is comparable or superior to the ensemble model, which outperforms CNN, which in turn is superior to ANN-BP.

Figure 9. Best yearly line-graph using the EMD dataset.

Figure 10. Scatter plots of yearly streamflow simulation using the EMD dataset in the validation period.

Simulation results using P and T

Monthly time-series

Table 6(a–c) summarises the performance metrics of the models in the training and testing phases. In both phases, the ensemble model consistently outperforms the other models, with the highest R, RMSE, NSE, and MARE values across the Grand, Winnipeg, and Moosonee Rivers. For the Grand River basin, the ensemble model shows an RMSE of 36.27 m³s⁻¹ in the validation phase, approximately 15% lower than that of the SVR (42.44), indicating reduced prediction errors. Moreover, its NSE is nearly two times higher, indicating an increase of almost 100%, and its MARE is approximately 8% lower than that of the SVR model. Table 6(b) reveals that the trends for the Winnipeg River are similar: The RMSE is 329.46 m³s⁻¹, which is approximately 40% lower than that of the SVR (543.86 m³s⁻¹). The NSE of the ensemble model (0.44) is significantly higher than that of the SVR (−0.53), and its MARE is nearly 53% lower (72.29). In Table 6(c), similar observations can be made for the Moosonee River. The ensemble model exhibits a 14% lower RMSE and 56% higher NSE compared with the SVR model.

Table 6. (a) Monthly model performance indicators during the calibration and validation phases for the Grand River basin (using P and T). (b) Monthly model performance indicators during the calibration and validation phases for the Winnipeg River basin (using P and T). (c) Monthly model performance indicators during the calibration and validation phases for the Moosonee River basin (using P and T).

(a)
	Calibration				Validation
Model	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR	32.23	0.74	0.53	33.72	42.44	0.51	0.18	46.73
ANN-BP	22.02	0.88	0.78	29.56	54.76	0.31	−0.36	61.01
CNN	39.03	0.65	0.32	40.42	43.23	0.62	0.15	51.4
Ensemble	20.17	0.91	0.82	19.6	36.27	0.65	0.4	42.94
(b)
SVR	323.72	0.69	0.46	31.34	543.86	0.60	−0.53	72.29
ANN-BP	238.26	0.84	0.71	26.12	485.98	0.22	−0.23	48.71
CNN	349.49	0.61	0.37	40.46	490.70	0.62	−0.25	64.34
Ensemble	286.20	0.77	0.58	34.19	329.46	0.66	0.44	34.03
(c)
SVR	68.08	0.76	0.53	53.42	84.31	0.59	0.32	122.38
ANN-BP	46.81	0.90	0.78	109.95	111.97	0.50	−0.20	183.62
CNN	78.98	0.67	0.37	84.24	79.03	0.68	0.40	80.03
Ensemble	42.53	0.91	0.82	34.29	72.42	0.71	0.50	85.59

Figure 11 compares the observed and predicted streamflows by the SVR and ensemble models, highlighting the superior performance of these two models. As shown in Figure 12, during the validation phase, the ensemble model outperforms the CNN, SVR, and ANN-BP across all river basins. In the validation phase (Figure 12), the ensemble model outperforms the CNN, SVR, and ANN-BP across all river basins. Therefore, the ensemble model is noted to be the most reliable for predicting streamflow.

Figure 11. Best monthly line-graph using P and T.

Figure 12. Scatter plots of monthly streamflow simulation using P and T in the validation period.

Yearly time-series

Table 7(a–c) indicate that during the calibration phase, the ANN-BP model exhibits excellent performance across all basins. However, its performance deteriorates in the validation phase, particularly for the Grand River basin (Table 7(a)). This deterioration is evidenced by the R value, which decreases from 1 to 0.55, corresponding to a 45% reduction, and its RMSE, which increases from 0 to 17.3 m³s⁻¹. Similar trends are observed for the Winnipeg River basin (Table 7(b)): The RMSE of ANN-BP increases by nearly 195%, from 101.17 to 298.07 m³s⁻¹. Similarly, in the case of the Moosonee River basin (Table 7(c)), the model performance deteriorates, with the RMSE increasing from 0 to 24.18 m³s⁻¹ and the R value decreasing from 1 to 0.55, corresponding to a 45% reduction. These findings highlight the inconsistency of the model performance in the calibration and validation phases.

Table 7. (a) Yearly model performance indicators during the calibration and validation phases for the Grand River basin (using P and T). (b) Yearly model performance indicators during the calibration and validation phases for the Winnipeg River basin (using P and T). (c). Yearly model performance indicators during the calibration and validation phases for the Moosonee River basin (using P and T).

(a)
	Calibration				Validation
Model	RMSE (m^3s^−1)	R	NSE	MARE	RMSE (m^3s^−1)	R	NSE	MARE
SVR	7.37	0.81	0.64	8.62	13.62	0.88	0.34	14.52
ANN-BP	0	1	1	0	17.3	0.55	−0.06	20.2
CNN	8.8	0.73	0.48	11.2	12.33	0.71	0.46	16.82
Ensemble	7.76	0.79	0.6	9.03	11.43	0.77	0.54	14.03
(b)
SVR	155.75	0.85	0.72	11.65	196.04	0.72	0.50	19.60
ANN-BP	101.17	0.94	0.88	4.67	298.07	0.53	−0.17	20.35
CNN	212.12	0.70	0.48	23.02	232.46	0.58	0.29	23.06
Ensemble	142.28	0.91	0.77	9.78	187.45	0.79	0.54	17.48
(c)
SVR	10.82	0.82	0.66	8.55	15.19	0.72	0.47	13.09
ANN-BP	0.00	1.00	1.00	0.00	24.18	0.55	−0.35	23.30
CNN	14.33	0.73	0.41	11.70	21.99	0.67	−0.11	17.91
Ensemble	11.95	0.77	0.59	11.94	15.79	0.68	0.43	13.21

In the validation phase, the SVR and ensemble models outperform the ANN-BP across all basins. In the case of the Grand River basin, the ensemble records a 30% lower RMSE (11.43 m³s⁻¹) and a 40% higher R value (0.77). Similarly, in the case of the Winnipeg River basin, the validation RMSE of the ensemble model is 37% lower (187.45 m³s⁻¹) and its R value is 49% higher (0.79). Finally, in the case of the Moosonee River basin, the SVR model exhibits superior performance with a 37% lower RMSE (15.19 m³s⁻¹) and 31% higher R value (0.72). The SVR model, which has an RMSE only 4% lower than that of the ensemble model, can be considered slightly inferior to the ensemble model.

Figures 13 and 14 show the results of simulations using the three main variables as the predictor variables. Figure 13 shows that the performance of most models in the case of the Grand River basin is higher than those for the Winnipeg and Moosonee River basins. Figure 14 shows the scatter plots for each model for the three Canadian basins in the validation period. For all three basins, all four models, especially ANN-BP, overestimate or underestimate the streamflow variable.

Figure 13. Yearly line-graph using P and T.

Figure 14. Scatter plots of yearly streamflow simulation using P and T in the validation period.

Conclusions

This paper proposes a hybrid model based on EMD and four ML models: ensemble, SVR, CNN, and ANN-BP models, to simulate the monthly and yearly runoff time-series and increase the simulation accuracy of long-term runoff. ML models based on the original variables (monthly and yearly time-series of the precipitation and maximum and minimum temperatures) were developed as a comparative standard. Monthly and yearly runoff data from the Grand, Winnipeg, and Moosonee Rivers in Canada were used, and four statistical metrics (RMSE, MARE, R, and NSE) were adopted to assess the model performances. The results demonstrated that the EMD increases the simulation precision, and the proposed EMD-ML models are superior to the standalone ML models in monthly and annual runoff time-series modelling. The proposed hybrid approach can be applied in future research for simulating monthly and annual runoff. The advantages of the proposed technique can be summarised as follows: First, despite the simplicity of the EMD, it can offer valuable insights into the characteristics of the monthly and yearly runoff time-series. Second, the monthly scale is associated with a lower accuracy than the yearly scale. Third, the ensemble model outperforms the SVR, CNN, and ANN-BP for the Grand, Winnipeg, and Moosonee River basins. Finally, the proposed models do not necessitate intricate decision-making regarding the explicit form for each instance. Overall, a hybrid simulation model with EMD can yield precise and consistent simulation results, and it is thus a valuable tool for studies that focus on hydrological time-series simulations to address various problems associated with reservoir management.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the Korea Environmental Industry and Technology Institute (KEITI) through the Project for Developing Innovative Drinking Water and Wastewater Technologies, funded by Ministry of Environment (MOE) (No.2020002700015); the Ministry of Science and ICT through the National Research Foundation of Korea (No. NRF-2022R1A4A3032838 and RS-2023-00222333); and the Chung-Ang University Research Grants in 2021.