Considering the three intrinsic components (of autoregressive, seasonality, and error) of streamflow time series, the overall performance of the streamflow modeling tool is associated with the correct estimation of these components. In this study, a new hybrid method based on the wavelet transform (WT) as a multiresolution forecasting tool and exponential smoothing (ES) method, with two presented scenarios (WES1 and WES2), was introduced. To this end, the performance of the proposed method was investigated versus four conventional methods of the autoregressive integrated moving average (ARIMA), ES ad-hoc, artificial neural network (ANN), and wavelet-ANN (WANN) for daily and monthly streamflow modeling of West Nishnabotna and Trinity River watersheds with different hydro-geomorphological conditions. In the presented WES technique, firstly, WT is employed for decomposing the observed signal to one approximation (deterministic trend) and more diverse components of subseries (each at a specific frequency). Then, for the first scenario (WES1), only two subseries are introduced to the model as input parameters; however, for the second scenario (WES2), decomposed subseries are separately used as the inputs of ES models. The obtained results indicated that combining WT with the ES method and ANN led to more accurate modeling. The proposed methodology (WES2) that used all decomposed subseries separately improved the efficiency of models up to 30% and 10% for the daily dataset and up to 88% and 57% for the monthly dataset, respectively, for the West Nishnabotna and Trinity Rivers.

1. Introduction

Forecasting streamflow has been investigated by several researchers [15] as it is a fundamental subject in hydrological modeling. As a result, many researchers are developing new models to improve streamflow modeling. For instance, recently, new models were developed for streamflow forecasting by improving artificial intelligence (AI) models [6, 7]. Accurate modeling of the streamflow is a significant step in any study for a river improvement and management of the related watershed and includes highly nonlinear and interacting components that cause complexity. Difficulties in the estimation of each component of the process causing the estimation of the streamflow become notably difficult. As a result, the black-box models can be more efficient than the physical base models due to the complexity of streamflow time series [812]. For instance, linear and stationary time series can be properly modeled by the classic black-box methods by using seasonal auto-regressive integrated moving average (SARIMA) [13, 14]. But such models (e.g., SARIMA) may not be a suitable model for hydrological time series that are highly nonlinear and complex. As an alternative, nonlinear AI-based methods, especially the artificial neural network (ANN), have achieved real success in modeling hydrological time series because of its significant advantages that are as follows [15]:(i)Being black-box tools, it can be easily used in the absence of prior knowledge about the physical concepts of the process(ii)Having inherent nonlinearity, provided through activation functions and neurons, it can handle the nonlinearity of the process(iii)Possibility to have different input variables, thus it can consider time-space variability

Among widely adopted nonlinear models, the used time series are considered to be merely nonlinear but they may not be always purely nonlinear [16]. For practical matters, it is a challenging step to realize whether a time series is generated from a composition of some linear or nonlinear processes. So the application of such techniques with the capability of considering linear and nonlinear properties sounds a reasonable manner.

Exponential smoothing (ES) contributed to the various successful modeling techniques (considering linear and nonlinear methods). It was originally classified by McCormick’s [17] taxonomy, then developed by Gardner [18], adjusted by Hyndman et al. [19], and again expanded by Taylor [20] and Gardner [21], which resulted in thirty different models. While the successful applications of the ES method have been reported widely in the field of financial studies and business areas, the available literature in the hydro-environmental applications is really limited. It is notable that only the linear models of the ES method can be taken as an example of the ARIMA approach. The linear ARIMA model may not be able to handle the nonlinear behavior of the process; in this case, the use of ES models that consider both linear and nonlinear behaviors of the process seems to be a reliable alternative. Forecasts produced using the ES method are weighted averages of past observations, with the weights decaying exponentially as the observations get older [22]. Ord et al. [23] and Hyndman et al. [19] have stressed that all ES models (including nonlinear models) could lead to optimal forecasts from the innovation of state-space models. More information can be found in Gardner [18, 21]. Demand for extensive data samples for developing nonlinear models is a well-known issue, and the pure linear ARIMA models may not also be a good choice for handling the behavior of those processes [24]. The application of ES models seems to be a more logical alternative since it is capable of considering the linear and nonlinear nature of the datasets [25].

On the contrary, due to the nonstationarity and climate variability on many time scales (diurnal, subseasonal, seasonal, low-frequency, etc.), the ability of mentioned models (i.e., ANN, SARIMA, and ES) to forecast highly nonstationary and seasonal hydrological time series (e.g., streamflow) may be restricted due to the multifrequency nature of the real hydrological process. To handle the mentioned nonstationary problem, the application of wavelet-based data preprocessing has been already proposed and used successfully in hydro-environmental modeling [2629]. For example, Jamei et al. [27] developed wavelet-multigene genetic programming for the simulation of surface water total dissolved solids. Ahmadianfar et al. [26] developed the locally weighted linear regression with the wavelet transform (WT) for the prediction of electrical conductivity in surface water. Ni et al. [29] used a hybrid wavelet and long short-term memory (LSTM) network for monthly streamflow and rainfall forecasting and indicated the superiority of the hybrid model against conventional LSTM.

In the presence of highly nonstationary and seasonal variations, data preprocessing methods may be used to get a better model performance. The WT is a suitable approach since it is capable of discerning underlying frequencies and extracting more information from the large datasets. Thus, the hybrid wavelet-ANN (WANN) model has shown good performance and is widely applied for hydrological modelings such as rainfall runoff, groundwater, and river flows [10, 3032].

Therefore, to the simultaneous use of WT and ES method capabilities for streamflow modeling as a complex process, the wavelet-ES (WES) method is to model the river streamflow. WES1 and WES2 were considered as two approaches of the WES method. The aggregated details’ subseries and approximation have been used for modeling via WES1, while the WES2 method employed all subseries as inputs of the ES method.

This study is the first investigation into the capability of the WES technique in streamflow modeling, and therefore, the purpose of this research is to investigate the capability of the wavelet-based models (i.e., WES1, WES2, and WANN), ANN, ES, and SARIMA (as benchmark models) in streamflow modeling. To this end, datasets from two different watersheds (in terms of geomorphological conditions) have been used at daily and monthly scales, and the obtained results were compared together.

2. Study Area and Dataset

In the current research, the data from 2 case studies, West Nishnabotna River, a subbasin in southwestern Iowa, and Trinity River, a subbasin in California, United States, were adopted for implementing the presented techniques (see Figure 1).

The used streamflow time series were obtained from the USGS (https://waterwatch.usgs.gov/). The streamflow time series data of 35 years from 1981 to 2015 were employed in the modeling process, with 75% divided as the training dataset and the rest 25% considered as the test dataset. Table 1 presents the statistical characteristics of streamflow data of watersheds (e.g., mean, maximum, coefficient of variation, etc.).

Trinity River is one of the important tributary of the Klamath River which flows through the Coast Ranges and Klamath Mountains (in northwestern California) and located at 41° 11′ 5″ North latitude and 123°42′ 31″ West longitude. The area of the watershed is about 7800 km2 (Figure 1(a)). The altitude of the river fluctuates between 58 m and 2709 m above sea level and contains oak, fir, and pine forests (about 92%).

The West Nishnabotna River is one of the two streams of the Nishnabotna River that rises in southwestern Carroll County, flows about 190 km long and located at 40° 34′ 00″ North latitude and 95°39′ 51″ West longitude (Figure 1(b)). The watershed area is approximately 7600 km2 with the field land cover. Broad and open valleys with rolling uplands form the topography of the Nishnabotna River watershed.

The studied watersheds are located at a distance of about 1460 km from each other and have different geomorphological conditions. As shown in Table 1, it is clear that the Trinity River has a greater mean monthly streamflow than that of the West Nishnabotna River. The Trinity River Watershed consists of mixed conifer and evergreen brush at the lower elevations with true fir and lodgepole pine at the higher elevation. It consists of 76 percent mixed conifer, 10 percent shrub, 6 percent mixed fir, 5 percent nonforested, and 3 percent hardwoods (Upper Trinity River Watershed Analysis, USFS, 2005). In contrast to the forest land cover of the Trinity River Watershed, the West Nishnabotna River Watershed is characterized by its field land cover within only 10 percent forest or grassland [33].

The stated autocorrelation function (ACF) shows the seasonal frequency and autoregressive properties of watersheds’ time series (see Figure 2). It indicated a more regular pattern for the Trinity River since the denser land cover contributes to the regular pattern of the flow regime.

Additionally, considering the statistical characteristics of the watersheds reported in Table 1, the coefficient of variation for the West Nishnabotna River Watershed is higher than that for the Trinity River Watershed which denotes the fact that the West Nishnabotna River Watershed can be considered as a wild watershed in regards with Trinity River Watershed. It can be inferred that the West Nishnabotna River Watershed experiences a more irregular and nonlinear pattern in its involved hydro-environmental processes rather than well-dominated seasonal weather. Therefore, it is expected that the modeling results for the Trinity River Watershed will be better than the Nishnabotna River. The reported results of this study confirm this claim (see Section 4).

3. Methodology

3.1. Wavelet Transform

The continuous form of WT (W) for x (t) is shown as follows [34]:where n stands for the temporal translation of the function that helps to study the signal around n, m indicates the dilation and corresponds to the complex conjugate, and M (t) denotes the mother wavelet (wavelet function). Localizing the time scale of the process is the core objective of this transformation. For pragmatic cases, discrete-time signal approaches are more preferred by hydrologists in comparison with continuous one. A discrete mother wavelet can be expressed as follows [34]:where a as an integer plays a regulatory role for the wavelet dilation and b as another integer denotes the wavelet translation. The location is specified by the parameter n0 that is a positive value, and the fined dilation step is represented by parameter m0 that is a positive score higher than one. and are widely chosen scores. This logarithmic scaling for dilation and translation is known as dyadic grid arrangement, and the dyadic wavelet function is specified as follows [34]:

The dyadic WT for xt as a discrete-time series becomes [34]where stands for the discrete wavelet coefficient based on the scale and localization parameter ( and is an integer power of ). The overall trend of the signal (indicated by the smoothed part of the signal) is indicated by T. The formula for renovating the time series xt by the inverse discrete transform is [34]where is termed as the approximation series at level . Also, indicates the detail series with decomposition level a () and time dimension t ().

3.2. Autoregressive Integrated Moving Average

As the benchmark model, ARIMA and seasonal ARIMA (SARIMA) models were examined for the streamflow modeling of watersheds, in order to compare the results of linear and nonlinear models.

SARIMA models have been employed to analyze and forecast univariate hydro-environmental time-series data [10]. The SARIMA model is a developed form of the ARIMA model, which reflects the seasonal variation of the time series, and the original time series employs a lag operator to process SARIMA . The SARIMA model can be formulated as follows [35]:

In equation (6), B is the lag operator (defined as ); and are polynomials of orders and , respectively, and are polynomials in of degrees and , respectively, and d denote the order of nonseasonal auto-regression and the number of regular differences, indicates the order of nonseasonal moving average, , and show the order of seasonal auto-regression, the number of seasonal differences, and the order of seasonal moving average, respectively, and finally, stands for the length of the season. Thus, for the ARIMA model, it means that the seasonality parameters (, and ) are considered to be equal to zero.

In this paper, the ARIMA and SARIMA models were applied for streamflow forecasting using the software program. ARIMA and SARIMA models were first trained via the training dataset and then verified using the validation set.

3.3. Artificial Neural Network

ANN models are frequently utilized in dynamic nonlinear datasets’ modeling and achieve accurate results, particularly for complex hydrological phenomena in which the physics of problem’s variables are not entirely comprehended. In other words, ANNs are suitable for modeling the data-driven nature time series. An ANN is a composition of a set of neurons (or nodes, i.e., interconnected simple processing elements) with an attractive set of information processing characteristics such as nonlinearity, parallelism, noise tolerance, and learning and generalization capabilities.

Between the ANNs, the Feedforward Backpropagation (FFBP) Neural Network is the most common method utilized in solving different engineering problems [30, 36]. The term ‘‘Feed-Forward’’ means that a neuron in the input layer is only connected to other neurons in the hidden layer and neurons in the hidden layer are connected to neurons in the output layer. So there are no interconnections within each layer’s neurons.

An output value’s explicit expression form of a three-layered FFBP neural network is specified by [37]where denotes the activation function for the output neuron, is a weight vector in the output layer connecting the th neuron in the output layer and the th neuron in the hidden layer, indicates the activation function of the hidden neuron, is the weights applied in the hidden layer connecting the th node in the hidden layer and the th node in the input layer, is the th input variable for the input layer, is the bias for the th hidden neuron, is the bias for the th output neuron, and is the computed output variable. Finally, and are the number of nodes in the input and hidden layers, respectively. It should be noted that the data are normalized to become between 0 and 1 before feeding into the network.

Also, it is notable that, in this study for modeling streamflow via ANN tool, the MATLAB® software and its nntool were used.

3.4. Exponential Smoothing

ES is a practical approach for forecasting whereby the prediction is based on an exponentially weighted average of past observation [38]. Hence, by adopting the ES state-space technique, better results could be expected from the models developed for the datasets which pose nonstationary or nonlinear behavior [39].

Table 2 denotes the ES different models. For example, consider the cell and cell as an example; they represent the “additive Holt-Winters” and the “multiplicative Holt-Winters” model, respectively. According to the literature, the Holt-Winters (HW) model which has easier formulation and simultaneously handles the autoregressive and seasonal components of time series, among different forms of ES approach, is a commonly applied model [18, 40]. There are two possible state-space models for each of the fifteen models in Table 2, one for a model with additive errors and the other corresponding to the model with multiplicative errors. In our work, the additive and multiplicative errors have been applied for streamflow simulation, and thirty class of ES models (15) × (2) have been used, and the most reliable results have been presented in the Results and Discussion section.

The () states the components of error, trend, and seasonality [39]. Therefore, e.g., the ETS () represents a model with multiplicative errors and a multiplicative trend with no seasonality. It should be mentioned that additive or multiplicative determines linear or nonlinear attitude.

The equations of the HW model can be expressed as follows (equations (8)–(17)), and readers can find more details about other classes of ES on Hyndman and Khandakar [39].

The HW approach includes trend and seasonal () factors with , , and as smoothing parameters where the trend factor formed from growth and level terms. The additive HW model is as follows [22]:and the multiplicative HW model can be defined as follows [41]:

As already stated, and denote the smoothing parameters and, respectively, are affected by level, slope, and seasonality (). The parameters and are calculated by applicable error functions such as root mean square error (RMSE). For obtaining the seed values, the following formula can be used [41]:since

where states the mean value of , N states the cycles number, and C states the length of the seasonal cycle.

Also, in this study for modeling streamflow via the ES method, the R software was used. To this end, the “ETS” function from the “forecast” package was used. The “ETS” function tests all possible ES models (see Table 2) and selects the best model by optimizing parameters based on the mean square error (MSE) criterion [39].

3.5. Proposed Methodology

To use the advantages of ES and WT methods at the same time, the hybrid WES model is adopted in the presented study. As a temporal preprocessor, the WT is employed for detecting and extracting the existed seasonality and temporal features in the streamflow time series. The WT decomposes the observed streamflow time series into an approximation subseries and detail subseries, , where states the decomposition level, and each subseries implies an individual level of seasonality. Using the WT tool, the resulted subseries are fed into the univariate ES method. According to the univariate nature of the ES method, two scenarios were designed to use the WT resulted subseries. Since the WT decomposes the time series into two components, an approximation subseries and fluctuation subseries (which include several subseries representing the seasonal periods), the designed scenarios must consider both components separately. Hence, two scenarios were designed to model the approximation and fluctuation components separately, with the difference that the WES1 models the approximate subseries and aggregated fluctuation subseries (two ES models) but WES2 models the approximate subseries and all fluctuation subseries individually (i+1 ES models).

The developed models of WES1 and WES2 are stated below.

The outline of the WES1 is described in three steps as follows (Figure 3(a)):(i)The fluctuation series () is computed through (ii)Employing two distinct ES models, both of approximation and fluctuation series are modeled individually; consequently, both of and are obtained as computed series(iii)In the end, the overall computed streamflow () is achieved by

The outline of the WES2 is described in two steps as follows (Figure 3(b)):(i)All of the resulted subseries are modeled individually by employing , ES models (). Consequently, the computed approximation series and detail series are obtained.(ii)In the end, the computed streamflow time series is obtained using .

As it can be seen in Figure 3, the major contrast between the two WES1 and WES2 models is the number of utilized ES models. The WES1 uses two ES models for predicting the streamflow but the WES2 uses models. In other words, WES1 uses the and as the inputs, while the inputs of WES2 are and subseries.

3.6. Efficiency Criteria

In this study, RMSE and Nash–Sutcliffe (NSE) were used as efficiency criteria. The RMSE (equation (18)) and Nash–Sutcliffe (equation (19)) are broadly employed as evaluating criteria of the efficiency performance of hydrological models [42]:

In equations (18) and (19), N, Oi, Ci, and , respectively, indicate the number of samples, observed values, calculated values, and average of the observed values. Legates and McCabe [43] revealed that RMSE and NSE can evaluate the performance of the hydro-environmental models appropriately. Furthermore, considering the importance of peak values in modeling as well, equation (20) may be employed to assess the performance of the model in terms of capturing the extreme values [42]:where indicates the Nash–Sutcliffe criterion of peak values, stands for the number of peak values, indicates the observations, and and are computed and average of peak values, respectively.

4. Results and Discussion

In this research, the ES, ARIMA, ANN, and WANN as benchmark models were developed beside the WES method (with two considered scenarios) to one-time step ahead modeling of daily and monthly streamflow time series (see Figure 4). In the following, the results of the WES method are presented and discussed; then, a comparison with benchmark models is performed.

For each daily data and monthly data, the resolution levels of 2 to 8 and 2 to 5 were adopted, respectively. In early studies, the optimum decomposition level was usually determined through a trial-and-error process, but afterward, a formula which relates the minimum level of decomposition, L, to the number of data points within the time series Ns was introduced in the literature [30]:

So, in this study, the initial point of view to select of the decomposition level was taken from L but since many seasonal characteristics may be embedded in hydrological signals, 2–8 resolution levels for the daily and 2–5 resolution levels for the monthly modeling were examined via the proposed WANN and WES models which, respectively, denote to the 22-day mode and 23-day mode (which is nearly weekly mode), 24-day mode (which is nearly semimonthly mode), 25-day mode (which is nearly monthly mode), 26-day mode, 27-day mode (which is nearly semiyearly mode), and 28-day mode (which is nearly yearly mode) in the daily scale and 22-month mode, 23-month, 24-month, and 25-month mode in the monthly scale. Besides, the Daubechies 4 wavelet (db4) that has been frequently assessed in hydrological modeling was considered as the mother wavelet in this study. Based on the literature, in contrast to the importance of the decomposition level, the type of the mother wavelet poses a minor effect since the model can adjust the output [30]. Finally, by applying the models (ES and ANN) to the inputs obtained by wavelet decomposition, the output of hybrid models (WES and WANN) as can be computed. It should be mentioned that, for the daily and monthly scale, the 4th decomposition level offered better results. Thus, the outcomes of the 4th level are presented in the tables.

To increase the model efficiency regarding the multifrequency feature of the WT, the preprocessed data by WT were introduced to the ES method (with thirty types of models). The results of WES1 and WES2 modeling scenarios were, respectively, presented in Tables 3 and 4 . It is obvious that, for both watersheds, the second scenario (WES2) offers more precise results in comparison with the first scenario (WES1) because the WES2 uses all resulted subseries by WT as inputs but in modeling by WES1, by aggregating the details’ subseries and using only two inputs (approximation and fluctuation subseries), the effect of WT becomes weaker.

A noticeable point about the ES method parameters in Tables 3 and 4 (also in Table 5) is the values of the as level (trend), as the slope (trend growth), as seasonal, and as damping factors (). The level index () gives an estimate of the local mean, the slope index () denotes the change between successive time points, the seasonality index () estimates the deviation from the local mean due to seasonality, and finally, the damping index () is related to the stationary state of the time series, meaning that if the time series do not converge to a stable value as t increases limitlessly, the time series is nonstationary and should get a high value (near 1). In other words, a small value of means that older observations in the time series are weighted markedly but values near 1 mean that the latest observations will get higher weights [25].

As aforementioned, parameter indicates the seasonal component of the time series. However, the WT already extracted the seasonal component (as detail subseries), and the resulted subseries by WT do not have extractable seasonal components. Therefore, as can be seen in Tables 3 and 4, the value of the parameter is zero or very low (especially in daily models and the approximation component in monthly models). But the approximation subseries indicates the major trend of the original time series, so the and parameters of its ES model take their highest value (very close to 1).

The most efficient ES models of the two first levels of the resulted subseries (by WT) express the high-frequency seasonal characteristics with an insignificant trend element, so contains no and a lower score (see Table 4) so that modeling the first and second detail subseries within notable variations was considered a challenging step.

As the benchmark models, ANN, WANN, SARIMA, and ES were also examined for the streamflow modeling of the watersheds.

To develop the model with the ANN, up to n lag times (Qt, Qt 1, ..., Qt  n) of past observations were considered as possible inputs of ANN to forecast the Qt +1. Due to the selection of best inputs to modeling via AI-based models (multivariate models), the Pearson correlation coefficient (CC) and mutual information (MI) were used [44], and according to the efficiency criteria (RMSE and NSE), only the best results were reported in Table 6. Since the proper selection of input parameters and the network structure (i.e., number of layers, activation function, epoch number, and hidden neurons) deeply influence the AI-based model, problems such as underfitting and overfitting can be avoided by suitable selection of the network structure. The Levenberg–Marquardt backpropagation algorithm was utilized to train the network. To this end, tangent sigmoid has been used as the activation functions, and up to quadruple of the inputs were considered as the hidden neurons, and up to 500 were examined as iteration epochs’ number [37]. The most efficient ANNs’ outcomes are stated in Table 6.

The ability of WT in the multiresolution analysis could be used in dealing with seasonal patterns. Since the combination of WT and ANN methods reinforce the capability of the model in dealing with seasonality, in this study, for daily and monthly modeling, the WANN method was utilized. So, by using the feature selection methods (CC and MI) to select inputs and constructing the networks, Q (t + 1) can be predicted. The tabulated results in Table 6 show the efficient performance of WANN.

For both watersheds, ARIMA and SARIMA performed better, respectively, for daily and monthly scales (see Table 7). However, at monthly modeling for the low NSE value (<0.5), the considered models could not make a good fit on the streamflow and thus are not reliable models. The streamflow time series intrinsic nonlinear nature and its multifrequency seasonality patterns may be the reason for the lack of accuracy in modeling by linear and unifrequency analysis (SARIMA).

The performance results of the ES method are stated in Table 5. As can be seen, for all of the case studies, the ES model offered better outcomes in comparison to the ARIMA and SARIMA. But regarding the low value of NSE (≈0.5) in monthly modeling, this model did not present a good capability to model the streamflow of these rivers. The multifrequency nature of the signal may be the reason for the resulted deficiency.

The hybrid WES2 model performed better than other single models up to 30% and 10% (daily scale) and up to 88% and 57% (monthly scale) for the West Nishnabotna and Trinity, respectively. Especially, WES2 was more efficient in the face of the peak values for both watersheds (Figures 58 ).

The WANN and the WES2 are the potential methods to deal with autoregressive features and seasonality fluctuations, which are widely seen components at both daily and monthly streamflow time series (see Figure 2). Hence, the WANN and WES2 offered more reliable results in comparison with other methods such as the ARIMA (an approach that is capable of handling with the autoregressive only), the ES and SARIMA (uni-frequency analysis tools), and the WES1 model (the detail subseries are aggregated, and two inputs are introduced to the model with a lower impact of WT), that is, because the WT method is capable of handling complicated streamflow processes due to its potential in the multiresolution analysis. Figures 58 depict this superiority for both watersheds at both scales. In other words, it means introducing of all the resulted subseries from decomposing the main signal using WT, as the inputs would increase the efficiency of the model since each subseries would represent a specific seasonal scale and the various seasonality characteristics would be covered in the model. As the seasonality features become more highlighted and common at the monthly scale, this process becomes more sensitive [25].

The computational cost for WANN was more than that of the WES2 model because ANN has more parameters than the ES model (e.g., input numbers, weights, biases, hidden neurons number, and transfer functions’ type), while the performance of the WANN and WES2 model was more or less in the same range; ES models offered more reliable results compared to ARIMA and SARIMA due to the comprising linear and nonlinear kernels and their capability in analyzing time series regarding the primary factors of time series (trend, seasonality, and error).

The daily ARIMA model led to reliable results (because of the autoregressive characteristic of streamflow daily time series, see Table 7); however, for the monthly time series, the SARIMA model could not gain acceptable results (see Table 7), that is, maybe because of the existence of the multifrequency pattern of the streamflow monthly time series.

For a more comprehensive presentation of the performance of the WES2 model in streamflow modeling, a Taylor diagram (Figure 9) was developed by considering the CC, standard deviation, and RMSE between the observed and computed streamflow time series. In the Taylor diagram, the standard deviation of the pattern is proportional to the radial measured from the origin, the correlation between the two fields is given by the azimuthal position of the test field, and the centered RMSE value is proportional to the distance between the actual and the estimated fields with the same units as the standard deviation [45]. When the predicted value is closer to the observed value in the matter of CC and standard deviation, one can expect more accurate prediction results. The value of the RMSE decreases with an increase in correlation. From the result shown in Figure 9, it can be clearly seen that the WES2 and WANN models demonstrated higher capability than the other used models for the streamflow modeling in terms of RMSE and pattern correlation (also standard deviation).

It can be concluded that, while for the time series in the monthly scale of the watershed, the WANN, WES2, ANN, ES, WES1, and SARIMA presented a more reliable performance; WES2, WANN, ANN, ES, WES1, and ARIMA models led to more acceptable outcomes in the daily scale, respectively.

From the aspect of spatial assessment, the studied watersheds’ land covers have different conditions. As stated in Tables 37, it is obvious that, by employing various time scales, West Nishnabotna River efficiency’s increase was more than Trinity River. Also, for all modeling procedures, the Trinity River got more reliable outcomes because the Trinity River Watershed has a more dense land cover (covered by forests and is not much affected by anthropogenic influences), and its reaction to rainstorm (or precipitation) is not so complex (nonlinear) as much as West Nishnabotna. Add to that, as it is clear in Figure 2, the Trinity River posed a more powerful “autoregressive characteristic” than the West Nishnabotna River. In fact, time-series within moderate “autoregressive characteristic” show relatively weaker performance in the modeling process. However, it should be mentioned that WT contributes to the more enhancement of performance for the cases that have a moderate autoregressive. So, in this study, WT led to the enhancement of efficiency for West Nishnabotna River streamflow modeling, which was more than that for the Trinity River.

The hydrological processes such as streamflow include several lumped and physical-based parameters, and it will be quite impossible to perform exact estimation of each parameter of the process. As a result, the black-box modeling can be more effective than the fully distributed models because of complex unknown and abstruse factors involved in streamflow nature. The black-box modeling (such as the present modeling) is case sensitive, and although the same methodology can be applied to other cases, the results and calibrated parameters cannot be similarly used for other cases. However, as presented in this study by applying to two different cases, it will be possible to apply the presented methodology to other watersheds by calibrating the parameters of the models. Although it seems the proposed methodology may include a high degree of freedom (the ES parameters selection, mother wavelet selection, level of decomposition, etc.), no need to calibrate all because some of them (e.g., wavelet type, mother wavelet type, and decomposition level) may be determined by prior knowledge (according to the physics of the process) or by referring to the literature [30].

Regarding the amount of data for streamflow modeling, this can be seen from two different aspects: diversity of datasets (having data from several cases and watersheds) and the length of the data itself (for each case). Referring to the aim of this study (which is to introduce a newly proposed methodology in streamflow modeling), it is revealed that at least one case study is needed to verify the applicability of the proposed model, and then, others may apply to other cases. In this way, the methodology was applied on two distinct case studies to see the performance and validity of the proposed model at different conditions. In spite of the different hydrological natures of these two case studies, results showed the priority of the model (WES2) and its capability in streamflow modeling. However, regarding the length of the used dataset, the dataset includes about 35 years of daily record; referring to literature in streamflow modeling, it is revealed that the length of data is suitable for the hydrological modeling [28, 46, 47]. It is noticeable that the longer period of historical records would be usually better for modeling hydrological phenomena, but the used historical records are taken due to their quality and the minimum gap existence in the records (particularly as the longest available period).

5. Conclusions

Due to the importance of hydrological time series and their complex intrinsic components (autoregressive, seasonality, and error), the performance of the hydrological models depends on their ability in dealing with these components. In this research, the capability of ARIMA (a classic autoregressive model), SARIMA (a classic autoregressive model which also considers uni-periodicity of the process), ES (a nonstationary but uni-frequency modeling approach), ANN (a nonlinear model), WANN (a multifrequency preprocessed nonlinear model), and proposed WES with two different scenarios were examined for daily and monthly streamflow modeling of two different watersheds.

In the presented WES technique, after using WT for decomposing the observed time series, for the WES1, only two input subseries were introduced to the model as input parameters, while the WES2 used all decomposed subseries separately as the inputs of developed ES models. Finally, by aggregating the outputs of all ES models, the outputs of the streamflow modeling were obtained.

The WES2 hybrid model performed better than other used models (single models) up to 30% and 10% (daily time scale) and up to 88% and 57% (monthly time scale), respectively, for West Nishnabotna and Trinity rivers. Also, the WES2 in comparison with other methods presented better results in estimating the extreme points.

Comparing the results, it can be concluded that merging the capability of the WT method with ES (WES2) and ANN (WANN) contributed to more accurate models to predict the streamflow. The performances of the WANN and WES2 model were more or less in the same range but the computational cost for the WANN model was more that of WES2 because ANN has more parameters than the ES model (e.g., input numbers, weights, biases, hidden neurons number, transfer functions’ type, etc.).

To spatial assessment, modeling enhancement for the West Nishnabotna River was more than that of Trinity River due to its different geomorphological conditions, which was deliberated in Section 2.

The depicted reliable outcome of WES2 application for streamflow forecasting recommends its employment for other processes such as wastewater treatment, air pollution, and groundwater. Moreover, applying the stated methodology of streamflow modeling on the other watersheds would be helpful for evaluating the impact of geomorphological conditions on the efficiency of the modeling.

Data Availability

The used streamflow time series were obtained from the USGS (https://waterwatch.usgs.gov/).

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.