*2.4. Wavelet Transforms*

The wavelet transform is a mathematical function that has an adjustable time-frequency window and can decompose time series into multiple resolution levels by controlling the scaling and shifting factors of a mother wavelet [46]. A mother wavelet needs to be determined before applying a wavelet analysis. The wavelet transform of time series data generates sets of wavelet coefficients for different scales and provides a time-scale localization of processes [7]. The wavelet transform has two forms: the continuous wavelet transform (CWT) and the discrete wavelet transform (DWT). In this paper, we applied DWT method to build the hybrid SWAT-WSVR model. The DWT discretizes the parameters of scales and positions before implementing the wavelet transform to decrease the redundancy. The DWT of signal *f*(*t*) is defined as [47]:

$$\mathcal{W}\_f(j,k) = a\_0^{-j/2} \int\_{-\infty}^{\infty} \psi^\* \left(a\_0^{-j}t - kb\_0\right) f(t)dt\tag{4}$$

where the dilation parameter *a* and temporal translation parameter *b* of the CWT are discretized as *a* = *a j* 0 , *b* = *kb*0*a j* 0 , *a*<sup>0</sup> > 0 *and a*<sup>0</sup> 6= 1, *b*<sup>0</sup> ∈ *R*. In most cases, parameter *a*<sup>0</sup> = 2 and *b*<sup>0</sup> = 1. Then, a discrete wavelet can be expressed as [47]:

$$
\psi\_{j,k}(t) = a\_0^{-j/2} \psi \left( a\_0^{-j} t - kb\_0 \right) \quad j, k \in Z \tag{5}
$$

The DWT obtains wavelet details (D) and approximations (A) of the original hydrological time series through high-pass and low-pass filters, respectively. Approximations

at various resolution levels can be further decomposed by high-pass and low-pass filters (Figure 3). Commonly used DWT wavelets have 'Daubechies', 'Symlets', 'Coiflets', and 'Biorthogonal'. More details about wavelet transform can be found in Labat (2005) [48] and Mallat (2009) [49]. at various resolution levels can be further decomposed by high-pass and low-pass filters (Figure 3). Commonly used DWT wavelets have 'Daubechies', 'Symlets', 'Coiflets', and 'Biorthogonal'. More details about wavelet transform can be found in Labat (2005) [48] and Mallat (2009) [49].

ି/ଶ൫

The DWT obtains wavelet details (D) and approximations (A) of the original hydrological time series through high-pass and low-pass filters, respectively. Approximations

ି − ൯ , ∈ (6)

*Water* **2022**, *14*, x FOR PEER REVIEW 8 of 19

,() =

**Figure 3.** The discrete wavelet transforms (DWT) decomposition procedure of flow time series. **Figure 3.** The discrete wavelet transforms (DWT) decomposition procedure of flow time series.

We examined Daubechies wavelet family with different wavelengths, such as extremal phase filter with length 1 ('haar' or 'd1′), filter with length 2 ('d2′), filter with length 4 ('d4′), filter with length 6 ('d6′), and Least Asymmetric filter with length 8 ('la8′), and found that 'haar' wavelet is suitable for this study. To date, there is not a standard method to determine the optimum DWT decomposition levels (*D*) for a specific time series. Some applied the equation as = (()) [50]; others used the formula as = (/(2 − 1))/(2) [27], where *N* is the length of monthly time series, and *m* is We examined Daubechies wavelet family with different wavelengths, such as extremal phase filter with length 1 ('haar' or 'd1'), filter with length 2 ('d2'), filter with length 4 ('d4'), filter with length 6 ('d6'), and Least Asymmetric filter with length 8 ('la8'), and found that 'haar' wavelet is suitable for this study. To date, there is not a standard method to determine the optimum DWT decomposition levels (*D*) for a specific time series. Some applied the equation as *D* = *int*(*Log*(*N*)) [50]; others used the formula as *D* = *Log*(*N*/(2*m* − 1))/*Log*(2) [27], where *N* is the length of monthly time series, and *m* is the number of vanishing moments of a Daubechies wavelet.

the number of vanishing moments of a Daubechies wavelet. In this study, the maximum and minimum length of the monthly flow time series of 12 sites is 228 and 42 (Table 1), respectively. Regardless of applying either the formula mentioned above to calculate *D*, the maximum decomposition levels of monthly flow for all sites were between 1.62 and 7.83. Therefore, wavelet decomposition levels of 1 to 7 were tested to obtain the optimal resolution levels. The results indicated that the decom-In this study, the maximum and minimum length of the monthly flow time series of 12 sites is 228 and 42 (Table 1), respectively. Regardless of applying either the formula mentioned above to calculate *D*, the maximum decomposition levels of monthly flow for all sites were between 1.62 and 7.83. Therefore, wavelet decomposition levels of 1 to 7 were tested to obtain the optimal resolution levels. The results indicated that the decomposition level of 3 and 2 attained the best model performance.

#### position level of 3 and 2 attained the best model performance. *2.5. Model Performance Evaluation*

*2.5. Model Performance Evaluation*  We applied four statistics such as *NSE* (Nash–Sutcliffe efficiency), *PBIAS* (percent bias), *RMSE* (root mean square error), and *RSR* (*RMSE*-observation's standard deviation ratio) to evaluate the model performance in calibration and validation. Table 3 listed these statistical indicators, their mathematic expressions, and their value range. Preferred sta-We applied four statistics such as *NSE* (Nash–Sutcliffe efficiency), *PBIAS* (percent bias), *RMSE* (root mean square error), and *RSR* (*RMSE*-observation's standard deviation ratio) to evaluate the model performance in calibration and validation. Table 3 listed these statistical indicators, their mathematic expressions, and their value range. Preferred statistics combination is the lower *RSR*, *PBIAS*, and *RMSE* but the higher *NSE*, which present the better the model prediction performance. We used the 'hydroGOF' package [51] in R to calculate the mentioned statistical indicators.

tistics combination is the lower *RSR*, *PBIAS*, and *RMSE* but the higher *NSE,* which present the better the model prediction performance. We used the 'hydroGOF' package [51] in R to calculate the mentioned statistical indicators. To further compare the performance of SWAT-CUP, SWAT-SVR, and SWAT-WSVR on the entire time series (i.e., combined calibration and validation together), we plotted hydrography for each site and applied the Taylor diagram [52] to examine the relative importance of different statistics such as *r*, *RMSE*, and *NSD* between the observed and simulated flow for three models and twelve sites. The advantage of the Taylor diagram is that it can highlight the goodness-of-fit of multiple models and compare their difference from observed data at the same graph.


**Table 3.** Evaluation indicators of the model performance and their mathematic expressions.

Note: *y<sup>i</sup>* is the observed data series, *y* 0 *i* is the simulated results series, the overbar represents the mean value of data series, and *n* is the sample number.

#### **3. Results and Discussion**

We developed a total of 72 models for the flow prediction at 12 sites by three methods: SWAT-CUP, SWAT-SVR, and SWAT-WSVR. Each method included 12 calibrated models and 12 validated models. The difference between SWAT-SVR and SWAT-WSVR is that model inputs of SWAT-SVR had only the flow outputted from SWAT-CUP (a calibrated SWAT model) and precipitation data. Instead, we replaced the simulated flow with its wavelet components at different resolution levels in SWAT-WSVR. Table 4 lists the statistical performance of the three above-mentioned models in calibration, validation, and the whole time series data combining calibration and validation time series data.

#### *3.1. Flow Prediction by SWAT-CUP*

Table 4 summarizes the average *RSR*, *NSE*, *PBIAS*, and *RMSE* for twelve sites from SWAT-CUP, and the corresponding values are 1.67, 0.22, 57.57, and 11.50 m<sup>3</sup> s −1 in calibration and 0.84, 0.26, 35.98, and 11.04 m<sup>3</sup> s −1 in validation, respectively. The *RSR*, *NSE*, and *RMSE* had approximately similar performances between calibration and validation, but the value of *PBIAS* in validation was lower than one in calibration, which indicated that the predicted discrepancy from validation was less than one from calibration. SWAT-CUP overestimated monthly flow in both calibration and validation. The low average *NSE* value (≤0.26) indicated that SWAT-CUP has a poor goodness-of-fit between the observed and simulated flow for both calibration and validation. Additionally, the 07195430 site had the best performance among all sites in validation with lowest *PBIAS* (−6.7) and *RSR* (0.58) and the highest *NSE* value (0.66). The model performance of 07195500 and 07196500 were also acceptable in validation. After combined calibration and validation data together, the averages of *RSR*, *NSE*, *PBIAS*, and *RMSE* for 12 sites are 0.85, 0.13, 50.67, and 11.38 m<sup>3</sup> s −1 , respectively. Simulations from SWAT-CUP greatly overestimated the observed flow according to *PBIAS* (i.e., a positive mean of 50.67 for 12 sites) and presented a low fitting degree due to a low average value (0.13) of *NSE*. Overall, SWAT-CUP had a poor simulation performance for most of sites during both calibration and validation periods, with only few exceptions.

#### *3.2. Flow Prediction by SWAT-SVR*

Due to the unsatisfactory overall performance of SWAT-CUP, we developed the SWAT-SVR and SWAT-WSVR model integrating the simulated flow (or its wavelet components) and precipitation to improve the prediction accuracy. Table 5 listed the model structure and optimal parameter sets for SWAT-SVR and SWAT-WSVR. In 24 SVR models, the value of *C* is inconstant from 2.015625 to 255.015625, the value range of *γ* is from 1 to 14, and *ε* keeps a constant value of 0.00390625.


**Table 4.** Performance of flow calibration, validation, and combined data series on each site by SWAT-CUP, SWAT-SVR, and SWAT-WSVR.

† Note: The data series combined calibration and validation time series.


† Note: Flow comes from SWAT-CUP simulated discharge output. Ds and As are wavelet components from the simulated flow of SWAT-CUP.

The average *RSR*, *NSE*, *PBIAS*, and *RMSE* of twelve sites from SWAT-SVR are 0.58, 0.64, <sup>−</sup>9.48, and 7.34 m<sup>3</sup> s −1 in calibration and 0.80, 0.32, <sup>−</sup>20.68, and 11.30 m<sup>3</sup> s −1 in validation (Table 4), respectively. Compared with SWAT-CUP, the performance of the SWAT-SVR model had the lower *RSR* and the absolute value of *PBIAS* and the higher *NSE* in calibration and validation, particularly for the calibrated simulations, as SVR has a strong learning ability for training data. The average *RMSE* in SWAT-SVR calibration is lower: only 7.34 m<sup>3</sup> s −1 compared with the average one of 11.50 m<sup>3</sup> s −1 in the SWAT-CUP calibration. The results showed that the SWAT-SVR calibration on all sites had lower deviation and higher *NSE* value in comparison with SWAT-CUP, but still underestimated monthly flows on most sites. SWAT-SVR was generally superior to SWAT-CUP on all sites. However, only a few sites (e.g., 07195500, 07195430) had lower *RSR*, *PBIAS*, and higher *NSE* values, which indicated that SWAT-SVR could greatly improve the model performance in calibration but did not possess good generalization capability, which means it failed to keep this prediction ability with high accuracy while it was applied in validation. From the perspective of the whole data series, the average *RSR*, *NSE*, *PBIAS*, and *RMSE* are 0.65, 0.57, <sup>−</sup>13.03, and 8.84 m<sup>3</sup> s −1 , respectively. Clearly, compared with SWAT-CUP, the performance of the SWAT-SVR model were improved but limited, although SWAM-SVR generally had a low *RSR*, absolute value of *PBIAS*, and higher *NSE* value for calibration, validation, and both periods.

## *3.3. SWAT-WSVR Development and Evaluation*

#### 3.3.1. Development of SWAT-WSVR

We conducted DWT of the simulated flow from SWAT-CUP at twelve hydrological sites using a 'haar' wavelet filter to obtain the flow series structure, trend, and temporal characteristics. The wavelet decomposition was implemented at three resolution levels: 2 months, 4 months, and 8 months. Figure 4 showed an example of the flow DWT at the 07195430 site, including temporal features of the flow 2-month mode (D1), 4-month mode (D2), 8-month mode (D3), and approximate mode (A3). From Figure 4a, we can observe a large deviation between the original observed (blue line) and SWAT-CUP simulated (orange dash-line) flow signals. In Figure 4b, D1, D2, and D3 modes indicated high-frequency details, and A3 mode revealed a low-frequency trend and the slowest flow changing of the simulated flow series at the 07195430 site. These wavelet components would be used to build the SWAT-WSVR model afterward.

To reduce the computational load and find the most related wavelet coefficients to participate in the construction of SWAT-WSVR, we analyzed Pearson's correlation coefficient matrices between the observed monthly flow, monthly precipitation, and its sub-time series D1, D2, D3, and A3 from wavelet decompositions (Figure 5). In Figure 5, the upper 'Pie' graphs are a corresponding display of the lower 'numeric' *r* in the diagonal direction where the flow, precipitation, and wavelet components are shown. Flow time series of 07196090 and 07196973 were decomposed into two resolution levels since their data lengths are shorter: 42 and 96, respectively.

The correlation analysis indicated that the flow and precipitation have a strong positive correlation, and the average value of *r* is 0.53 on twelve sites. Most wavelet components D1, D2, D3, and A3 on the flow also showed positive correlations. The average of *r* of D1, D2, D3, and A3 (or A2) on flow is 0.35, 0.47, 0.3, and 0.32. The wavelet coefficient of D2 had the strongest correlation with the flow. We chose wavelet components with *r* greater than 0.2 to participate in SVR prediction so that we can keep the intrinsic nonlinear features in wavelet components as much as possible while avoiding large computational burden. The specific model structure of SWAT-WSVR for each site is listed in Table 5.

the simulated flow series at the 07195430 site. These wavelet components would be used

**Figure 4.** An example of the flow DWT at 07195430 site: (**a**) Observation and SWAT−CUP simulation, (**b**) wavelet decomposition. **Figure 4.** An example of the flow DWT at 07195430 site: (**a**) Observation and SWAT-CUP simulation, (**b**) wavelet decomposition.

#### 3.3.2. Statistical Evaluation of SWAT-WSVR

to build the SWAT-WSVR model afterward.

To reduce the computational load and find the most related wavelet coefficients to participate in the construction of SWAT-WSVR, we analyzed Pearson's correlation coefficient matrices between the observed monthly flow, monthly precipitation, and its subtime series D1, D2, D3, and A3 from wavelet decompositions (Figure 5). In Figure 5, the upper 'Pie' graphs are a corresponding display of the lower 'numeric' *r* in the diagonal direction where the flow, precipitation, and wavelet components are shown. Flow time series of 07196090 and 07196973 were decomposed into two resolution levels since their data lengths are shorter: 42 and 96, respectively. 07195800 07195855 07196000 Table 4 showed the statistical performance of SWAT-WSVR with the average *RSR*, *NSE*, *PBIAS*, and *RMSE* of 0.02, 1.00, <sup>−</sup>0.15, and 0.27 m<sup>3</sup> s −1 in calibration; 0.14, 0.98, −1.88, and 2.91 m<sup>3</sup> s −1 in validation; and 0.08, 0.99, <sup>−</sup>0.74, and 1.61 m<sup>3</sup> s −1 in the whole data series. Compared with SWAT-SVR and SWAT-CUP, SWAT-WSVR had the lowest *RSR*, the absolute value of *PBIAS*, and *RMSE* but the highest *NSE* value in validation. This result clearly indicated that SWAT-WSVR could effectively decrease the discrepancy of the simulation and obtain the best prediction accuracy for validation in comparison with SWAT-SVR and SWAT-CUP. Based on the value of *PBIAS*, SWAT-WSVR slightly underestimated the monthly flow in calibration. SWAT-WSVR also presented the best performance on the whole data series, along with the lowest *RSR*, *PBIAS*, and *RMSE* but the highest *NSE* in comparison with SWAT-CUP and SWAT-SVR. By comparison, the SWAT-WSVR model outperformed the SWAT-CUP and SWAT-SVR model in calibration, validation, and both periods.

07195500 07195430 07196090


data lengths are shorter: 42 and 96, respectively.

the simulated flow series at the 07195430 site. These wavelet components would be used

**Figure 4.** An example of the flow DWT at 07195430 site: (**a**) Observation and SWAT−CUP simula-

To reduce the computational load and find the most related wavelet coefficients to participate in the construction of SWAT-WSVR, we analyzed Pearson's correlation coefficient matrices between the observed monthly flow, monthly precipitation, and its subtime series D1, D2, D3, and A3 from wavelet decompositions (Figure 5). In Figure 5, the upper 'Pie' graphs are a corresponding display of the lower 'numeric' *r* in the diagonal direction where the flow, precipitation, and wavelet components are shown. Flow time series of 07196090 and 07196973 were decomposed into two resolution levels since their

to build the SWAT-WSVR model afterward.

tion, (**b**) wavelet decomposition.




The correlation analysis indicated that the flow and precipitation have a strong positive correlation, and the average value of *r* is 0.53 on twelve sites. Most wavelet components D1, D2, D3, and A3 on the flow also showed positive correlations. The average of *r*

greater than 0.2 to participate in SVR prediction so that we can keep the intrinsic nonlinear features in wavelet components as much as possible while avoiding large computational burden. The specific model structure of SWAT-WSVR for each site is listed in Table 5.

Table 4 showed the statistical performance of SWAT-WSVR with the average *RSR*, *NSE*, *PBIAS*, and *RMSE* of 0.02, 1.00, −0.15, and 0.27 m3 s−1 in calibration; 0.14, 0.98, −1.88, and 2.91 m3 s−1 in validation; and 0.08, 0.99, −0.74, and 1.61 m3 s−1 in the whole data series. Compared with SWAT-SVR and SWAT-CUP, SWAT-WSVR had the lowest *RSR*, the absolute value of *PBIAS*, and *RMSE* but the highest *NSE* value in validation. This result clearly indicated that SWAT-WSVR could effectively decrease the discrepancy of the simulation and obtain the best prediction accuracy for validation in comparison with SWAT-SVR and SWAT-CUP. Based on the value of *PBIAS*, SWAT-WSVR slightly underestimated the monthly flow in calibration. SWAT-WSVR also presented the best performance on the

3.3.2. Statistical Evaluation of SWAT-WSVR

#### *3.4. Taylor Diagram and Hydrographic Comparison between Different Models*

To further compare the overall performance of three models, we combined the flow observations and the calibrated and validated simulations from SWAT-CUP, SWAT-SVR, and SWAT-WSVR; recalculated statistical indicators including *r*, *RMSE*, and *NSD*; and re-estimated the model performance on the whole time series (i.e., calibration and validation periods are considered together). The Taylor diagrams (Figure 6) depicted the overall performance of different models for each site by identifying the pattern correlations, variability, and *RMSE* between observations and simulations. In Figure 6, both the x-axis and y-axis denote *NSD*; black dashed lines represent the *r* between observations and simulations; the normalized *RMSE* of the simulation is proportional to the distance from the *x*-axis identified as "observation" (green contours); the *NSD* of the simulation is proportional to the radial distance from the origin point (black contours). The modeling results in Figure 6 demonstrated that the developed SWAT-WSVR had the best performance in comparison with the other two models since it had the lower *RMSE* and *NSD*, but the higher *r* with observed flows in most cases. For example, SWAT-WSVR (redpoint in Taylor diagram) at all sites is closer to the reference point (observation) located on the *x*-axis than the other two models, which illustrated the SWAT-WSVR fitted best with the observed flow for most sites. The rank of the overall performance of three models from high to low follows SWAT-WSVR > SWAT-SVR > SWAT-CUP. The proposed SWAT-WSVR model that uses wavelet components as inputs of SVR presented a more satisfactory flow prediction than the SWAT-SVR with a single pattern input (i.e., the simulated flow from SWAT-CUP). Possible reasons include a periodical feature of sub-series represented by wavelet components is more obvious than those directly obtained from SWAT-CUP [53], and SVR captured the intrinsic nonlinear features between SWAT-CUP simulations and observed flows and built a mapping relationship at a higher dimension space successfully. This kind of fitting relationship is typically determined by using the trial-and-error method and a large number of iterations of many parameters sets in physically based hydrological models.

To investigate the entire and continuous performance of monthly flow prediction in the IRW, we plotted flow hydrography on the whole data series for each site (Figure 7). Here, we only labeled statistics of SWAT-WSVR for clarity, and other details related to estimate indicators can be found in Table 4. This figure reflected where the developed SWAT-WSVR model performed better than SWAT-CUP and SWAT-SVR methods. An oval region at the 07196090 site showed clear evidence that SWAT-WSVR agreed well with the observed flow, but SWAT-SVR missed the peak flow during this time window. Although SWAT-WSVR had desirable modeling performance, it missed few peak flows. For example, the SWAT-WSVR simulated flow (105.24 and 107.68 m<sup>3</sup> s −1 ) was 27.2% and 27.9% lower than observed flow (144.61 and 149.42 m<sup>3</sup> s −1 ) in April 2011 at 07195430 and 07195500 site, respectively. Overall, the other two models more or less capture the rising and recession of the observed monthly flow over time at all sites, but the SWAT-WSVR is more efficient at fitting with the observation and corrected errors compared to SWAT-CUP. This result is in line with others' conclusions that the application of wavelet transform in data-driven models can improve the accuracy of flow prediction [53,54]. The developed SWAT-WSVR model fit the observations well at all sites of the IRW based on the statistical results, Taylor diagram, and hydrography analysis.

In this study, we applied wavelet transforms on the simulated flow time series from SWAT-CUP rather than on the observed flow data to show that the developed SWAT-WSVR model can be applied in practice or future scenario prediction where observed data are impossible to access or has limited availability. An SVR coupled to a wavelet transform model approach based solely on observed monthly flow data could produce greater accuracy. This is because wavelet decompositions with less noise can represent the periodical and trend characteristics of flow series structure better than the original data series [25,53], and SVR has a strong learning capability to capture the corresponding relationship between wavelet decompositions and its original data series. However, the limitations of a purely mathematical or machine learning approach prior to constraint by a

physically-based model could fail to capture the relationship between rainfall and runoff and exhibit less predictive or unrealistic behavior [28]. method and a large number of iterations of many parameters sets in physically based hydrological models.

*Water* **2022**, *14*, x FOR PEER REVIEW 14 of 19

*3.4. Taylor Diagram and Hydrographic Comparison between Different Models* 

periods.

comparison with SWAT-CUP and SWAT-SVR. By comparison, the SWAT-WSVR model outperformed the SWAT-CUP and SWAT-SVR model in calibration, validation, and both

To further compare the overall performance of three models, we combined the flow observations and the calibrated and validated simulations from SWAT-CUP, SWAT-SVR, and SWAT-WSVR; recalculated statistical indicators including *r*, *RMSE*, and *NSD*; and reestimated the model performance on the whole time series (i.e., calibration and validation periods are considered together). The Taylor diagrams (Figure 6) depicted the overall performance of different models for each site by identifying the pattern correlations, variability, and *RMSE* between observations and simulations. In Figure 6, both the x-axis and y-axis denote *NSD*; black dashed lines represent the *r* between observations and simulations; the normalized *RMSE* of the simulation is proportional to the distance from the *x*axis identified as "observation" (green contours); the *NSD* of the simulation is proportional to the radial distance from the origin point (black contours). The modeling results in Figure 6 demonstrated that the developed SWAT-WSVR had the best performance in comparison with the other two models since it had the lower *RMSE* and *NSD*, but the higher *r* with observed flows in most cases. For example, SWAT-WSVR (redpoint in Taylor diagram) at all sites is closer to the reference point (observation) located on the *x*-axis than the other two models, which illustrated the SWAT-WSVR fitted best with the observed flow for most sites. The rank of the overall performance of three models from high to low follows SWAT-WSVR > SWAT-SVR > SWAT-CUP. The proposed SWAT-WSVR model that uses wavelet components as inputs of SVR presented a more satisfactory flow prediction than the SWAT-SVR with a single pattern input (i.e., the simulated flow from SWAT-CUP). Possible reasons include a periodical feature of sub-series represented by wavelet components is more obvious than those directly obtained from SWAT-CUP [53], and SVR captured the intrinsic nonlinear features between SWAT-CUP simulations and observed flows and built a mapping relationship at a higher dimension space successfully. This kind of fitting relationship is typically determined by using the trial-and-error

**Figure 6.** Taylor diagram of three models in comparison with the observation for each site. **Figure 6.** Taylor diagram of three models in comparison with the observation for each site.

To investigate the entire and continuous performance of monthly flow prediction in the IRW, we plotted flow hydrography on the whole data series for each site (Figure 7). Here, we only labeled statistics of SWAT-WSVR for clarity, and other details related to estimate indicators can be found in Table 4. This figure reflected where the developed SWAT-WSVR model performed better than SWAT-CUP and SWAT-SVR methods. An

with the observed flow, but SWAT-SVR missed the peak flow during this time window. Although SWAT-WSVR had desirable modeling performance, it missed few peak flows. For example, the SWAT-WSVR simulated flow (105.24 and 107.68 m3 s−1) was 27.2% and 27.9% lower than observed flow (144.61 and 149.42 m3 s−1) in April 2011 at 07195430 and 07195500 site, respectively. Overall, the other two models more or less capture the rising and recession of the observed monthly flow over time at all sites, but the SWAT-WSVR is more efficient at fitting with the observation and corrected errors compared to SWAT-CUP. This result is in line with others' conclusions that the application of wavelet transform in data-driven models can improve the accuracy of flow prediction [53,54]. The developed SWAT-WSVR model fit the observations well at all sites of the IRW based on the

statistical results, Taylor diagram, and hydrography analysis.

**Figure 7.** Comparison of SWAT-CUP, SWAT-SVR, SWAT−WSVR, and observed flow time series on each site. (Note: a vertical blue line separated the 70% training and 30% testing data. Only statistics from SWAT-WSVR were labeled for clarity.). **Figure 7.** Comparison of SWAT-CUP, SWAT-SVR, SWAT-WSVR, and observed flow time series on each site. (Note: a vertical blue line separated the 70% training and 30% testing data. Only statistics from SWAT-WSVR were labeled for clarity).

In this study, we applied wavelet transforms on the simulated flow time series from SWAT-CUP rather than on the observed flow data to show that the developed SWAT-WSVR model can be applied in practice or future scenario prediction where observed data are impossible to access or has limited availability. An SVR coupled to a wavelet transform model approach based solely on observed monthly flow data could produce greater accuracy. This is because wavelet decompositions with less noise can represent the periodical and trend characteristics of flow series structure better than the original data series [25,53], and SVR has a strong learning capability to capture the corresponding relationship Moreover, the current proposed SWAT-WSVR model can be regarded as a compromise method when we cannot attain a desirable prediction accuracy even after conducting extensive SWAT-CUP iterative computation, although we still encourage researchers to directly obtain a satisfactory prediction from SWAT-CUP if possible. In this work, the construction of SWAT-WSVR heavily depends on the procedure of SWAT-CUP parameter calibration. Alternatively, we can also build SWAT-WSVR based on wavelet decomposition of the initial output from SWAT without the procedure of calibration based on the methodology proposed if SWAT-CUP is not accessible.

between wavelet decompositions and its original data series. However, the limitations of

#### a purely mathematical or machine learning approach prior to constraint by a physically-**4. Summary and Conclusions**

based model could fail to capture the relationship between rainfall and runoff and exhibit less predictive or unrealistic behavior [28]. Moreover, the current proposed SWAT-WSVR model can be regarded as a compromise method when we cannot attain a desirable prediction accuracy even after conducting extensive SWAT-CUP iterative computation, although we still encourage researchers to directly obtain a satisfactory prediction from SWAT-CUP if possible. In this work, the construction of SWAT-WSVR heavily depends on the procedure of SWAT-CUP parameter calibration. Alternatively, we can also build SWAT-WSVR based on wavelet decomposition of the initial output from SWAT without the procedure of calibration based on the methodology proposed if SWAT-CUP is not accessible. **4. Summary and Conclusions**  This study developed the SWAT-WSVR monthly flow prediction model which was built on the basis of the SWAT and SVR with discrete wavelet transforms and investigated the performance and effectiveness of this model. Precipitation and wavelet components of flow outputted from SWAT-CUP were served as input variables into SWAT-WSVR. The methodology loosely integrated the physically based model and the data-driven model. The proposed SWAT-WSVR model had the best statistical performance with the lower *RSR*, the absolute value of *PBIAS*, *RMSE*, and higher *NSE* in comparison with regular SWAT-SVR and SWAT-CUP, which indicated that SWAT-WSVR possessed the lower discrepancy and higher goodness-of-fit between the simulated and observed flow. The rank of the overallperformance of the three models on the entire study period was SWAT-WSVR > SWAT-SVR > SWAT-CUP. The SWAT-WSVR model can predict monthly flow more accurately than the other two models for all sites in the IRW.

This study developed the SWAT-WSVR monthly flow prediction model which was built on the basis of the SWAT and SVR with discrete wavelet transforms and investigated the performance and effectiveness of this model. Precipitation and wavelet components of flow outputted from SWAT-CUP were served as input variables into SWAT-WSVR. The strength of the SWAT-WSVR is its ability to capture the intrinsic nonlinear and non-stationary features between rainfall and runoff while considering physical processes by integrating SWAT. It could be regarded as a compromise method when one cannot directly obtain a desirable accuracy from SWAT-CUP simulation or when one applies SWAT into

a region with limited data available. In future work, more predictors (e.g., temperature, evaporation, relative humidity) will be considered as the model input variables to raise the forecasting accuracy of SWAT-WSVR further and increase its generalization ability.

**Author Contributions:** Conceptualization, L.Y. and K.J.F.; methodology, L.Y.; software, L.Y.; validation, L.Y.; formal analysis, L.Y.; investigation, L.Y.; resources, K.J.F.; data curation, K.J.F.; writing original draft preparation, L.Y.; writing—review and editing, L.Y. and K.J.F.; visualization, L.Y.; supervision, K.J.F.; project administration, K.J.F.; funding acquisition, K.J.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project is funded by the EPA Safe and Sustainable Water Resources National Program.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** https://doi.org/10.23719/1527817.

**Acknowledgments:** The simulation experiments and original draft were completed while Lifeng Yuan worked as an NRC Research Associate at Robert S. Kerr Environmental Research Center, Ada, OK, USA. The final draft was completed while Lifeng Yuan worked at Homeland Security and Materials Management Division, NC, USA. This work does not reflect the views of the U.S. EPA, and no official endorsement should be inferred. It has been reviewed by the Agency but does not necessarily reflect the Agency's views. No official endorsement should be inferred. EPA does not endorse the purchase or sale of any commercial products or services.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Jiashuai Yang <sup>1</sup> , Chaowei Xu 1,\* , Xinran Ni <sup>2</sup> and Xuantong Zhang <sup>3</sup>**


**Abstract:** The imbalance of water supply and demand forces many cities to transfer water across basins, which changes the original "rainfall–runoff" relationship in urban basins. Long-term hydrological simulation of urban basins requires a tool that comprehensively considers the relationship of "rainfall–runoff" and the background of inter-basin water transfer. This paper combines the rainfall–runoff model, the GR3 model, with the background of inter-basin water transfer to simulate the hydrological process of Huangtaiqiao basin (321 km<sup>2</sup> ) in Jinan city, Shandong Province, China for 18 consecutive years with a 1 h time step. Twenty-one flood simulation results of different scales over 18 years were selected for statistical analysis. By comparing the simulation results of the GR3 model and the measured process, the results were verified by multiple evaluation indicators (the Nash–Sutcliffe efficiency coefficient, water relative error, the relative error of flood peak flow, and difference of peak arrival time) at different time scales. It was found that the simulation results of the GR3 model after inter-basin water transfer were considered to be in good agreement with the measured data. This study proves the long-term impact of inter-basin water transfer on rainfall–runoff processes in an urban basin, and the GR3-ibwt model can better simulate the hydrological processes of urban basins, providing a new perspective and method.

**Keywords:** basin water balance; GR3 model; hydrological model; long-term series; single flood process; urban hydrology simulation
