A Hybrid Approach to Improve Flood Forecasting by Combining a Hydrodynamic Flow Model and Artiﬁcial Neural Networks

: Climate change is driving worsening ﬂood events worldwide. In this study, a hybrid approach based on a combination of the optimization of a hydrodynamic model and an error correction modeling that exploit different aspects of the physical system is proposed to improve the forecasting accuracy of ﬂood water levels. In the parameter optimization procedure for the hydrodynamic model, Manning’s roughness coefﬁcients were estimated by considering their spatial distribution and temporal variation in unsteady ﬂow conditions. In the following error correction procedure, the systematic errors of the optimized hydrodynamic model were captured by combining the input variable selection method using partial mutual information (PMI) and artiﬁcial neural networks (ANNs), and therefore, complementary information provided by the data was achieved. The developed ANNs were used to analyze the potential non-linear relationships between the considered state variables and simulation errors to predict systematic errors. To assess the hybrid forecasting approach (hydrodynamic model with an ANN-based error correction model), performances of the hydrodynamic model, two ANN-based water-level forecasting models (ANN1 and ANN2), and the hybrid model were compared. Regarding input candidates, ANN1 considers the historical observations only, and ANN2 considers not only the historical observations that used in ANN1 but also the prescribed boundary conditions required for the hydrodynamic forecast model. As a result, the hybrid model signiﬁcantly improved the forecasting accuracy of ﬂood water levels compared to individual models, which indicates that the hybrid model is able to take advantage of complementary strengths of both the hydrodynamic model and the ANN model. The optimization of the hydrodynamic model allowing spatially and temporally variable parameters estimated water levels with acceptable accuracy. Furthermore, the use of PMI-based input variable selection and optimized ANNs as error correction models for different sites were proven to successfully predict simulation errors in the hydrodynamic model. Hence, the parameter optimization of the hydrodynamic model coupled with error correction modeling for water level forecasting can be used to provide accurate information for ﬂood management.


Introduction
Floods are one of the most severe natural disasters, and more frequent floods occur due to climate change and changes in land use [1][2][3]. River flooding can significantly disrupt important services, including energy, infrastructure, water provision, and transport [4]. Therefore, flood forecasting is a critical issue for the security of individuals and infrastructures and water resources management. An accurate flood simulation and forecasting can provide useful information for flood management and control. A flood routing model can be used to estimate potential damage in a specific reach in which a hydraulic structure may also be included.
Physically based hydraulic models are based on an understanding of the physical process that occurs within a natural system and reproduce the dominant processes through governing equations that are derived under a few simplifying assumptions. Typically, it is not possible to include each process within a simulation model owing to various simplifications and assumptions, and thus, the resulting model can only approximate a real-world system. Despite efforts focused on improving the performance of various open channel flow models, errors between the computed and corresponding observed values are inevitably obtained due to the uncertainties that arise in the model itself, hydrological boundary conditions (upstream flow or lateral inflow) and initial conditions, and model parameters (numerical parameters, geometric parameters, and hydraulic parameters) [5]. All these uncertainties simultaneously affect errors in the model output that produce fairly poor model performance despite a thorough understanding of the governing laws. Gupta et al. [6] emphasized that sophisticated methods are extremely essential to provide corrections for physically based models with respect to water control problems.
In order to provide more accurate information about flood flow, simulation errors produced during the modeling process can be corrected by using postprocessor methods. The selected error correction model should extract useful information about the physical process that was not considered during the construction of a physically based hydrodynamic model such that additionally obtained information can remove systematic biases of the calibrated hydrodynamic model. In contrast to physically based models, data-driven models are not based on either a pre-conceived conceptualization of the behavior of the system or an explicit representation of discrete physical processes [7]. Conversely, in data-driven models, system response is characterized by exploiting statistical information in a set of time-series data [8]. Furthermore, data-driven models can use the advantages of whatever relevant data are available, and this allows such models to represent particular processes [9]. The fact that data-driven models exploit the relationship between data, whereas hydrodynamic flow models are based on physical principles, gives them a complementary nature [7,10,11].
Several studies were performed to improve the river flow forecasting of rainfall runoff models by coupling a conceptual model with a data-driven error correction model [7,[12][13][14][15][16][17]. Konapala et al. [10] used a machine learning algorithm known as a long short-term memory (LSTM) network to simulate the process-based hydrologic model residual. Cho and Kim [18] also employed LSTM to predict the residual errors of WRF-Hydro. The complementary modeling approach was also applied to groundwater problems [19][20][21] and reservoir water level prediction [22]. Nevertheless, there is a paucity of research attention on the use of complementary modeling approaches for hydraulic models, especially in open channel unsteady flow simulations for flood forecasting. Torres-Rua et al. [8] used data-driven models to correct errors in hydraulic simulation models for canal flow schemes in irrigation water management. Bruen and Yang [23] and Preis et al. [24] applied complementary modeling to improve the model performance for urban drainage systems and urban water networks, respectively. Abebe and Price [25] indicated that it is possible to effectively use a parallel data-driven model to predict errors in a flood routing model (a particular form of the 1D diffusive wave equation). The data-driven models used in the aforementioned studies include autoregressive, autoregressive moving average, genetic algorithm, and artificial neural networks (ANNs).
In order to improve the accuracy of flood forecasting, it is significant to integrate or hybridize a physically based hydraulic model and a data-driven model such that it is possible to extract complementary strengths and eliminate weaknesses in the respective modeling methodologies as opposed to continuing to choose between the individual techniques and employing them in isolation [9,26,27]. In this study, ANN models were conjunctively used as error correction models to complement a flood routing model based on the idea that it is possible to significantly improve the overall description of flood flow [25]. A looped network unsteady flow model was constructed for the Han River in South Korea. Unlike Ayvaz [28] who only considered the variability in roughness parameters along a flow reach, the variability in Manning's roughness coefficient, which is related to both the discharge and the flow reach, was considered by adopting a power function as the relationship between Manning's roughness coefficient and discharge in each sub-reach such that the characteristic of the physical process can be fully considered (hereafter, the term "Manning's n" is used instead of "Manning's roughness coefficient"). Specifically, ANNs were employed as error correction models in the study as this modeling approach was applied in several studies and was demonstrated as an optimal technique for developing a model of residual errors for a hydrodynamic model [25]. A nonlinear input variable selection algorithm based on the estimation of partial mutual information (PMI), which was proved as efficient and highly suited for the development of ANN models [29], was applied to ensure that appropriate state variables (that are neither under-specified nor over-specified for estimating simulation errors) are selected as model inputs. The approach integrating the hydrodynamic model and ANNs was estimated by comparing the performances of the hydrodynamic model, two ANN-based forecast models with different input sets, and the hybrid model. The goal of this study is to improve the accuracy of flood water level estimation by coupling the optimized hydrodynamic model with ANN-based error correction modeling, so that different aspects of the physical system from the two types of models are fully exploited.

Study Area
The flow reach in the study corresponds to the 69 km between Paldang Dam and Junryu Gauging Station and is the main stream of the Han River in South Korea ( Figure 1). The water levels are measured at five stations, namely, Paldang Bridge, Banpo Bridge, Hangang Bridge, Heangju Bridge, and Junryu Gauging Station ( Figure 2). Two submerged weirs in the channel, namely, the Jamsil and Singok submerged weirs (Figure 1), are partially gated structures. The gates of movable weirs are fully opened during flood periods, and thus, two different flows occur at the fixed weir and movable weir.
both the discharge and the flow reach, was considered by adopting a power function as the relationship between Manning's roughness coefficient and discharge in each subreach such that the characteristic of the physical process can be fully considered (hereafter, the term "Manning's n" is used instead of "Manning's roughness coefficient"). Specifically, ANNs were employed as error correction models in the study as this modeling approach was applied in several studies and was demonstrated as an optimal technique for developing a model of residual errors for a hydrodynamic model [25]. A nonlinear input variable selection algorithm based on the estimation of partial mutual information (PMI), which was proved as efficient and highly suited for the development of ANN models [29], was applied to ensure that appropriate state variables (that are neither under-specified nor over-specified for estimating simulation errors) are selected as model inputs. The approach integrating the hydrodynamic model and ANNs was estimated by comparing the performances of the hydrodynamic model, two ANN-based forecast models with different input sets, and the hybrid model. The goal of this study is to improve the accuracy of flood water level estimation by coupling the optimized hydrodynamic model with ANN-based error correction modeling, so that different aspects of the physical system from the two types of models are fully exploited.

Study Area
The flow reach in the study corresponds to the 69 km between Paldang Dam and Junryu Gauging Station and is the main stream of the Han River in South Korea ( Figure  1). The water levels are measured at five stations, namely, Paldang Bridge, Banpo Bridge, Hangang Bridge, Heangju Bridge, and Junryu Gauging Station ( Figure 2). Two submerged weirs in the channel, namely, the Jamsil and Singok submerged weirs ( Figure  1), are partially gated structures. The gates of movable weirs are fully opened during flood periods, and thus, two different flows occur at the fixed weir and movable weir.

Hydrodynamic Flow Model
A looped-network unsteady flow model was developed for the Han River. This is owing to the unique features of the Han River with respect to the hydraulic structures in it. Governing equations for the looped-network model consist of equations for the link and the node. With respect to the links, equations are further divided into one for fluvial flow and another for weir-type flow. Governing equations for fluvial links include following continuity and momentum equations: in which Q(x, t) denotes the discharge, h(x, t) denotes the water surface elevation, x and t denote the space and time variables, respectively, A(x, h) represents the cross-sectional area of water, K(x, h) denotes the conveyance, α denotes the energy correction factor assumed to be unity, and g denotes the gravitational acceleration. The conveyance is given as follows: in which R denotes the hydraulic radius, and n denotes Manning's roughness coefficient.
Governing equations for weir-type links are expressed as follows: where s and f denote discharge coefficients, b denotes the effective width of the weir, and hw denotes the weir crest elevation. Upstream and downstream water levels are denoted as hu and hd, respectively, and Qu or Qd denotes the discharge of weir overflow. The nodal continuity equation at any node j is expressed as follows:

Hydrodynamic Flow Model
A looped-network unsteady flow model was developed for the Han River. This is owing to the unique features of the Han River with respect to the hydraulic structures in it. Governing equations for the looped-network model consist of equations for the link and the node. With respect to the links, equations are further divided into one for fluvial flow and another for weir-type flow. Governing equations for fluvial links include following continuity and momentum equations:

∂A ∂t
in which Q(x, t) denotes the discharge, h(x, t) denotes the water surface elevation, x and t denote the space and time variables, respectively, A(x, h) represents the cross-sectional area of water, K(x, h) denotes the conveyance, α denotes the energy correction factor assumed to be unity, and g denotes the gravitational acceleration. The conveyance is given as follows: in which R denotes the hydraulic radius, and n denotes Manning's roughness coefficient.
Governing equations for weir-type links are expressed as follows: where µ s and µ f denote discharge coefficients, b denotes the effective width of the weir, and h w denotes the weir crest elevation. Upstream and downstream water levels are denoted as h u and h d , respectively, and Q u or Q d denotes the discharge of weir overflow. The nodal continuity equation at any node j is expressed as follows: where J denotes the total number of nodes, L j denotes the total number of links connected to node j, and Q j,k denotes the discharge from node k. Q ext (j, t) denotes the external inflow to the node, such as inflows from upstream boundaries or tributaries of the channel. Equal water level is assumed as a nodal energy equation as follows: This implies that the water level h j,k for each computational point on link k is identical to the nodal water level, h j .
Equations (1) and (2) are discretized by the Preissmann's four-point implicit scheme (e.g., [31,32]). The discretized equations in conjunction with Equations (4)-(7) constitute a nonlinear system in which stages and discharges at computational points and nodal water surface elevations are unknown. This system is solved by using the Newton-Raphson method involving the application of the looped double sweep algorithm. The general solution algorithm comprises of four phases for each iteration as follows: link forward sweep, node matrix loading, node solution, and link backward sweep. Details of the solution algorithm are reported in a study by Holly et al. [33].

Model Set Up
The Paldang Dam is the upstream boundary of the model where discharge is known. Water levels are gauged at the Paldang Bridge, Banpo Bridge, Hangang Bridge, and Jeonryu station. Observed water levels at the Paldang Bridge, Banpo Bridge, and Hangang Bridge were used for model calibration and validation, and those at the Jeonryu station were used for downstream boundary conditions. Measured tributary inflows were given as model inputs. The model consists of 10 nodes (1 at either end of the river reach, an upstream/downstream pair at each of the 2 submerged weirs, and 1 at each of the 4 tributary-inflow junction points) and 11 links ( Figure 2). Nodes located upstream and downstream of the submerged weirs are connected by two links while the remaining nodes are connected by one link. There are 156 computational points, and the water levels computed at 3 sites (the Paldang, Banopo, and Hangang Bridges) were used for comparison with the corresponding observed levels.
As previously stated, the submerged weirs at Jamsil and Singok are composed of a fixed weir and a movable weir, and the gates of the movable weirs were fully opened during flood periods. Therefore, channel-and weir-type flows separately occurred on either side. In order to simulate these two flows, a fluvial-and a weir-type link were used for the channel-and weir-type flows, respectively, in the model. The two links connect the upstream and downstream nodes of submerged weirs by forming a loop ( Figure 2) [30].
The initial conditions of each flood event were obtained by modeling 200 h of steady state flow. Boundary conditions at the beginning of an event were used for the inputs of the steady state flow. A steady condition was assumed to be achieved after performing 200 h of steady state flow simulation and obtained flows and water levels were used for the initial conditions of the unsteady state flow model.

Basic Idea
The main procedure used in the study is described as follows: (1) Develop an unsteady flow model that allows for time-varying distributed parameters; (2) Identify model parameters by using an optimization technique; (3) Compute model errors between computed and observed data; (4) Select the variables (at corresponding lag times) that are closely related to the model errors by analyzing the relationship between the model residuals and different state variables (e.g., observed discharge and water level, previous residuals, etc.) at varying lag times by using the PMI-based input variable selection method; Water 2022, 14, 1393 6 of 22 (5) Construct data-driven models (artificial neural networks) that estimate the residuals of the hydrodynamic model by using the best-related variables selected in the previous step; (6) Improve the outputs of the hydrodynamic model by adding the estimated errors as corrections. Figure 3 shows the framework used in the study. The procedure shows that the hydrodynamic flow model is initially executed, and this is followed by the error correction model. The forecasts are improved by combining the results of the two models. The hydrodynamic model used in the study is the looped-network unsteady flow model that considers variable hydraulic roughness. An ANN model was used to estimate residual errors of the hydrodynamic model by considering time series that are closely related to the model errors. The PMI-based approach was used to identify optimal input variables for ANN development based on the dependence between input and output variables. The underlying hypothesis of the proposal is that the data-driven model can capture the systematic patterns presented in the residuals of the hydrodynamic model based on the learning relationship between model residuals and a few data series. It is assumed that the structural problems in the hydrodynamic model and data can be reflected in the detected patterns. (4) Select the variables (at corresponding lag times) that are closely related to the model errors by analyzing the relationship between the model residuals and different state variables (e.g., observed discharge and water level, previous residuals, etc.) at varying lag times by using the PMI-based input variable selection method; (5) Construct data-driven models (artificial neural networks) that estimate the residuals of the hydrodynamic model by using the best-related variables selected in the previous step; (6) Improve the outputs of the hydrodynamic model by adding the estimated errors as corrections. Figure 3 shows the framework used in the study. The procedure shows that the hydrodynamic flow model is initially executed, and this is followed by the error correction model. The forecasts are improved by combining the results of the two models. The hydrodynamic model used in the study is the looped-network unsteady flow model that considers variable hydraulic roughness. An ANN model was used to estimate residual errors of the hydrodynamic model by considering time series that are closely related to the model errors. The PMI-based approach was used to identify optimal input variables for ANN development based on the dependence between input and output variables. The underlying hypothesis of the proposal is that the data-driven model can capture the systematic patterns presented in the residuals of the hydrodynamic model based on the learning relationship between model residuals and a few data series. It is assumed that the structural problems in the hydrodynamic model and data can be reflected in the detected patterns.

Parameter Identification
A hydrodynamic flow model that numerically solves Saint-Venant equations requires Manning's roughness coefficient as a model parameter to represent flow

Parameter Identification
A hydrodynamic flow model that numerically solves Saint-Venant equations requires Manning's roughness coefficient as a model parameter to represent flow resistance. Most previous studies only focused on a lumped-parameter model in which Manning's n is considered as a constant value (e.g., [34][35][36]). Few studies, including Fread and Smith [37] and Ayvaz [28], consider the spatial and/or temporal variation in Manning's n, although it is extremely significant for practical problems. The entire flow reach was divided into two sub-reaches based on the differences in the channel characteristics of the upstream and downstream reaches of the Wangsook stream junction. The downstream reach of the Wangsook stream junction is a channelized reach, while the upstream reach is more similar to a natural river, in which the cross sections are highly irregular. A power function (Equation (8)) that describes the relationship between Manning's n and discharge was adopted for each sub-reach ( Figure 4) [30,38]. Therefore, the model has two different Manning's n-discharge relationships for two sub-reaches that describe the variation in Manning's n with space and discharge, i.e., n = n(x, Q(x, t)), as follows: where I denotes the total number of sub-reaches, α and β denote the coefficients of the power function, and n i denotes Manning's n for sub-reach i. Two power functions are used for the model, and thus it is possible to identify four parameters (α 1 , α 2 , β 1 , β 2 ) in the model calibration. resistance. Most previous studies only focused on a lumped-parameter model in which Manning's n is considered as a constant value (e.g., [34][35][36]). Few studies, including Fread and Smith [37] and Ayvaz [28], consider the spatial and/or temporal variation in Manning's n, although it is extremely significant for practical problems. The entire flow reach was divided into two sub-reaches based on the differences in the channel characteristics of the upstream and downstream reaches of the Wangsook stream junction. The downstream reach of the Wangsook stream junction is a channelized reach, while the upstream reach is more similar to a natural river, in which the cross sections are highly irregular. A power function (Equation (8)) that describes the relationship between Manning's n and discharge was adopted for each sub-reach ( Figure  4) [30,38]. Therefore, the model has two different Manning's n-discharge relationships for two sub-reaches that describe the variation in Manning's n with space and discharge, i.e., n = n(x, Q(x, t)), as follows: where I denotes the total number of sub-reaches, and denote the coefficients of the power function, and denotes Manning's n for sub-reach i. Two power functions are used for the model, and thus it is possible to identify four parameters ( 1 , 2 , 1 , 2 ) in the model calibration. The criterion used for the parameter optimization involves determining model parameters that minimize the sum of squares of the errors between computed and observed stages. The objective function is expressed as follows: (8) in which k t H and k t h denote observed and computed water levels, respectively, at time t at the k-th station, n denotes a parameter vector that represents a set of model parameters, and S denotes the sum of squares of errors between computed and observed water levels.
A parameter optimization tool, PEST, that applies a Gauss-Marquardt-Levenberg algorithm, was used to estimate model parameters [39].

ANN-Based Error Correction Method
Error correction methods mainly focused on correcting the model output time-series of interest [40]. An ANN, which is a data-driven modeling approach, was used to develop an error correction model in this study to reproduce the errors of the hydrodynamic model. Therefore, the output of hydrodynamic model is modified by passing it through the ANN. The criterion used for the parameter optimization involves determining model parameters that minimize the sum of squares of the errors between computed and observed stages. The objective function is expressed as follows: in which H k t and h k t denote observed and computed water levels, respectively, at time t at the k-th station, n denotes a parameter vector that represents a set of model parameters, and S denotes the sum of squares of errors between computed and observed water levels. A parameter optimization tool, PEST, that applies a Gauss-Marquardt-Levenberg algorithm, was used to estimate model parameters [39].

ANN-Based Error Correction Method
Error correction methods mainly focused on correcting the model output time-series of interest [40]. An ANN, which is a data-driven modeling approach, was used to develop an error correction model in this study to reproduce the errors of the hydrodynamic model. Therefore, the output of hydrodynamic model is modified by passing it through the ANN.
An ANN is a powerful tool to address various practical problems, and it is extensively used for simulation and forecasting in diverse areas such as water resources, power generation, finance, and environmental science [41]. Extant studies that focus on its applications in water resources include [23,42]. An ANN learns from data and their potential to accurately describe the behavior of complex nonlinear systems [23]. The ANN model used in this study involves determining the relationship between simulation errors of the hydrodynamic model and a few state variables.
An ANN consists of several layers of neurons (nodes) that are interconnected by links using weights (parameters) [23]. It processes information in a manner similar to the information processes in the human brain [43]. Specifically, ANNs are developed to represent the relationships that are inherent within the data by training the network. A network receives data through input neurons arranged in an input layer and subsequently passes it to the successive layers through a transfer or activation function such as a sigmoid or logistic curve. The flexible architectures enable ANNs to represent complex nonlinear processes without explicitly using a theoretical algorithm to describe the nonlinearity [43].
A three-layer neural network was used in the study to estimate residuals by training with the gradient-based method. The inputs of the ANN were selected in the procedure of input variable selection by using the PMI method, and the output corresponded to the modified residuals. Generally, more than one hidden layer is used for the multi-layer feedforward network. However, the use of a single hidden layer is recommended [44], and this recommendation was adopted by several studies that apply multi-layer feedforward networks in hydrology and water resources [45][46][47][48]. Hence, a hidden layer was used in this study, and the appropriate numbers of hidden neurons were determined by trial and error during the training phase. Figure 5 shows the basic structure of the ANN used in the study. Each neuron in a neural network receives and processes input from other neurons with the exception of neurons at the input layer. The weighted sum of the input data is processed in the hidden layer by using a sigmoid function in this study and is given as follows: where x i denotes the input data at the i-th input neuron; y i denotes the output of the j-th hidden neuron; w ji denotes the weighting value connecting x i and y i ; and L denotes the number of neurons in the input layer.
power generation, finance, and environmental science [41]. Extant studies that focus on its applications in water resources include [23,42]. An ANN learns from data and their potential to accurately describe the behavior of complex nonlinear systems [23]. The ANN model used in this study involves determining the relationship between simulation errors of the hydrodynamic model and a few state variables. An ANN consists of several layers of neurons (nodes) that are interconnected by links using weights (parameters) [23]. It processes information in a manner similar to the information processes in the human brain [43]. Specifically, ANNs are developed to represent the relationships that are inherent within the data by training the network. A network receives data through input neurons arranged in an input layer and subsequently passes it to the successive layers through a transfer or activation function such as a sigmoid or logistic curve. The flexible architectures enable ANNs to represent complex nonlinear processes without explicitly using a theoretical algorithm to describe the nonlinearity [43].
A three-layer neural network was used in the study to estimate residuals by training with the gradient-based method. The inputs of the ANN were selected in the procedure of input variable selection by using the PMI method, and the output corresponded to the modified residuals. Generally, more than one hidden layer is used for the multi-layer feedforward network. However, the use of a single hidden layer is recommended [44], and this recommendation was adopted by several studies that apply multi-layer feedforward networks in hydrology and water resources [45][46][47][48]. Hence, a hidden layer was used in this study, and the appropriate numbers of hidden neurons were determined by trial and error during the training phase. Figure 5 shows the basic structure of the ANN used in the study. Each neuron in a neural network receives and processes input from other neurons with the exception of neurons at the input layer. The weighted sum of the input data is processed in the hidden layer by using a sigmoid function in this study and is given as follows: where denotes the input data at the i-th input neuron; denotes the output of the jth hidden neuron; denotes the weighting value connecting and ; and L denotes the number of neurons in the input layer. The activation function in the output layer is a simple linear function that calculates the weighted sum of the hidden layer outputs that is given as follows: The activation function in the output layer is a simple linear function that calculates the weighted sum of the hidden layer outputs that is given as follows: where z k denotes the output of k-th output neuron; w kj denotes the weighting value connecting y j and z k ; and M denotes the number of neurons in the hidden layer.
In the training phase, the Levenberg-Marquardt algorithm, which is a blend of gradient decent and Gauss-Newton iteration, was used to optimize the weights connecting neurons by minimizing the mean absolute error (MAE) between the targeted data and the  [49,50] describe the ANN model details. The modified residuals were subsequently added to the outputs of the hydrodynamic model to obtain improved results.

PMI-Based Input Variable Selection Method
The first and an important step in ANN development involves identifying the variables that exhibit the strongest relationship with simulation errors and can replicate error behavior. Specifically, the accuracy of the artificial neural network is significantly related to effective input variables. The challenge of selecting appropriate input variables involves selecting the fewest variables that are neither over-specified nor under-specified while adequately describing the underlying input-output relationship [51].
Several input variable selection techniques were developed by ANN modelers. The two main types of input variable selection techniques include wrapper and filter algorithms. Wrapper algorithms are model-based techniques that treat the variable selection as an optimization of the model structure. A weakness of using a wrapper algorithm is the large number of potential computations for trial calibration and evaluation [52]. Furthermore, the applicability of the selected input variables obtained by using a wrapper technique is restricted because the appropriateness of selected variables for a certain model structure is not guaranteed for another structure [53]. In contrast to model-based wrapper techniques, model-free filter techniques use a statistical measure. The separation of the variable selection tasks from the model development tasks makes the technique more effective and also allows the wider applicability of the selected input set to different model architectures [54]. May et al. [54,55] presented a review of recent input variable selection methods for ANNs.
An efficient technique reduces the need for extensive trial-and error and arbitrary judgement during model development. The PMI technique is a model-free filter approach that considers inter-dependencies between candidates and deals with nonlinear modeling without the assumption of linearity. An input variable selection algorithm based on the estimation of PMI was initially introduced by Sharma [56] and further developed by Bowden et al. [29] and May et al. [54,55]. Some studies indicated that PMI-based input variable selection is a highly suitable measure for input variable selection during ANN development since the ANN development implicitly assumes non-linear relationships between variables [43,54,55,57].
Input variable selection using PMI-based algorithms can significantly reduce the amount of trial-and error during ANN model development and improve the performance of an ANN model. The PMI-based selection algorithm is briefly described as follows [43,54,55]: (1) Determine potential candidates based on previous knowledge; (2) Initialize the selected input variable set S that corresponds to a null set since no input variables are selected initially, and the candidate set C; (3) Compute the residual of output (u) and that of candidate variables (v) by subtracting relationship with currently selected inputs as follows: where, Y denotes the output variable; X and Z denote the selected and unselected input variables, respectively, up to the current step;m Y (X) denotes the non-parametric regression estimator by applying a Kernel function that is given as follows: where o denotes the number of observed values (y i , x i ); E[y|X = x ] denotes the conditional expectation of y given x; K b denotes the Gaussian kernel that is given as follows: Here, d denotes the dimension of X, b denotes the kernel bandwidth (or smoothing parameter), and σ denotes the sample covariance matrix.
The kernel bandwidth b is as follows: where sf denotes the scaling factor that controls the degree of kernel smoothing; b r denotes the Gaussian reference bandwidth that is given as follows: Harrold et al. [58] empirically determined that a bandwidth of approximately 1.5 b r provided stable estimates of mutual information; (4) Estimate the PMI I(v; u) and determine candidate ∈ C that maximizes PMI as follows: where f (v i ) and f (u i ) denote the marginal probability densities, and f (v i , u i ) denotes the joint probability density.
The estimator for f at a given location x is given as follows: The kernel density estimation (KDE) that is one of non-parametric density estimation techniques were used in step (3) and (4) since it is typically considered as accurate and suitably robust [55]; (5) Compute the values of the considered termination criterion.
The termination criterion is a key consideration in the input variable selection as it significantly affects computational efficiency and selection accuracy [59]. The Akaike information criterion (AIC) [60] can be used to measure the trade-off between the size of the input set and the accuracy of the regression filter [54]. The AIC is given as follows: where o and p denote the number of observations and model parameters, respectively, and r i represents residuals. During the input variable selection, the AIC initially decreases with the increase in the number of parameters due to a reduction in the residuals and reaches a minimum value. Subsequently, the AIC increases as the 2p term penalizes the selection of additional variables. Therefore, the selection procedure is terminated when the AIC reaches its minimum point; (6) Based on the values of the considered termination criterion, a decision is adopted to either select a candidate and continue or terminate the selection.

Parameter Identification for Hydrodynamic Model
The hourly observed data of 12 flood events at Paldang Bridge, Banpo Bridge, and Hangang Bridge were selected for the calibration and validation of the model. The peak discharges of these flood events at the upstream boundary range from 5147 m 3 /s to 16,146 m 3 /s. The flood events were numbered in the decreasing order of the peak discharge at the Paldang Dam and split into two data sets. Flood events corresponding to 1, 3, 5, 7, 9, and 11 were used for calibration, and the rest were used for model validation (Table 1). The model parameters were determined by averaging the optimized parameters for each flood event in the calibration. The estimated power functions in Figure 6 that were derived by considering the geometric average of calibrated coefficients show that Manning's n decreases when the discharge increases in both sub-reaches and that the values of Manning's n for the upstream reach tend to be greater than those for the downstream reach given a uniform discharge. The identified power functions obtained in the calibration were used to simulate flood events corresponding to 2, 4, 6, 8, 10, and 12, as a validation procedure. The average root mean square errors (RMSEs) for calibration and validation are 0.205 and 0.308, respectively, and this indicates that the model that uses a power function as the Manning's n-discharge relationship can estimate flood water levels with acceptable accuracy.
During the input variable selection, the AIC initially decreases with the increase in the number of parameters due to a reduction in the residuals and reaches a minimum value. Subsequently, the AIC increases as the 2p term penalizes the selection of additional variables. Therefore, the selection procedure is terminated when the AIC reaches its minimum point; (6) Based on the values of the considered termination criterion, a decision is adopted to either select a candidate and continue or terminate the selection.

Parameter Identification for Hydrodynamic Model
The hourly observed data of 12 flood events at Paldang Bridge, Banpo Bridge, and Hangang Bridge were selected for the calibration and validation of the model. The peak discharges of these flood events at the upstream boundary range from 5,147 m 3 /s to 16,146 m 3 /s. The flood events were numbered in the decreasing order of the peak discharge at the Paldang Dam and split into two data sets. Flood events corresponding to 1, 3, 5, 7, 9, and 11 were used for calibration, and the rest were used for model validation (Table 1). The model parameters were determined by averaging the optimized parameters for each flood event in the calibration. The estimated power functions in Figure 6 that were derived by considering the geometric average of calibrated coefficients show that Manning's n decreases when the discharge increases in both sub-reaches and that the values of Manning's n for the upstream reach tend to be greater than those for the downstream reach given a uniform discharge. The identified power functions obtained in the calibration were used to simulate flood events corresponding to 2, 4, 6, 8, 10, and 12, as a validation procedure. The average root mean square errors (RMSEs) for calibration and validation are 0.205 and 0.308, respectively, and this indicates that the model that uses a power function as the Manning's n-discharge relationship can estimate flood water levels with acceptable accuracy.

PMI-Based Input Variable Selection
In order to test the approach that combines the hydrodynamic model and ANNs, a performance comparison of the hydrodynamic model, ANNs, and hybrid model (hydrodynamic model with ANN-based error correction model) was carried out. An ANN model, labeled here as ANN1, was developed to forecast water levels at specific locations using historical observations as inputs. The considered input candidates and output of ANN1 are shown in Figure 7. Furthermore, an ANN model, denoted as ANN2, was constructed for water level forecasting with input candidates including not only the inputs used in ANN1 but also prescribed boundary conditions at t + 1, t + 2, . . . , t + l that were required for the hydrodynamic forecast model (Figure 8). Finally, an ANN for error correction of the hydrodynamic model, labeled as ANN3, was developed with candidate inputs including the predicted boundary conditions of the hydrodynamic model, historical observations, calculated water levels by the hydrodynamic model, and recorded errors of the hydrodynamic model ( Figure 9). In Figures 7-9, Q D and Q D represent the observed and prescribed discharges at the Paldang Dam, respectively; Q Trib and Q Trib are the observed and prescribed tributary inflows, respectively; h J and h J are the observed and prescribed water levels at the Junryu gauge station; h Sta is the observed water level at gauging stations; HM_h Sta is the computed water level by the hydrodynamic model; ∂h Sta ∂t is the time derivative of the observed water level; HM_ ∂h Sta ∂t is the time derivative of the computed water level by the hydrodynamic model; h k is the forecasted water level at the k-th station; and Err k is the errors of the hydrodynamic model at the k-th station.
ANN model, labeled here as ANN1, was developed to forecast water le locations using historical observations as inputs. The considered input c output of ANN1 are shown in Figure 7. Furthermore, an ANN model, den was constructed for water level forecasting with input candidates includin inputs used in ANN1 but also prescribed boundary conditions at + 1, that were required for the hydrodynamic forecast model (Figure 8). Finall error correction of the hydrodynamic model, labeled as ANN3, was d candidate inputs including the predicted boundary conditions of the model, historical observations, calculated water levels by the hydrodynam recorded errors of the hydrodynamic model (Figure 9). In    performance comparison of the hydrodynamic model, ANNs, and hybrid mo (hydrodynamic model with ANN-based error correction model) was carried out. ANN model, labeled here as ANN1, was developed to forecast water levels at spec locations using historical observations as inputs. The considered input candidates a output of ANN1 are shown in Figure 7. Furthermore, an ANN model, denoted as ANN was constructed for water level forecasting with input candidates including not only inputs used in ANN1 but also prescribed boundary conditions at + 1, + 2, …, that were required for the hydrodynamic forecast model (Figure 8). Finally, an ANN error correction of the hydrodynamic model, labeled as ANN3, was developed w candidate inputs including the predicted boundary conditions of the hydrodynam model, historical observations, calculated water levels by the hydrodynamic model, a recorded errors of the hydrodynamic model ( Figure 9). In Figures 7-9, D and represent the observed and prescribed discharges at the Paldang Dam, respectively; T and Trib ′ are the observed and prescribed tributary inflows, respectively; ℎ J and ℎ J ′ the observed and prescribed water levels at the Junryu gauge station; ℎ Sta is the observ water level at gauging stations; HM_ ℎ Sta is the computed water level by hydrodynamic model; ℎ Sta is the time derivative of the observed water level; HM_ ℎ is the time derivative of the computed water level by the hydrodynamic model; ℎ ′ is forecasted water level at the k-th station; and Err is the errors of the hydrodynam model at the k-th station.   For the considered ANNs, the number of potential inputs can be quite large, shown in Figures 7-9; however, given their correlated nature, many may be redunda The travel time of flow from upstream boundary (Paldang Dam) to downstre boundary (Junryu gauging station) approximately corresponds to 1~2 h for most histori flood events. However, it is difficult to precisely measure the travel time from upstream boundary to the downstream boundary owing to the tidal effect. Therefo lagged variables from current time t to the previous 5 h (the maximum lag time d) w assumed as sufficient in the study.
The derivative of water levels can be used to describe the flow situation. For examp a zero or low derivative indicates normal or base flow while a high positive derivat indicates the rising limb of a flood event, and a high negative derivative indicates recession limb of a flood event. Therefore, the rate of flow changed and its lagged variab were also considered as possible input variables in ANN models.
From the input candidates, input variables for ANN models were selected applying the PMI technique. Data of 12 flood events (Table 1) were used for input varia selection using PMI. Based on the AIC criterion, input sets that satisfy AIC min w selected. The selected input variables for ANNs at different lead times and at differ sites are listed in Tables 2-4. In the tables, subscribes P, B, H, and W represent the Palda Bridge, Banpo Bridge, Hangang Bridge, and Wangsook Stream, respectively.
Selected input variables for ANN1 show that water level forecasting at the Palda Bridge is mainly determined by observed water levels at the Paldang Bridge, a discharges at the Paldang Dam. The inflow from the Wangsook Stream also affects water level forecasting as the lead time increases (Table 2). Table 2 shows that 1-h le time forecasting at the Banpo Bridge and Hangang Bridge is related to the observed wa levels at four gauge stations, i.e., Paldang Bridge, Banpo Bridge, Hangang Bridge, a Junryu gauge station. Table 3 shows the selected input variables for ANN2 wh considers not only historical observations but also the predicted boundary conditio used in the hydrodynamic model. The selected variables show that water level forecast at the Paldang Bridge is determined by the latest observations at the Paldang Bridge a predicted discharges at the Paldang Dam. For the 1-h lead time forecasting at the Ban Bridge and Hangang Bridge, observations at four gauge stations are selected as inp variables. As the lead time increases, water levels at the Banpo Bridge and Hanga Bridge are forecasted by including additional predicted data as inputs such as water lev at the Junryu station and discharges from the Paldang Dam and Wangsook Stream. Ta For the considered ANNs, the number of potential inputs can be quite large, as shown in Figures 7-9; however, given their correlated nature, many may be redundant. The travel time of flow from upstream boundary (Paldang Dam) to downstream boundary (Junryu gauging station) approximately corresponds to 1~2 h for most historical flood events. However, it is difficult to precisely measure the travel time from the upstream boundary to the downstream boundary owing to the tidal effect. Therefore, lagged variables from current time t to the previous 5 h (the maximum lag time d) were assumed as sufficient in the study.
The derivative of water levels can be used to describe the flow situation. For example, a zero or low derivative indicates normal or base flow while a high positive derivative indicates the rising limb of a flood event, and a high negative derivative indicates the recession limb of a flood event. Therefore, the rate of flow changed and its lagged variables were also considered as possible input variables in ANN models.
From the input candidates, input variables for ANN models were selected by applying the PMI technique. Data of 12 flood events (Table 1) were used for input variable selection using PMI. Based on the AIC criterion, input sets that satisfy AIC min were selected. The selected input variables for ANNs at different lead times and at different sites are listed in Tables 2-4. In the tables, subscribes P, B, H, and W represent the Paldang Bridge, Banpo Bridge, Hangang Bridge, and Wangsook Stream, respectively.
Selected input variables for ANN1 show that water level forecasting at the Paldang Bridge is mainly determined by observed water levels at the Paldang Bridge, and discharges at the Paldang Dam. The inflow from the Wangsook Stream also affects the water level forecasting as the lead time increases (Table 2). Table 2 shows that 1-h lead time forecasting at the Banpo Bridge and Hangang Bridge is related to the observed water levels at four gauge stations, i.e., Paldang Bridge, Banpo Bridge, Hangang Bridge, and Junryu gauge station. Table 3 shows the selected input variables for ANN2 which considers not only historical observations but also the predicted boundary conditions used in the hydrodynamic model. The selected variables show that water level forecasting at the Paldang Bridge is determined by the latest observations at the Paldang Bridge and predicted discharges at the Paldang Dam. For the 1-h lead time forecasting at the Banpo Bridge and Hangang Bridge, observations at four gauge stations are selected as input variables. As the lead time increases, water levels at the Banpo Bridge and Hangang Bridge are forecasted by including additional predicted data as inputs such as water levels at the Junryu station and discharges from the Paldang Dam and Wangsook Stream. Table 4 shows that all the selected input sets for the error correction models (i.e., ANN3) at different locations and different lead times include Err(t), suggesting that Err(t) significantly contributes to the forecasting of simulation errors. Several studies indicated that Err(t) includes most of the information on Err(t + T) (e.g., [7]). Therefore, a good error prediction can be obtained by considering the Err(t) as the only input variable of the model. However, Torres Rua [61] performed internal tests and indicated that the performance of the single variable is less robust than the variable combination for an error correction model. Table 4 shows that 1-h lead time forecasting of errors for the Paldang Bridge is mainly related to Err p (t), while 1-h lead time forecasting of errors for the Banpo Bridge or Hangang Bridge is determined by the latest errors, time derivatives of forecasted water levels at their own locations, and the latest observed water levels at the Junryu gauge station. Table 4 also shows that relatively high numbers of input variables were selected for the 3-h lead time forecasting of errors at the Banpo Bridge and Hangang Bridge, which shows that the error forecasts generated by an ANN model are somewhat inexplicable, and this is a shortcoming of the ANN modeling approach. May et al. [55] stated that the lack of interpretability is not surprising because the PMI variable selection method for the ANN development is somewhat holistic and it does not consider the contribution of individual input variables to the model.

Application of ANN Models
The data from six flood events that were used for the calibration of the hydrodynamic model were used for the ANN training, and six events that were used for the validation of the hydrodynamic model were used for the ANN validation (Table 1). Therefore, approximately 55% and 45% of the data were used for training and validation, respectively.
ANNs for water level forecasting at different sites were constructed based on the selected input variables and the number of hidden neurons, as determined by trial and error. The appropriate number of hidden neurons was determined by testing the training data set with hidden nodes ranging from 1 to 10. Therefore, an ANN was optimized by selecting the optimal number of hidden neurons given the selected inputs. The ANN that exhibited optimal performance in the validation was selected as the optimal network. The network architectures for different ANNs are listed in Table 5. Given an ANN architecture, the weights and biases of the network were determined in the model training. Figures 10-12 show the observed and predicted stage hydrographs at three sites at different lead times for flood event 2. As shown in the figures, the combination of the hydrodynamic model and ANNs provides a highly accurate prediction of flood water levels and adequately captures the behavior of water flow. Even with respect to the 2-h or 3-h lead time forecasting, the hybrid approach results in an estimation that agrees well with the observations.    Two statistics, i.e., the Nash-Sutcliffe efficiency coefficient (NSE) and RMSE, wer used to evaluate the accuracy of different forecasting models. The values of two statistic for training periods are summarized in Table 6, and the best results for each locatio according to each performance metric are highlighted in bold. NSE and RMSE indicat that the hydrodynamic model alone has relatively low accuracy compared to othe forecasting models. Forecasting the accuracy of ANN1 decreases with the longer lead tim of forecasting, which is naturally expected as ANN1 only uses the historical observation to forecast water levels. Forecasts by the hybrid approach are most accurate at the Paldan Bridge, while ANN2 also shows good performance at the Banpo Bridge and Hangan Bridge. Two statistics, i.e., the Nash-Sutcliffe efficiency coefficient (NSE) and RMSE, were used to evaluate the accuracy of different forecasting models. The values of two statistics for training periods are summarized in Table 6, and the best results for each location according to each performance metric are highlighted in bold. NSE and RMSE indicate that the hy-drodynamic model alone has relatively low accuracy compared to other forecasting models. Forecasting the accuracy of ANN1 decreases with the longer lead time of forecasting, which is naturally expected as ANN1 only uses the historical observations to forecast water levels. Forecasts by the hybrid approach are most accurate at the Paldang Bridge, while ANN2 also shows good performance at the Banpo Bridge and Hangang Bridge. The validation results are shown in Table 7 and the best results are also highlighted in bold. It is immediately obvious from the table that the hybrid approach is the most accurate based on the RMSE and NSE when applied to all three locations, and the accuracy reduces with the increasing lead time of forecasting. ANN2 showed better accuracy in the model training at the Banpo Bridge and Hangang Bridge, while it shows worse performance than the hybrid approach in the model validation. The prediction accuracy of ANN1 decreases as lead time increases, and it performs worse than the hydrodynamic model for 3-h lead time forecasting at the Hangang Bridge. The reduction in the simulation errors of the hydrodynamic model is evident after error correction by ANNs, which implies that coupling a hydrodynamic model and an ANN (as opposed to only using a hydrodynamic model or an ANN) results in highly accurate predictions of water levels. For 1-h lead time forecasting, the hybrid approach achieved 70%, 80%, and 76% reduction in the RMSE at the Paldang Bridge, Banpo Bridge, and Hangang Bridge, respectively, when compared with the hydrodynamic model. In addition, the reduction in the RMSE by applying the hybrid approach are 18%, 13%, and 33% at 3 locations when compared with the ANN2, which indicates the better performance of a hybrid model compared to an ANN model. For the 2-h lead time forecasting, the hybrid model achieved 65%, 69%, and 65% reduction in the RMSE at the Paldang, Banpo and Hangang Bridges, respectively, when compared with the hydrodynamic model, and the reduction is 26%, 16%, and 33% when compared to the ANN2. ANN1 shows less accuracy than the hybrid model and ANN2. The NSE for the hybrid model is over 0.99 when the model is applied for the 1-h or 2-h forecasting at all locations, and it is over 0.98 when it is applied to the 3-h forecasting. The results of this research indicate that the hybrid model performs better than a data-driven model (i.e., ANN1 or ANN2), and a data-driven model has better performances when compared to the physically based model (i.e., hydrodynamic model), which is consistent with the research by Cho and Kim [18].

Conclusions
This study presents the combined application of a physically based hydrodynamic model and a data-driven model for the improvement of flood water level prediction. In order to enhance the performance of the hydrodynamic model, calibration was initially performed by using power functions that describe the relationships between Manning' s n and discharge for two sub-reaches. The optimized hydrodynamic model performed well when compared to the observed data at gauging stations. Subsequently, the impact of the systematic errors in the optimized hydrodynamic model was significantly reduced by using a complementary data-driven error correction model.
To test the performance accuracy of the hybrid model which combines the hydrodynamic model and ANNs, a comparison of the hydrodynamic model, two water level forecasting models based on ANNs (i.e., ANN1 and ANN2), and the hybrid model (HM+ANN3) was carried out. It was demonstrated that the water level forecasted by the hydrodynamic model is generally less accurate than those by two ANN models (i.e., ANN1 and ANN2). ANN1 is a data-driven forecasting model, which was developed based on the historical observations. Since the model only uses the historical data, the accuracy of the model decreases as the lead time of forecasting increases. As the hydrodynamic forecast model need predicted boundary conditions as inputs, ANN2 also included the predicted boundary conditions as potential inputs. By providing the same types of data for the hybrid model and ANN2, it is reasonable to test if a hybrid model performs better than an ANN. Though the same data were provided for constructing these two models, the outputs of the hydrodynamic model, e.g., computed water levels and their time derivatives, and historical errors between computed and observed water levels were additionally included in the input candidates of the error correction model (ANN3).
A PMI-based input variable selection method was employed to identify suitable inputs to ANN models. The selected variables for predicting water levels or simulating errors can play a crucial role in theoretically describing physical mechanisms. The ANNs developed based on the PMI-based input variable selection provide efficient and accurate predictions by mapping specific relationships between the inputs and output. ANN-based error correction approach for a hydrodynamic model presented in this study is simple to implement and provides the uncertainty reduction in water level forecasting. This method can implicitly account for a range of uncertainties that are contained in the deterministic model output time series derived offline prior to application.
The comparison of the model performances showed that the hybrid model performs better than the individual hydrodynamic model and ANNs in terms of RMSE and NSE. Though ANN2 is more accurate than the hybrid model during the training period at the Banpo Bridge, it does not show better performance for validation data. The predictions produced by the hybrid model can capture the variation pattern of the forecasted water level by describing the rising and falling tendencies of the water levels in a fairly accurate manner. The results of this study indicate that aggregating a physically based hydrodynamic model and a data-driven model rather than only applying a hydrodynamic model or a data-driven model can fully exploit different aspects of the physical system to minimize the prediction errors between estimated and observed water levels. The findings of this study suggest that the combination of a hydrodynamic flow model and an ANN can be applied for water level forecasting with high accuracy. It should be noted that this study focused on error correction in the forecasting mode, given the assumption that the hydrodynamic model is executed in the forecasting mode with the same accuracy as in the simulation mode. The accuracy of water level forecasting can be different depending on the accuracy of the predicted boundary conditions of the hydrodynamic model.