An improved framework to predict river flow time series data

Due to non-stationary and noise characteristics of river flow time series data, some pre-processing methods are adopted to address the multi-scale and noise complexity. In this paper, we proposed an improved framework comprising Complete Ensemble Empirical Mode Decomposition with Adaptive Noise-Empirical Bayesian Threshold (CEEMDAN-EBT). The CEEMDAN-EBT is employed to decompose non-stationary river flow time series data into Intrinsic Mode Functions (IMFs). The derived IMFs are divided into two parts; noise-dominant IMFs and noise-free IMFs. Firstly, the noise-dominant IMFs are denoised using empirical Bayesian threshold to integrate the noises and sparsities of IMFs. Secondly, the denoised IMF’s and noise free IMF’s are further used as inputs in data-driven and simple stochastic models respectively to predict the river flow time series data. Finally, the predicted IMF’s are aggregated to get the final prediction. The proposed framework is illustrated by using four rivers of the Indus Basin System. The prediction performance is compared with Mean Square Error, Mean Absolute Error (MAE) and Mean Absolute Percentage Error (MAPE). Our proposed method, CEEMDAN-EBT-MM, produced the smallest MAPE for all four case studies as compared with other methods. This suggests that our proposed hybrid model can be used as an efficient tool for providing the reliable prediction of non-stationary and noisy time series data to policymakers such as for planning power generation and water resource management.


INTRODUCTION
The economic development of any country is directly linked to the proper management of their water resources operations that can minimize the effects of various natural disasters such as floods and droughts. Therefore, river flow time series prediction is an imperative task, which plays a significant role for effective and appropriate water resource planning and management, early flood warning, irrigation, and hydropower generation (Yaseen et al., 2018a;Yaseen et al., 2018b). Several algorithms have been used for prediction and estimation for river flow time series data (Aichouri et al., 2015;Shathir & Saleh, 2016;AlMasudi, 2018;Gjika, Aurora & Arbesa, 2019). The Box-Jenkins methodology technique (Box & Jenkins, 1970) is commonly used in the literature (Shathir & Saleh, 2016) because it can be used in a wide class of models i.e., Autoregressive (AR), Moving Averages (MA), Autoregressive Integrated Moving Averages Model (ARIMA), Seasonal ARIMA. But the problem with such statistical models is that they only considered the stationary and linear behavior of the process (Di, Yang & Wang, 2014). However, river flow time series data is non-linear, non-stationary, multi-scale and noise-corrupted (Di, Yang & Wang, 2014) as stochastic nature of several factors (e.g., rainfall, evaporation, and temperature). This complex non-stationary, multi-scale and noise-corrupted characteristics make the prediction a challenging task. During the past decades, this issue has been addressed by developing some data-driven models i.e., Artificial Neural Network (ANN), model tree, Support Vector Machine (SVM), Adaptive Inference-Based Neural Network (AIBNN), Extreme Learning Machine (ELM) (Ali et al., 2018). Yaseen et al. (2018a) and Yaseen et al. (2018b) demonstrated the viability of ELM and enhanced version of ELM method to forecast river flow data as compared to other statistical models. However, the data-driven models ignore the time-varying and noise characteristics of river flow processes which deprives the researcher to efficiently predict such data. To address the drawback of such data-driven models, several hybrid models are introduced to extract the time-varying information and reduce noises which ultimately increase the prediction accuracy of data-driven models (Toth, Brath & Montanari, 2000;Di, Yang & Wang, 2014;Su et al., 2016;Kang et al., 2017;Hadi & Tombul, 2018;Yaseen et al., 2018a;Yaseen et al., 2018b). The hybrid model uses data pre-processing methods such as Singular Spectrum Analysis (SSA), Wavelet Analysis (WA), Empirical Mode Decomposition (EMD) and Empirical-EMD (EEMD) with data-driven models also called intelligence models. An advantage of such data decomposition methods is that they used not only for decomposing the data into time-frequency components but also used to separate noises from data. Wu & Chau (2011) coupled data pre-processing techniques i.e., SSA with ANN to accurately model the rainfall and river flow data. Azadeh et al. (2011) demonstrated the ability to use the data pre-processing method to enhance the precision of data-driven models. They used various data processing techniques and reported that the processed non-linear data is efficiently forecasted with simple statistical and intelligent models. Gjika, Aurora & Arbesa (2019) found that proper consideration of the volatility nature of non-linear data through data pre-processing methods could improve the prediction quality. Santos et al. (2019) proposed a model by coupling WA and ANN techniques. They used WA to transform the daily flow time series data to enhance the precision of the ANN model. The data pre-processing methods used to decompose the non-stationary and non-linear data into physical modes of time-frequency components (Han & Liu, 2009;Azadeh et al., 2011;Wang et al., 2015;Peng et al., 2017). The derived time-frequency components, also called multi-scale components (Nazir et al., 2019), are predicted through a mixture of data-driven models accurately (Azadeh et al., 2011). Further, the time-frequency components, are filtered out by using appropriate thresholds. The reason for using a filter is to preserve significant features of original time series data while removing noises or sparsity from multi-scale components. A growing body of research on denoising found that the extracted high and low multi-scale components in different fields can be found by using linear and non-linear thresholds. Moreover, many other methods of thresholds including the Stein Unbiased Risk Estimator (SURE) (Candes, Sing-Long & Trzasko, 2012;Hansen, 2017), the fixed and soft threshold (Di, Yang & Wang, 2014;Nazir et al., 2019), and minimax algorithms (Hansen, 2017) are also used to remove noises and to preserve important information from complex data. Peng et al. (2017) developed a novel hybrid model by employing an empirical wavelet transform estimator to remove the redundant noises from river flow data. Further, the denoised data are predicted through Particular Swarm Optimization based-ANN model (PSO-ANN). They demonstrated the efficiency of their proposed model over simple PSO-ANN model without denoising. Di, Yang & Wang (2014) proposed a hybrid model by considering two data pre-processing methods i.e., EMD and WA with a soft and hard threshold to find the denoised time-varying information to decrease the complexity of hydrological series. Holzfuss & Kadtke (1993) suggested that the noise reduction methods coupled with the radial basis functions may enhance the quality of traditional statistical and data-driven models.
However, WA-, EMD-, and EEMD-based noise removal techniques have their own drawbacks in extracting the optimal multi-scale components. WA-based hard and soft threshold first comprised of the choice of mother wavelet, which is subjectively selected among many wavelet basis functions. The subjective selection of mother wavelet may also cause errors which decrease the performance of hard and SURE thresholds. Moreover, the EMD is a purely data-driven technique, used for pre-processing data, which is effected through its own mathematical property of mode mixing problem resulting in spurious time-frequency information. To overcome the drawback of EMD, the EEMD is introduced which is an improved version of EMD to solve the mode mixing problem. The EEMD added Gaussian white noise successively to solve noise-assisted (Hadi & Tombul, 2018;Yaseen et al., 2018a;Yaseen et al., 2018b). Although an improved EEMD has proved useful for denoising hydrological time series data (Jiao, Guo & Ding, 2016), it also has some drawbacks that may deprive the researchers in extracting accurate multi-scale components i.e., Intrinsic Mode Function (IMF) by simple averaging them. However, to cope with the simple averaging problem of EEMD, Complete Ensemble Empirical Mode Decomposition with Adaptive Noise (CEEMDAN) have been proposed by Torres et al. (2011). The application of CEEMDAN is successfully applied to derive the timescale components of hydrological time series data (Antico, Schlotthauer & Torres, 2014). Johnstone & Silverman (2005) have considered an empirical Bayesian approach to threshold the multi-scale components derived from WA. They argued that the empirical Bayesian method efficiently models the sparsity and noises of complex data by considering multiple priors for each level. One would hope such methods which possibly estimate thresholds that stably reflect the noises from sparse multi-scale to enhance the accuracy and reliability of prediction performance of complex river flow data.
In this study, we aimed to improve denoising stage of the hybrid model by the novel framework i.e., CEEMDAN-Empirical Bayesian Threshold estimator (EBT) estimator to optimally reduce noises from IMFs which are further used as inputs to get a precise prediction of river flow data which plays a decisive role in the accurate prediction. The principal motivation of choosing EBT is that it is purely a data-based method which deals different level of noises efficiently because of some high multi-scale components derived from CEEMDAN relatively sparse than the lower time-scale components.

PROPOSED METHOD
The proposed method is comprised of four stages such as decomposed, denoised using novel threshold, prediction, and aggregation. The CEEMDAN method is used as a decomposition tool to handle the non-linear and non-stationary data by extracting IMFs. The extracted IMFs are further divided into two parts; one part is comprised of noisiest IMFs which contains errors and sparsity and the second part is noise free IMFs. The noisiest IMF's coupled with novel Empirical Bayesian Threshold (EBT) estimator to get a ride from noises and sparsity to denoise IMF's. Then the denoised IMF's are predicted through the complex data-driven model and noise free IMF's are predicted through a simple stochastic model. Finally, all the predicted IMF's are aggregated to get the final prediction. The proposed framework of deriving multi-scale IMF's through CEEMDAN coupled with optimal denoising method i.e., EBT plays a vital role in the accurate prediction of river flow data. To our convenience, the proposed strategy is labialized as; CEEMDAN (decomposed), EBT (denoised) and Multi Models MM (data-driven and stochastic model) i.e., CEEMDAN-EBT-MM and the proposed scheme is illustrated in Fig. 1.

Decomposition stage
EMD: To decompose non-linear and non-stationary data, the EMD method introduced by Huang et al. (1998) which decomposed data into IMF's by satisfying two conditions as follows: (a) The number of zero crossings and extreme, from complete data, must be equal or differ to one at most; (b) At all levels, mean value of envelope must be zero.
The complete EMD procedure is defined as: 1. Find all local maxima and minima from data y (t ), ( t = 1,2,..,N ). Use cubic splines interpolation to find an upper envelope from maxima e max(t ) and lower envelope from minima e min(t ) . 2. Calculate the average of upper and lower envelope m(t ) = (e max(t ) + e min(t ) )/2. Take the deviations between original time series data and calculated envelope mean as: (1) 3. Match the requirements of g (t )which is defined in (a) and (b) as an IMF, if the conditions are satisfied then mark this g (t ) as ith IMF. 4. In the next step replace original time series y (t ) by r (t ) = y (t ) − g (t ), if g (t ) does not meet the (a) and (b) then just replace y (t ) with g (t ). 5. The process of (1-4) is repeated until no IMF is being extracted from the residue r (t ). In the end, original time series data will be written as the summation of all extracted IMF's and residue as: where g i (t ) is the ith IMF and r(t ) is the trend of the signal. However, EMD techniques have two drawbacks. The one is the endpoint problem that is the extreme values of two ends of the series can't be determined properly which distort the IMF's and the other is mode mixing aliasing in which same IMF contains even more than two frequencies.
EEMD: To resolve the mode mixing problem of EMD, (Wu & Huang, 2004) proposed EEMD. The EEMD decompose the non-linear signals into IMF's as follows: a) Add a white Gaussian noise series to the original data set as follows: (3) b) Decompose the new y 2 (t ) with EMD and obtains the IMFs. c) Repeat step (a) and (b) mth time i.e., (j = 1,2,. . . ,m) with different white noises to get IMF's from new series. d) Find the ensemble means of all IMFs obtained as mthensemble time as follows: where k = 1,2,..,K is kth IMF.
where k = 1,2,..,K is kth IMF. CEEMDAN : Although the mode mixing problem is alleviated with EEMD, taking the simple averages of IMF's without considering the independent addition of white noises could not completely remove noises from IMF's. To tackle the simple averaging problem of EEMD, the CEEMDAN function is introduced by Torres et al. (2011). Here in our study, CEEMDAN is used to decompose the river flow data which is briefly described as follows: 1) The CEEMDAN add Gaussian noises like EEMD in signals as follows: (5) where w 0 is the amplitude of the added white noises and (j = 1,2,. . . ,m). Find the first IMF using simple EEMD defined as: 2) Compute the deviation of original signals from the first IMF as: Decompose r 1 (t ) + w 0 n j (t ) to get first IMF and find the second IMF as: 4) Repeat the (2-3) until stoppage criteria are met and the residual contains not more than two extremes. Finally, the residual is defined as: However, the selection of a number of ensemble and amplitude of white noise is still an open challenge but here in our study we used a number of ensemble members as 100 and standard deviation of white noise is settled as 0.2 according to Di, Yang & Wang (2014).
Identification of noisiest IMFs: After deriving the IMF's, next step is to screen out the noise only IMF's and noise free IMFs (Wei et al., 2016). Two types of IMF's are derived from CEEMDAN/EEMD; the IMF's contains high frequency which is corrupted with sparsity and noises (Wei et al., 2016) and the second part of IMF's comprised on low frequencies which are free from noises (Wei et al., 2016). To get numerical validation, cross-correlation is calculated between all extracted IMF's and original river flow data. The low cross-correlation implies that the high-frequency IMF's are overwhelmed with noises.
Then, the noise only IMF's are further denoised through the appropriate threshold to get the important features from them and to get rid of noises.

Denoising stage
After decomposing the non-linear and non-stationary data into IMF components, an appropriate estimator is chosen to remove noises and sparsities from extracted IMFs. The reason for selecting an appropriate estimator is to find an optimal value of threshold as the highest threshold value would lead to biases, whilst the lowest threshold value would increase the noise variance. The IMFs which are extracted are mostly empirical Bayesian threshold estimator is adapted to denoised the noisiest IMFs. Later, the existing soft and hard threshold and improved threshold functions are used for comparison purposes. Details of all estimators are given below: CEEMDAN/EEMD-based Empirical Bayesian Threshold: To estimate the sparseness and noises from decomposed IMF's (Wei et al., 2016), an Empirical Bayes Threshold (EBT) inspired from wavelet denoising (Johnstone & Silverman, 2005;Jansen & Bultheel, 2014) is used. For the successful application of EBT, first, all the data is scaled transformed to efficiently select the prior distribution for noises and sparsities. After the scaled transformation data follows N (θ i ,1). Then a mixture of priors for θ i are considered as follows: where δ 0 (θ ) is a zero part of scaled data and γ (θ ) is a density of non-zero part. The density of prior should be chosen in such a way that it must belong to a family of distributions whose tails decays at polynomial rates. The parameters and weights of a mixture of prior distributions are estimated through maximum likelihood approach. The reason for using this method to estimate unknowns is that it estimates weights and parameters to be proportional to the likelihood function evaluated at the estimators based on data (Hossain, Kozubowski & Podgórski, 2018). Finally, the posterior medianθ i (imf ,w) is calculated from a mixture of prior distribution is given as follows: which is used as a threshold rule forμ given data (Johnstone & Silverman, 2005). In general, an estimation rule comprised on η (imf,t ) defined for all t > 0, is a thresholding rule if and only if for all t > 0, η imf ,t is an antisymmetric and increasing function of data and η (imf,t )=0 if and only if imf ≤ t where t is defined as a median value which is calculated through the Eq. (11).
Other traditional threshold estimators: To compare the EBT estimator with other nonlinear estimators to suppress the noises and sparsities from noisiest IMFs, Soft Threshold (ST), Hard Threshold (HT) and Improved Threshold Function (ITF) are used which are most widely used in literature (Jeng et al., 2007;Candes, Sing-Long & Trzasko, 2012;Jansen & Bultheel, 2014) listed as follows, respectively; and where, T k is the threshold value calculated as T k = a √ 2E k ln(N ), where k = 1,2,. . . ..,K and a is constant takes the values between 0.4 to 1.4 with a step of 0.1 and E k = median(|IMF t ,k |,t = 1,2,..,N )/0.6745 is median deviation of kth IMF.

Prediction and aggregation
In the prediction stage, the decomposed-denoised IMFs of river flow data are predicted through some data-driven and statistical models. Specifically, the denoised IMFs are predicted through a data-driven model whenever noise-free IMFs and the residual are predicted through a simple stochastic model. To train the model, 70%, 80%, and 90% data are used, and the remaining 30%, 20% and 10% data are used to test the accuracy of the models. The selected models are briefly described as follows: The denoised-IMF prediction with the neural network: A data-driven technique has been proved a powerful tool to model complex non-linear data (Campisi-Pinto, Adamowski & Oron, 2012). The Multi-Layer Perceptron (MLP) which is the most popular sub-model of NN (Talaee, 2014;Ali et al., 2017) consists of three layers of nodes used here for the prediction of denoised IMFs of river flow data. The complete layout of MLP is given in Talaee (2014). The structure of MLP comprised of NN structure. Three nodes are included in MLP as an input layer, hidden layer, and an output layer. First, the single output is calculated by using linear combinations of inputs which are further transferred to some non-linear activation functions mathematically defined as by where w i are the weights of inputs, x i are the inputs, b is the biased value for each layer and ϕ is a non-linear activation function which supplies output to the next layer. The mostly used activation function i.e., logistic function is used as activation function defined as follows: 1 (1 + e − n i w i x i +b ) For the optimization of neurons, the supervised learning algorithms called backpropagation and forward-propagation can be used.
The noise-free IMF prediction with ARIMA model: To predict the noise-free IMF's and the residual, an autoregressive moving average model (ARMA) is selected which is described as follows: where, IMF k t is the kth IMF, IMF k t −p is pth lag value of kth IMF ε k t is the residual of kth IMF, p and q are autoregressive and moving average lags. Moreover, in some cases time series data is not stationary. To make such stationary, differences at an appropriate degree are used (Box & Jenkins, 1970). If such a situation occurs, then the model is known as ARIMA (p,d,q) where d is the differenced value used to make the non-stationary data stationary.

Selection of study area
In this research, the largest water system of Pakistan is considered for the application of the proposed strategy. The Indus Basin System (IBS) is the largest river of Pakistan and it plays an important role specifically in power generation and irrigation department. The major tributaries of IBS i.e., River Jhelum, River Chenab, and River Kabul are selected for the present study. Pakistan is facing a large amount of frequent river flooding each year due to monsoon rain and melting snow or glaciers. As in Pakistan, glacier-covered 13,680 km 2 area which is estimated 13% of the mountainous areas of Upper Indus Basin (UIB). Melted water from these 13% areas adds the significant contribution of water in these rivers. which leads to complex characteristics in river flow data. Therefore, for sustainable economic development and efficient water resources planning, there is a need for analyzing such complex characteristics and predict the behaviors of river flow data at IBS and its tributaries.

Data
To investigate the improved framework, four daily river flow data comprised on (1st-January to 19th-June) for the period of 2015-2018 is used in this study. The main river flow of Indus at Tarbela is considered with its three principal tributaries: Jhelum at Mangla, Chenab at Marala and Kabul at Nowshera. Data is measured in 1,000 ft/s. The described river flow data was obtained from the website of Pakistan Water and Power Development Authority (WAPDA).

Evaluation criteria
Evaluation of noise reduction methods: The performances of denoised series needs comprehensive evaluation after using appropriate noise reduction threshold methods. In this research, to check the performances of CEEMDAN/EEMD-ST CEEMDAN/EEMD-HT, CEEMDAN/EEMD-ITF, and proposed framework i.e., CEEMDAN/EEMD-EBT, Signal-to-Noise Ratio (SNR), Mean Square Error (MSE) and Mean Absolute Error (MAE) (Nazir et al., 2019) are employed which are given as follows respectively; and where y ot is the t th observed value and y pt is the t th predicted value. Moreover, the performance of proposed model ( where, y ot and y pt is defined above, y mo is the mean value of observed values, y mp is the mean value of predicted values and y sp is the standard deviation of predicted values. Results of decomposition stage: First, to confirm the non-stationarity of river flow data, Augmented Dickey Fuller (ADF) unit root test (Khalili et al., 2013) is used for all selected case studies. The test is applied to data by taking the log in order to confirm non-stationarity. Results of ADF showed that all selected river flow data i.e., Indus river flow, Jhelum river flow, Chenab river flow, and Kabul river flow is significantly non-stationary with p-value as 0.3353, 0.4135, 0.333 and 0.414 respectively. Then, the non-linear and non-stationary data decomposed into different time scale oscillation called IMFs to reduce the non-stationarity by extracting the time-varying characteristics of daily river flow data from all selected four stations. The CEEMDAN decomposition technique is used to extract IMF's of river flow data. All selected river flow data is decomposed into thirteen IMFs and one residual. The starting IMF's represents the highest frequencies whereas the last half IMF's showed the low frequencies and the residual represent the overall trend. The decomposed results of Indus River system is depicted in Fig. 2. The amplitude of white noise is set 0.2 as in Di, Yang & Wang (2014) and numbers of ensemble members are settled as 1,000. By inspecting the IMFs, it is noticed that each IMF component represents the oscillation characteristic Further, the cross-correlation method is employed to find the noisiest IMFs from all thirteen IMFs. To do this, first the decomposed IMFs are further divided into two parts to find the noisiest IMFs by using the cross-correlation between IMFs and original river flow data. The low correlation implies the more uncertainty present in IMFs. The first ten IMFs showed the least correlation with original river flow data indicated that these IMFs are overwhelmed with noises. The graph of cross-correlation between Indus river flow data and first IMF and between Indus river flow and eleventh IMF is depicted in left and right corner of Fig. 3, which shows that the first IMF is filled with noises with low correlation at all lags and eleventh IMF is free from noises as it showed 0.75 correlation not only at lag zero but also at other lag values. For all four rivers, first tenth IMFs are characterized as noisiest IMFs and last three IMFs are labeled as noise free IMFs. Results of denoising stage: the next step is to denoise the noisiest and sparse IMFs. To eliminate noises from IMFs, EBT estimator is used as a filter which assumes a mixture of prior as defined in Eq. (10), for each IMF separately by considering the nature of IMFs. First, the scaled transformation is applied to get normal distribution so that each IMF follows N (θ i ,1). According to the nature of IMFs as depicted in Fig. 4, it is known that most of the coefficients in all IMFs are zero and some are non-zero in which fewer coefficients are either very low or high. By looking (Fig. 4) both zero and non-zero part of IMFs, a mixture of an atom of probability at zero and multiple distributions are considered for non-zero part (Johnstone & Silverman, 2005). Among all of them, Laplace distribution is configured out as prior distribution of θ i with maximum SNR. Finally, the important coefficients of IMFs are preserved with posterior median threshold estimator described in Eq. (11) by attaining highest SNR value with minimum MSE and MAE given in Table 1 for all selected case studies. For comparison purpose, the conventional denoising methods as ST, HT and ITF are implemented on all river flow data. We observed that SNR of CEEMDAN-EBT Results of prediction stage: The decomposed and denoised IMFs of all selected case studies are further predicted through data-driven and statistical model. The denoised IMFs are predicted through MLP-neural network model. Training is performed by using forward and back-propagation by setting the learning rate parameter between 0.1 to 1. The back-propagation method with the optimal learning rate is selected to test the model. The second part of the decomposed IMFs comprised of noise-free IMFs and the residual, which are predicted through simple traditional statistical model (i.e., ARIMA (p, d, q)) for all case studies. The river flow data of all four rivers are splitted, 70%, 80% and 90% for the training set and 30%, 20% and 10% for testing set. The results achieved by such splitting criteria are not much deviated from each other so the only values of 80% training

DISCUSSION
Decomposition and denoising results: In order to understand the applicability of our proposed CEEMDAN-EBT model, the other EEMD-based decomposition method and three different denoising methods (i.e., EEMD-HT, EEMD-ST and EEMD-ITF were employed. The evaluation of denoising of proposed CEEMDAN-EBT and benchmark models were carried out using SNR, MSE and MAE measures for all river flow data. From    proposed strategy (i.e., CEEMDAN-EBT) performs better than that of all other methods (i.e., CEEMDAN/EEMD-ST, CEEMDAN/EEMD-HT and CEEMDAN/EEMD-ITF and EEMD-EBT based decomposition and denoising methods). Final prediction model: To verify the superiority of proposed (i.e., CEEMDAN-EBT-MM) strategy to model the complex river flow data, we choose the CEEMDAN/EEMD-ST-MM, CEEMDAN/EEMD-HT, CEEMDAN/EEMD-ITF AND EEMD-EBT-MM models to analyze the prediction results of non-linear and noise-corrupted data. Our proposed prediction framework based on decomposition and novel strategy of denoising performs well as compared to all other decomposition and denoising methods. As shown in Table 2 and Fig. 7, the proposed model (i.e., CEEMDAN-EBT-MM) exhibits a good prediction performance for Indus river flow with least MAE, MSE and MAPE values. The MAE and MAPE values of Jhelum and Chenab river flow are also better than the other models. For Kabul river flow, our proposed model showed its efficiency with the least MAPE value. The other models are not consistent in terms of efficiency, their prediction behaviors for each river vary as CEEMDAN-HT-MM shows effective performance for Indus river flow but was poor in predicting Kabul river flow. However, overall the proposed CEEMDAN-EBT-MM model shows consistent and excellence prediction results which indicate that with appropriate decomposition and the novel improvement of denoising technique, one can overall improve the performance of existing data-driven models to handle the river flow data.
In this study, the river flow of the selected Indus River system exhibited seasonal persistence at several multi-scales which manifested through data decomposition method, and was used to enhance the prediction accuracy of river flow. Moreover, instead of overall season, Kalhoro et al. (2017) shown that there is a further sub-seasonal variation i.e., high flow, called Kharif season, and low flow, called Rabi season, present in Indus River system. The majority of Pakistan's population lives in rural areas and have the occupation of farming, which is directly influenced by the severity of both sub-seasonal variations. In the future, there is a need for establishing reliable sub-season identification and prediction methods that will be fruitful for efficient water allocation in both seasons, which can ultimately boost the economy of Pakistan.

CONCLUSION
In this paper, due to the non-linearity and noises complexity of the river flow data, we proposed a strategy to improve the prediction of data-driven models with appropriate decomposition and a novel denoising method. Our proposed denoising method improves the performance of the CEEMDAN-based decomposition method by improving the time-scale components which enhance the prediction accuracy of data-driven models. The proposed method comprised of four steps such as decomposition, denoising and prediction, and aggregation. The performance of the proposed CEEMDAN-EBT-MM model is evaluated using four daily river flow data of IBS. The CEEMDAN-EBT-MM having the smallest MAPE values for all four case studies compared to other benchmark models. The improved results suggest that the proposed hybrid model can be used as an efficient tool for providing a reliable prediction of river flow data to policymakers for planning power generation and water resource management.