Impact Assessment of Coupling Mode of Hy-Drological Model and Machine Learning Model on Runoff Simulation

1. Introduction

In recent years, global climate change has led to a significant increase in the frequency of droughts and floods worldwide, resulting in severe economic losses [1]. The issue of drought and flood early warning has garnered widespread attention from various coun-tries and sectors of society [2]. Hydrological models provide a mathematical abstraction of the complex hydrological processes[3] in nature. They are widely used for water re-sources management, flood and drought forecasting, and the assessment of eco-hydrological effects [4]. However, accurate streamflow prediction faces numerous challenges due to spatial variations in climate and land surface characteristics[5] and the inherent uncertainty of hydrological models[6]. Notably, in the context of significant climate and land surface changes, the predictive capability of hydrological models based on historical data calibration becomes unstable [7]. Therefore, enhancing the streamflow prediction capability of hydrological models is crucial in the current context.

With the advancement of computational capabilities and the widespread application of big data, many scholars are exploring using data-driven models for streamflow predic-tion [8]. Compared to hydrological models, data-driven models do not require an under-standing of the physical transformation process of precipitation-runoff[9]; they can di-rectly analyze and learn from data, fitting them to the target and offer the advantages of higher simulation accuracy and lower computational complexity. Machine learning (ML) models are a type of data-driven model, and employing machine learning for hydrolog-ical simulation is a current research focus[10], with widespread applications in stream-flow prediction [11]. Standard machine learning models include support vector ma-chines[12], random forests[13], and recurrent neural networks (RNNs)[14]. As hydro-logical events exhibit lag effects[15], most hydrological variables change over time. However, specific machine-learning models are unsuitable for considering the lag ef-fects between hydrological drivers and outputs in streamflow prediction [16]. Some machine learning models incorporating time series data structures, such as RNN mod-els[17], have demonstrated exemplary performance in time series prediction[18]. How-ever, their unique structure also leads to the problem of gradient explosion or disap-pearance [19]. Therefore, to overcome this issue, Hochreiter [20] developed the Long Short-Term Memory (LSTM) neural network based on the RNN structure in 1997. Compared to traditional hydrological models, the LSTM model can better learn the nonlinear relationships in the streamflow process and effectively capture and utilize long-term dependencies in time series[21], thereby improving streamflow prediction capability. However, the LSTM model has multiple parameters, including learning rate and regularization parameters, which must be adjusted according to different watersheds and hydrological conditions. Therefore, selecting appropriate parameters for the LSTM model is crucial in improving its performance. Some scholars have attempted to blend machine learning models with optimization algorithms and achieved promising results. For instance, Yuan et al. [22] combined the LSTM model with Ant Lion Optimizer (LSTM-ALO) for monthly runoff prediction. Adnan et al. [23,24] integrated Relevance Vector Machine (RVM) with Improved Moth-Flame Optimization (IMRFO) model for evaporation forecasting and also applied Firefly Algorithm Particle Swarm Optimization (FFAPSO) to Support Vector Machine (SVM) model for predicting dissolved oxygen concentration in rivers, all yielding favorable outcomes.

Coupling machine learning models in traditional hydrological models can enhance the hydrological model's ability to capture the nonlinear relationship between precipitation and runoff [25], and the hydrological model can provide data input for the machine learning model, thereby improving the runoff modeling capability [26]. Many scholars have researched different aspects of the coupling methods of hydrological and machine learning models. Yang et al. [27] combined a distributed hydrological model with an artificial neural network (ANN) to develop an ML-based hydrological model, finding that this combined model provided satisfactory accuracy for long-term daily runoff sim-ulation. Liu et al. [28] integrated traditional hydrological models with meteorological forecasts and machine learning, reducing the error in runoff prediction. Yu et al.[29], using the Loess Plateau semiarid region as an example, combined the HBV hydrologi-cal model with the LSTM model, proposing two different coupling methods, tight and loose, and found that the hybrid model significantly improved the accuracy and stability of runoff simulation. While many scholars have proposed different model combination methods, there needs to be more research evaluating the performance of different hybrid models and studying the performance of hybrid models in simulating extreme flows, which can confirm whether hybrid models can genuinely contribute to flood and drought early warning. Therefore, our study utilized data from 30 basins in the US CAMELS dataset to construct and evaluate two coupling methods of the hydrological model SIMHYD and the LSTM model. We also compared these two combined models with independent LSTM and SIMHYD models, validating the performance improve-ment of the hybrid models in simulating extreme flows, aiming to analyze whether the combined models enhance the runoff modeling capabilities compared to individual models. The research findings can provide new insights for improving the runoff mod-eling capabilities of traditional hydrological models.

2. Materials and Methods

2.1. Research Area and Data

This study evaluates the predictive capability of hydrological models in 30 watersheds in Washington, USA (Figure 1). Washington is located on the continent's west coast between latitudes 40 and 60 degrees. The western coastal region has a temperate marine climate with moderate temperatures and high precipitation, with an annual average rainfall of over 1000 millimeters. The eastern inland region has a semiarid climate with lower precipitation, with an annual average rainfall of around 250-500 millimeters. Due to the mountainous areas in Washington, winter precipitation often occurs in the form of snow, which melts in spring and becomes a critical water resource for irrigation, urban water supply, and hydropower.

Each watershed includes 25 years (1980-2004) of hydrometeorological data, which are divided into a model training period (1980-1999) and a validation period (2000-2004). The data are sourced from the US CAMLES dataset (https://gdex.ucar.edu/dataset/camels.html), which includes daily sunshine duration, daily precipitation, daily streamflow, daily maximum and minimum temperatures, as well as the elevation, latitude, and atmospheric pressure of the watersheds. Potential evapotranspiration (PET) is an essential input for hydrological models[30], but various methods exist to calculate it. In this study, the FAO-56 Penman-Monteith method is used as a reference, and daily reference crop evapotranspiration is calculated based on data such as daily sunshine duration, daily maximum and minimum temperatures, at-mospheric pressure, elevation, and latitude.

2.2. Models and Methods

2.2.1. Introduction to Hydrological Models

Chiew simplified the HYDROLOG and MODHYDROLOG models and proposed the daily conceptual hydrological model SIMHYD[31]. The advantages of this model are its few input parameters, simple principles and structure, and good accuracy and ap-plicability. The original SIMHYD model does not include a snowmelt module, which may lead to significant errors in simulating snowmelt-dominated watershed models.

The SIMHYD_Snow model incorporates a snowmelt module, Snowmelt[32], into the SIMHYD model (Figure 2). This module divides daily precipitation into rainfall and snowfall based on maximum and minimum temperatures and then calculates daily snowmelt using the degree-day factor method. The SIMHYD_Snow model has 11 pa-rameters, including maximum infiltration capacity, maximum soil moisture content, ac-tual evapotranspiration parameters, and groundwater recharge parameters. This study used a genetic algorithm for parameter calibration to achieve the best fit between the model and observed data, with the root mean square error (RMSE) set as the objective function.

2.2.2. Introduction to LSTM Model

LSTM is a recurrent neural network (RNN) model widely used to classify, process, and predict time series data. In terms of long-term memory, LSTM outperforms RNN due to its unique architecture, which overcomes the limited memory capacity issue in RNN.

The LSTM model consists of recursively connected memory cells, each equipped with a cell state, forget gate, input gate, output gate, and hidden state (Figure 3). The cell state stores and transmits information, while the forget gate is responsible for discarding certain information from the previous time step's cell state. The input gate updates and stores important information in the cell state. The output gate extracts information into the hidden state at the current time step, which serves as the model's output at each time step.

In this study, the LSTM model was constructed using the Deep Learning Toolbox in MATLAB. In order to optimize the structure and performance of the model and obtain the best prediction results, the model parameters were set as shown in Table 1 (The se-lected parameter set demonstrates exemplary performance in terms of computational efficiency and model stability, facilitating practical training and optimization of the model.).

2.2.3. Combined Model Hybrid

In this study, a novel approach is proposed to utilize the output of the SIMHYD model as input for the LSTM model to simulate streamflow. Figure 4 shows primary meteoro-logical data (including potential evapotranspiration PET, precipitation PRE, maximum temperature Tmax, and minimum temperature Tmin) that are input into the SIMHYD model to calculate the runoff. Then, the output of the SIMHYD model is used as input for the LSTM model to obtain the simulated streamflow of the combined model. De-tailed information on the variables is provided in Table 2.

2.2.4. Hybrid Model PE_Trend Based on Dynamic Prediction Effectiveness

2.2.4.1. Modeling Mechanism of Combination Prediction Model Based on Prediction Effectiveness

The combination forecasting model based on predictive validity utilizes the practical information provided by various forecasting methods, calculated in an appropriate weighted average form to obtain the combined model. In order to improve the accuracy of the combined forecasting model, this study focuses on two aspects: firstly, the deter-mination of the weighted average coefficients for each individual forecasting model in the combined forecasting model, and secondly, the examination of whether the com-bined forecasting results are superior to those of the individual forecasting models. The Predictive Validity-based Combined Forecasting Model (PE) selected in this study uti-lizes the mean of predictive accuracy and the mean square deviation reflecting the de-gree of dispersion to construct an optimal calculation model for weighted coefficients through linear programming (Figure 5). It overcomes the influence of the deviation in forecasting results due to the different dimensions of the indicator sequences in the or-dinary combined forecasting model, thereby improving the forecasting accuracy. Moreover, the model is intuitive, computationally concise, and possesses practical ap-plication value.

The measured runoff sequence is {xt, t=1,2,3…N}, and there are m individual prediction models predicting it. x_it is the water quality prediction value of the i-th single-item prediction method at time t, where i ranges from 1 to N. e_it is the relative error between the measured runoff and the predicted runoff by the i-th prediction method at time t, i.e., e_it=（x_i-x_it）/xi. Let Ait=e-|e_it|, and 0≤e_it≤1, then A_it is the prediction accuracy of the i-th runoff prediction model at time t, obviously 0≤A_it≤1.

Let _t be the combined prediction value of x_t, then we have:

ln is the weighted coefficient of a single prediction model and has:

If A_t and e_t are the combined prediction accuracy and relative error at time t, then there are:

The effectiveness of the i-th prediction model is:

Where Q_it represents the weight coefficient of the accuracy A_it of the i-th prediction method at time t on the sample interval and has:

Because prior information is not clear, Q is taken as 1/N. So there are:

The prediction validity sequence values for each individual model. The prediction accuracy sequence of the combined prediction model satisfies:

Variance σ The solution formula for σ²(A) is:

Among ρ_ij is the correlation coefficient between the prediction accuracy of the i-th runoff prediction model and the prediction accuracy of the j-th runoff prediction model, namely:

E(A) represents the average prediction accuracy of the combination prediction method at different times, with larger values leading to higher accuracy: σ (A) The prediction accuracy sequence of the combination prediction model is unstable, and the smaller the numerical value, the better the model.

The effectiveness index of the combination prediction model is defined as:

In the equation, S represents the fitting between the predicted model and the measured runoff. The larger S, the better the linear fitting between the two, and the higher the effectiveness of the predicted results.

The combination prediction model based on prediction effectiveness is:

In this system of equations, as long as the predicted values x_it of m prediction methods at different times are calculated, the definition equation can be used to calculate (E)A_i, (σ)A_i, (σ)²A_i and then to calculate ρ_ij. Using the Lagrange multiplier method, the calculated values of (E)A_i, (σ)A_i, (σ)²A_i, ρ_ij are used to calculate the weighting coefficients l₁, l₂,... for various prediction models, The functions of l_m and l_n are used to find l₁, l₂,... that meet the maximum value requirements of the formula maxS, The weight values of l_m and l_n.

2.2.4.2. A Combined Forecasting Model for Dynamic Prediction Effectiveness

In order to further enhance the prediction accuracy of the combination prediction model based on predictive validity, this study proposes an improved dynamic Predictive Effec-tiveness Trend hybrid model (PE_trend). As shown in Figure 6, we set 5 unit times as one-time step T and input T into the predictive validity hybrid model to calculate the individual model weight values corresponding to this time segment, which are then used to compute the combined model simulation results for the time point after T. (The short-er the set time interval, the higher the accuracy of the model, although it comes at the cost of increased computational time. In our testing, it was observed that setting the time interval to 5 units yielded more favorable simulation results while requiring less compu-tational time.) Subsequently, the initial and final values of period T are incremented by 1, and the combination prediction model is re-input for the calculation to obtain the weighted coefficients of the individual models for the time point after T. This process is iterated until the last time point, thus computing the individual model weight values corresponding to all time points.

2.2.5. Model Evaluation Indicators

This study selects the Nash-Sutcliffe Efficiency (NES) to evaluate model performance. The Nash-Sutcliffe Efficiency is a commonly used index proposed by hydrologist J.R. Nash in 1970. It measures the fit between model simulation results and actual observed values, assessing model simulation performance[33]. A value of NES closer to 1 indicates a better fit of the model to the observed values. NES is calculated using the following formula:

In the equation, Q_m represents the modeled runoff (m³/s), and Q_obs represents the observed runoff (m³/s).

The Taylor diagram is a graphical representation method used to compare the correlation, bias, and standard deviation between model outputs and observed data. It shows the performance of different model predictions in terms of correlation, variance, and standard deviation compared to observed data, helping to evaluate the accuracy and bias of the models. The Taylor diagram is typically presented in a polar plot, with the standard deviation and correlation coefficient of the observed values as axes and the standard deviation and correlation coefficient of the model predictions plotted as points, visually comparing the differences between the models and observed data. A Taylor diagram typically displays three evaluation metrics: Correlation Coefficient, Standard Deviation, and Root Mean Square Error (RMSE).

The Correlation Coefficient measures the strength and direction of the linear relationship between model predictions and observed values. It is calculated as:

In the equation, N is the total number of data points, xi and yi are individual data points for the two variables, x and y are the means of the two variables.

The Standard Deviation measures the variability or spread of data points around the mean. In the Taylor diagram, it represents the dispersion of model predictions and observed values. It is calculated as:

In the equation, N is the total number of data points, xi represents each data point, x is the mean value of the data points.

The Root Mean Square Error is a measure of the differences between values predicted by a model and the observed values. It is calculated as:

In the equation, N is the total number of data points, yi represents the observed values, yi represents the model-predicted values.

3. Results

3.1. Model Runoff Simulation Capability

Figure 7 presents the boxplots of evaluation metrics for the simulated runoff results from the individual models and two hybrid models. The results indicate that, compared to the independent SIMHYD model, the median NES of the validation period for the Hybrid model and the PE_trend model have improved by 16.5% and 21.2%, respectively. That demonstrates a significant enhancement in the runoff modeling capabilities of the two hybrid models, with the PE_trend model exhibiting more robust runoff modeling capabilities. Additionally, the runoff simulation stability of the two hybrid models has significantly improved compared to the SIMHYD model. Furthermore, the Taylor diagrams of the models reveal that the hybrid models and the PE_trend model exhibit higher correlation coefficients, lower root mean square errors, and more minor standard deviations than the SIMHYD model.

3.2. The Ability of the Model to Simulate Extreme Traffic

The key to whether the hybrid model can better forecast flood and drought disasters lies in whether the ability of the hybrid model to simulate extreme flows has been improved compared to that of the individual models [34]. In this study, the top 25% of observed runoff data during the validation period is selected as the relatively high-flow data, while the bottom 25% is selected as the relatively low-flow data for the basin. A boxplot of evaluation metrics for the model's simulation of extreme flows is then plotted (Figure 8).

The models exhibit significant differences in their ability to simulate high flows, with the PE_trend hybrid model demonstrating the best model performance. Compared to the SIMHYD model, the PE_trend model shows a 23.6% increase in the median NES for high-flow simulation. Compared to the LSTM model, the PE_trend model exhibits a 28.4% increase in the median NES of simulated high flows, and the length of the box-plot for this model is shorter, indicating more excellent stability in simulating high flows. In simulating low flows, both the individual models and the hybrid model perform very well. However, compared to the SIMHYD model, the PE_trend hybrid model shows significantly improved stability in simulating low flows. These results indicate that the PE_trend model outperforms the individual models in simulating extreme flows.

4. Discussion

4.1. The Predictive Performance of the Model Under Different Training Periods of Length

The training period length of the runoff prediction model is an important consideration affecting model performance [35]. An excessively long training period may lead to overfitting the model to the training data [36], resulting in poor performance on new data, while a too short training period may prevent the model from fully capturing the data's features, thereby affecting the model's accuracy and stability [37]. Determining the appropriate training period length is crucial for establishing an effective runoff pre-diction model. Therefore, this study analyzed the impact of three different training peri-od lengths on model performance and plotted the evaluation metrics in box plots for training periods of 6 years, 12 years, and 18 years (see Figure 9).

Selecting an appropriate training period length enhances the precision and efficiency of model predictions. As depicted in Figure 9, the predictive performance of the PE_trend model surpasses that of the standalone model in all three scenarios, further highlighting the superior predictive ability of this coupled model in runoff prediction. Comparing the model prediction outcomes for validation period lengths of 6 years and 12 years, it is evident that when the training period is set to 12 years, the performance of all hybrid models strengthens significantly. That indicates that the models needed more data for learning during the 6-year training period. In shorter training periods, models may struggle to capture the intricate patterns and trends within the data. Extending the train-ing period can provide the model with more data to learn from, thereby increasing the accuracy of predictions [38].

Comparing the model prediction results for training period lengths of 12 years and 18 years, the SIMHYD model and the combined prediction model PE_trend exhibit lower NES values when the training period is set to 18 years compared to the 12-year training period. That may be due to the models overfitting the training data, resulting in a de-creased generalization ability on new data [39]. A more extended training period may cause the model to overly rely on noise or specific patterns within the training data, po-tentially overlooking genuine trends and patterns, which could lead to poor predictive performance on new data [40].

In summary, setting an appropriate length for the training period helps the model learn the patterns and trends within the data, thereby enhancing its generalization ability on new data. Additionally, it can reduce the risk of the model overfitting the training data, improve training efficiency, and enable the model to make more accurate predictions.

4.2. The Runoff Simulation Ability of Individual and Combined Models

Compared to individual models, the combined model exhibits more robust runoff simu-lation capability, which confirms Mohammadi's [41] conclusion that combining hydro-logical models and machine learning models can enhance the accuracy of runoff simu-lation. The advantage of the combined model may stem from the complementarity of hydrological and machine learning models. Hydrological models utilize physical princi-ples and parameter estimation based on observed data to simulate runoff processes, of-fering high interpretability and reliability. Machine learning models enhance prediction accuracy by learning associations in data, particularly excelling in handling large amounts of nonlinear data. Combining the two can overcome their respective limitations and achieve a more accurate and robust runoff simulation.

Comparing our research with similar hydrological simulation studies in other regions, we can find some similar and different conclusions. For example, Zarei et al. [42] found that simulating high traffic presents challenges, and even with the use of composite models, the improvement obtained is very limited. This is consistent with our findings that although the composite model PE_trend exhibits good simulation ability for high traffic, underestimation of high traffic still occurs. Sezen and Sraj [43] also acknowl-edged the ability of combined models to simulate low traffic, but their research overes-timated low traffic. When conducting research on selecting training period lengths, Sharma and Kumari [44] found similar conclusions to ours, suggesting that an appro-priate training period length can yield better simulation results. Compared to previous research on runoff simulation combination models, such as the FT-LSTM model pro-posed by Sonali et al. [45] that combines Fourier transform and LSTM models (12.0%), and the Conv TALSTM model proposed by Liu et al. [46] that introduces convolution kernels and attention mechanisms into LSTM networks (5.1%), the PE_trend model has a greater improvement on individual models (12.8%).

Although the combined model performs well in runoff simulation, it still faces chal-lenges and limitations. Specifically, challenges primarily focus on model structure de-sign and prediction accuracy. Firstly, the reasonable combination of hydrological and machine learning models is crucial. According to Lee [47] and others, designing a com-bined model requires considering the integration of the physical principles of hydrolog-ical models and the data-driven capabilities of machine learning models to improve sim-ulation performance. Unreasonable model structure design may lead to reduced perfor-mance. For example, the Hybrid model only inputs primary hydrological data and the output of the SIMHYD model into the LSTM model, which may result in the model's inability to accurately capture the characteristics and patterns of hydrological processes, leading to inaccurate predictions. Secondly, the complexity of model inputs and the is-sue of overfitting. Research by Kavetski [48] and others indicates that machine learning models are prone to overfitting when there are too many input variables, meaning they excessively fit the training data and lose generalization ability. In such cases, the model's performance in simulating extreme flows may need improvement. Therefore, the com-plexity of model inputs may be one of the reasons for the decline in the performance of the Hybrid model during the validation period. In summary, a reasonable model struc-ture design for hybrid models and avoiding overfitting are directions for future im-provement.

In order to enhance the runoff simulation capability of the combined prediction model based on prediction effectiveness, this paper proposes the improved model PE_trend. The model can calculate the weighted coefficients of the individual model for the next step based on the runoff simulation results of the previous step. The results in Figure 7 and Figure 8 indicate that, compared to the individual model, the runoff simulation capability of the combined model PE_trend is more robust. In addition, the PE_trend model demonstrates superior performance in simulating extreme flow events, indicating its ca-pability to capture and forecast extreme flow conditions accurately. That further empha-sizes the robustness and reliability of the improved combination forecasting model in handling extreme flow scenarios, which is crucial for enhancing the overall predictive effectiveness of the model in hydrological applications.When evaluating the model per-formance, we utilized a subset of data from the US CAMELS dataset and did not assess the model's ability to simulate data from other regions. That may result in a slightly lower representation of the model evaluation outcomes. In future work, we plan to evaluate the model performance using data from various regions to enhance the gener-alizability of our findings.

4.3. Limitations and Future Challenges

Data reliability and continuity: The model has strict requirements for the accuracy of input data, and the model requires that the input data should be continuous and real. Continuous data can provide the model with better learning ability, and the model can only demonstrate better simulation ability during the validation period.

The ability to simulate high flow: Hydrological models provide better response to high flow, but still underestimate it. In the future, models for extreme flow need to have faster and more accurate response capabilities.

5. Conclusion

In order to improve the runoff simulation capability of traditional hydrological models, this study combines the machine learning model LSTM with the hydrological model SIMHYD. Two different hybrid models were constructed, and the runoff data of 30 watersheds in Washington State were simulated. The runoff simulation effects of the two model combination methods were evaluated, and further validation of the model's performance in simulating high and low flows was conducted. Based on the experi-mental results, the following main conclusions were drawn:

1. The hybrid model has better runoff simulation capability than traditional hydrological models, significantly improving the median NSE during the validation period.

2. The hybrid model PE_trend performs better in simulating extreme flow conditions than the individual model, contributing to better early warning of floods and droughts.

3. The overall improvement in the hybrid model's performance demonstrates the hybrid model's ability to improve runoff simulation accuracy. Although this study only selected river basins in Washington State, and the results may not be generalized to other basins, the excellent hybrid approach provided can be used as a reference for other regions. However, there are still some issues that need to be addressed in our future research:

(1) Further optimize the model combination method of the hybrid model to improve its learning ability.

(2) Explore the runoff simulation capability of hybrid models in different climatic re-gions.

(3) Enhance the model's ability to simulate high flows and improve its capability to forecast flood disasters.

In order to more intuitively demonstrate the enhancement of the runoff simulation capa-bility of the hybrid model over individual models, this study only selected the LSTM model and the SIMHYD model for combined analysis, which indeed needs a certain level of representativeness. This aspect will be further improved and modified in future work, where updated individual models and more efficient model combination methods will be used for evaluation.