Abstract
Dynamically dimensioned search (DDS) algorithm is a new-type heuristic algorithm which was originally developed by Tolson and Shoemaker in 2007. In this study, the DDS algorithm is applied to automate the calibration process of an unsteady river flow model in the Tamsui River basin, which was developed by Wu et al. (2007). Data observed during 2012 and 2013 are collected in this study. They are divided into three groups, one for the test case, one for calibration, and one for the validation. To prove that the DDS algorithm is capable of solving this research problem and the convergence property, a test simulation is first performed. In the studied area, the whole river systems are divided into 20 reaches, and each reach has two parameters ( and ) to be determined. These two parameters represent resistance coefficients for low- and high-water conditions. Comparing with another algorithm, it is shown that the DDS algorithm has not only improved on the efficiency but also increased the stability of calibrated results.
1. Introduction
The river resistance coefficient, also known as the value in the Gauckler–Manning formula, is an important hydraulic parameter used in flow calculations for planning and design in hydraulic engineering. Proper hydraulic parameters, coupled with appropriate flood simulation models, can simulate water level variations in the calculation for basins during typhoon floods and thereby help mitigate the potential loss of life and property. Nowadays, one-dimensional (1D) and two-dimensional (2D) hydraulic models seem to be crucial tools for the development of flood propagation and for supporting flood risk assessment [1–3]. Whether or not river resistance coefficients are correctly estimated, they affect the result of water level calculation, planning, design, and operation of channels, and accuracy of various hydraulic calculations including water transport efficiency. Since different watersheds have different channel characteristics, each reach of a river has its own resistance coefficients, which is generally computed by an experienced person.
Strickler [4] used an empirical formula with of the riverbed to estimate the value of . Cowan [5] proposed a formula to estimate the river resistance coefficient based on the material of the river bottom, windiness of the channel, distribution of plantation, and man-made structures. In addition, textbooks on canal water mechanics offer their own methods for estimating river resistance under various canal conditions, for example, Chow [6] and Henderson [7]. Yen [8] indicated that, because of vegetation and other obstacles, the local roughness of a typical composite channel cross section will be higher or lower for the floodplains due to different bed texture. Wilson [9] highlighted that Manning’s coefficient decreases with increasing flow depth, reaching an asymptotic constant at relatively high flow depths. Vegetation was also pointed out by other references [10–12] that will affect the water level. However, these types of estimates are subjective in nature and thus are difficult to pass on.
The trial and error method has also been used to adjust the river resistance coefficients so as to match the water levels with the observed values as closely as possible. However, calibration using this process is usually very time-consuming. Additionally, there is no guarantee that the resistance coefficients calibrated in this manner will reach the optimal values and are thus unable to pass the validation test. Therefore, research on automatic calibration has been conducted recently. Unlike traditional, labor-intensive trial and error methods, automated calibration is a systematic optimization process. The parameter to be calibrated is fed into the model, and then the numerically simulated value is compared to empirical measurements to determine the right direction for parameter adjustment. The iterative process is executed until the results converge.
Becker and Yeh [13] used the influence coefficient algorithm in an attempt to calibrate the river resistance coefficients and hydraulic radius index for a one reach river. The same researchers expanded their method to calibrate river resistance coefficients and hydraulic radius indexes for multiple reach rivers in 1973 [14].
Lal [15] investigated the use of singular value decomposition (SVD) to calibrate the river resistance coefficients for the upper reaches of the Niagara River. It was suggested that an underdetermined parameter validation problem (the number of parameters to be calibrated exceeds the number of water level or flow rate measurement stations) could be simplified by grouping parameters, where all parameters in a given group had the same value. However, the resistance coefficients calculated using different numbers of groups were not alike, thus suggesting that the method was not stable.
Currently, automatic calibration methods can be divided into two types: gradient search methods and heuristic algorithms. Although gradient search methods are based on a theoretical foundation rigorously, they are not suitable for solving multidimensional optimization problems, for two reasons: the solution space is limited to the neighborhood of the initial solution; if providing a good initial solution is not possible, the calculation mechanism is trapped by the local optimal solution in most cases; thus, the search mechanism lacks diversity; the gradient search method can only handle problems where the gradient exists. However, there are many real problems in nature that contain discontinuities. This further limits the possibility of finding a solution with the gradient method. For instance, Yang [16] investigated the use of the conjugate gradients method to calibrate the resistance coefficients of the Shimen Dazhen Channel and showed that the method was only suitable for simple river reaches.
Thus, Heuristic algorithm methods were developed to overcome these deficiencies. With the help of artificial intelligence (AI) technologies, these methods are mainly adopted effectively and efficiently to find an optimal solution within a feasible solution space. Researchers have developed various algorithms and calculation mechanisms by imitating various examples of natural intelligence.
Nowadays, heuristic algorithms are widely employed to solve optimization problems in many fields, such as electronics, IT, engineering, and economics [17]. The most well-known heuristic algorithms currently in use include tabu search, genetic algorithms, simulated annealing algorithms, neural networks, and fuzzy methods [18–21]. The most beneficial feature of heuristic algorithms is that they can find good solutions quickly. They also have built-in mechanisms to avoid falling into local optimum solution traps, such as the tabu list mechanism in the tabu search and the mutation operand in the genetic algorithm. Simulated annealing algorithms, in particular, escape local optimal solutions by adjusting the probability of the solution moving to neighboring regions based on temperature changes. These search methods are nothing more than intelligent trial and error methods. They combine several natural laws, the ability to learn, probability characteristics, fuzzy concepts, memory functions, and so on to construct calculation methods that are most capable of solving optimization problems.
Chan [22] and Huang [23] both used unsteady-flow simulation models for basin-wide river systems coupled with a real-value-coding of genetic algorithm or an additional simulated annealing algorithm (SARvcNGA, Simulated Annealing Real-value-coding Niche Genetic Algorithm). The root-mean-square error (RMSE) between the calculated water level and the measured water level was used as the objective function to ascertain the river resistance coefficients that could best describe the actual flow conditions. Though SARvcNGA has been proven to be more efficient and accurate than genetic algorithm with a real-value-coding, it still requires lots of time to complete the computations. Clearly, there is room for improvement in terms of efficiency, especially in support of flood mitigation efforts, where timely information is paramount.
The dynamically dimensioned search algorithm (DDS), which is a new type of heuristic algorithm, was developed by Tolson and Shoemaker [24]. It was designed for calibration problems with many parameters. Many studies related to this method have been conducted, such as [25–27]. Tu [28] used DDS algorithm in automatic calibration of an unsteady river flow model.
Present study focuses on optimizing the river resistance coefficients by DDS algorithm. It is the extension of the previous research efforts [28]. Because one of the branches of study area was dredged for flood control in 2010, the coefficients in the reaches need to be renewed; thus the flood forecasting model can keep its accuracy.
These resistant coefficients will be used in combination with an unsteady-flow simulation model for basin-wide river systems [29] to simulate water levels in flood research. The results will also be compared with those by using the algorithm proposed in Huang [23]. It is illustrated that the DDS algorithm will provide an accurate result and is a diverse and robust method for automatic calibration of the unsteady river flow model.
2. Dynamically Dimensioned Search Algorithm
This section is mainly drawn from “dynamically dimensioned search algorithm for computationally efficient watershed model calibration” [24], for completeness.
2.1. Overview
The dynamically dimensioned search algorithm (DDS) is mainly designed to solve multiparameter calibration problems. When applied to optimization problems requiring large amounts of computational time, this search method does not demand complicated parameter adjustment. That is, there are no additional control parameters in the algorithm that require separate adjustment, such as the initial temperature, temperature reduction factor in the simulated annealing algorithm, mutation probability, or number of ethnic groups in the genetic algorithm. Thus, the number of uncertainty factors decreases, and the greatest possible portion of the calculation time is used for searching the optimal solution. Other features of this algorithm include the option to limit the number of searches to be performed to suit the user’s time limit. The range of feasible solutions to be searched can also be adjusted based on the maximum number of searches to be performed. These design features do not affect the diversity or robustness of the search.
2.2. Calculation Mechanism
One special characteristic of the DDS algorithm is that, regardless of the number of searches set by the user, the search for candidate solutions proceeds from global to local. This is made possible by an adjustment mechanism that decreases the number of algorithm parameters to be determined based on the changes in the probability. Moreover, a new candidate solution is created only by modifying the best current solution. Both solutions are substituted into the objective function for comparison. The candidate solution only replaces the best current solution if it decreases the value of the objective function. Otherwise, it is rejected, and another candidate solution is created from the best current solution. This process is repeated until convergence conditions are achieved or the maximum number of searches has been performed.
The DDS algorithm includes three main calculation steps: setting up the initial solution; selecting a candidate solution; and updating the best current solution. Following are detailed descriptions of each of these steps.
(1) Setting up the Initial Solution. Consider that the model being calibrated contains parameters (called decision variables). The only three algorithm parameters to set in the DDS algorithm which must be initialized before are as follows: the neighborhood perturbation size parameter ; the maximum number of searches ; and the upper and lower limits of the parameters to be determined (). Next, a solution set can be chosen at random from the feasible solution space to serve as the initial solution for the calibration process. And the initial solution is substituted into the objective function as , see (1), to complete the setup.
(2) Selecting a Candidate Solution. During the search for the optimal solution, the solution currently being searched and modified is called the “current best solution.” The current best solution is modified according to calculation mechanism to search for the next candidate solution. During this process, all of the solutions that could possibly become the next candidate solution are known as “neighboring solutions.” The set of all neighboring solutions is called the “neighborhood” . The solution that is actually chosen from the neighborhood is called the “candidate solution.” The heuristic algorithm methods are distinguished from each other by comparing the manner through which they simultaneously satisfy both diversity and robustness while selecting the candidate solution. The following are the rules used by the DDS algorithms for selecting a candidate solution.
Totally there are parameters to be determined. Among them, parameters are selected in neighborhood, , based on the probability as shown in (2), where is the number of the current iterations. As the number of iterations increases, the number of chosen parameters decreases. This is equivalent to the gradual progression from a global to a local search. This continues until the convergence conditions are satisfied or the maximum number of searches has been performed.where is the maximum numbers of searches.
For example, when the number of parameters is 10 (i.e., ) and the maximum number of searches is 1000 (i.e., ), the probability of the first search () can be calculated as
Therefore, during the first search, the probability that each parameter to be determined will be chosen for updating is 100%, and so forth.
(3) Updating the Best Current Solution. The parameters chosen in the previous step are updated according towhere and is standard normal random variable.
If the value of any updated parameter is less than the lower limit or greater than the upper limit, it is adjusted according to
Further, the objective function is evaluated and the best current solution is updated if necessary, as shown in
3. Problem Description and Research Zone
3.1. Overview of Tamsui River Basin
This study attempts to ascertain the river resistance coefficients for the unsteady-flow model of the Tamsui River basin. The Tamsui River basin is located at the northern extremity of Taiwan. It has an area of 2,726 km2 and a main stream length of 158.7 km, making it the third largest river in Taiwan. The Tamsui River has three principal tributaries: the Dahan Stream in the south and the Sindian Stream in the middle, both converging at Jiangzihcuei near Taipei City; and in the north, the other tributary Keelung River converges with the Tamsui River at Guandu.
The Dahan Stream is 135 km long, and its watershed has an area of 1,163 km2 and an average gradient of 1/37. The upper reaches are the Shimen Reservoir catchment, which has an area of 759 km2. There are mountain valleys in the upper reaches of the Dahan Stream and terraces and alluvial plains in the middle and lower reaches. Further, there is well-developed transportation infrastructure within the Dahan Stream watershed boundaries. It is located in the core region of the greater Taipei metropolis.
The Sindian Stream is 82 km long, and its watershed has an area of 910 km2, 89% of which is mountain slopes. The river’s abundant water flow is an important source of water for Taipei.
The Keelung River originates at Jingtong Mountain in Pingci Village of Taipei County. Its main stream is 89.4 km long and its watershed has an area of 490.77 km2. The upper reaches of the Keelung River are between the Jieshou Bridge in Houdong and the Dahua Bridge in Cidu. Here, the average gradient is approximately 1/250. The middle reaches, from the Dahua Bridge to the Nanhu Bridge, have an average gradient of around 1/4900. The lower reaches are from the Nanhu Bridge to the mouth of the river, with an average gradient of around 1/6700. The riverbed within the basin is fairly flat and windy (data source: the 10th River Management Office, Water Resource Agency, Ministry of Economic Affairs, Taiwan).
3.2. Study Area
This research focuses on the Tamsui River basin because of its rather comprehensive profile data and the availability of several water level measurement stations for data collection. The research zone includes the Keelung River downstream of the Dahua Bridge, the Dahan Stream downstream of the Sinhai Bridge, the Sindian Stream downstream of the Zhongzheng Bridge, and the Tamsui River from the mouth of the river and up, as shown in Figure 1.
3.3. Construction of the River Hydraulic Model
The hydraulic analysis method used in this research is an unsteady-flow simulation model for basin-wide river systems, CCCMMOC (complex-compound channel flow modeling by the multimode method of characteristics) [30–32]. This method mainly uses the method of characteristics to solve Saint Venant’s equations. It is a practical, stable, accurate hydraulic simulation model that can simultaneously handle tidal flows, floods, and rapidly varied flows due to natural or man-made circumstances and be used in compound-complex channel systems [33–35].
The multimode method of characteristics of the second kind (MMOC-II) has been chosen as the numerical scheme of solution. This is a powerful numerical scheme belonging to the method of characteristics (MOC); it merges the temporal reach-back, spatial reach-back, spatial reach-out, temporal reach-out (only at the boundaries), and classical schemes into one large scheme. The scheme operates in the explicit mode, that is, for all the interior nodes it computes the unknown variables at one grid point at a time. It is thus deemed exceptionally suitable for the basin-wide unsteady-flow computation [36].
In an explicit form each explicit scheme is limited by the respect of the Courant number which is limited up to 1. The time increment Δt is set to 15 (s) for the stability. Because we roundly set Δ as 500 (m) and the water depth is about 5~20 (m), the Courant number is in the range of 0.21~0.42. On the other hand, the numerical scheme does not have shock-capturing properties so it could not deal with discontinuous phenomena like hydraulic jumps, backwater effects, and so on.
Schaffranek and Lai [37] indicated that if the effects of nonhomogeneous terms (other than the frictional-resistance term) are taken into consideration, it is possible to perceive that these terms would also affect the model performance. Thus, if a numerical model is calibrated by neglecting certain nonhomogeneous terms that should be present in the prototype function, the frictional-resistance term will be compelled to absorb any discrepancy, so that “resistance coefficient” does not merely mean “roughness factor/coefficient” in the unsteady-flow model.
Considering the highly variable stage while the flood must rise and return, the value in this study is programmed to vary with depth, , besides its ordinary variation along the reach; that is, . This variation in depth mainly emerges from the change in the cross section (e.g., the main channel and overbank section), and the change in vegetation and bank texture along the depth.
Generally, river hydraulic models deal with steady flows, which can execute calculation effectively by setting the resistance coefficients as fixed values. However, as water levels rise and fall at the natural state, flow profiles and surface roughness change as well. In order to truly grasp the water level, , changes and simulate entire flood processes to ensure that the model is just as accurate for high-water levels as well as low-water levels, the river resistance coefficients used must be automatically calibrated within the range of predicted water levels. This adjustment mechanism is shown inwhere is the average gradient of the river resistance coefficient over the range of allowable water levels; and are the allowable upper and lower water level limit in each reach of the river and can be measured from on-site observation, respectively; river resistance coefficients and correspond to water levels and , respectively [29, 34]. The calibration of the Gauckler–Manning coefficient, , is the objective of this research. In the following section, water level data from a severe typhoon flood event are used to calibrate . After calibration, water level data from two more typhoon flood events are used to verify the applicability of the calibration.
Because most regions in the research area are mountains and in the urban regions the three branches are all confined by the embankments, 1D modeling is enough for flood forecasting usage. One could find more explanation in [29] about why the 1D model could be employed for the unsteady-flow simulation in basin-wide river systems.
3.4. Setup of the Parameter Calibration Optimization Model
The river system for which this study intends to carry out parameter optimization is the entire Tamsui River basin. The research zone includes the Keelung River (downstream of the Dahua Bridge), the Dahan Stream (downstream of the Sinhai Bridge), the Sindian Stream (downstream of the Zhongzheng Bridge), and the Tamsui River (down to the mouth). In all, there are 130 segments (a segment is a length of channel between two profiles). According to channel characteristics, the research zone was separated into 20 reaches. Reaches 1 through 7 belong to the main stream—Dahan Stream and Tamsui River; Reach 8 belongs to the Sindian Stream tributary; and Reaches 9–20 belong to the Keelung River tributary system. See Figure 2.
3.4.1. Determining the Number of Parameters
Each reach has two parameters that need to be determined, and . In all, there are 40 parameters. According to the manner by which this method was initially set up, these parameters should have been calibrated for each segment. However, since there were far fewer measurement stations than the parameters to be determined (there were a total of 20 water level measurement stations in the river system), the parameters for each reach were determined instead of that for each segment. This lowered the uncertainty error due to insufficient measurement data.
3.4.2. Objective Function
The objective function of this optimization method minimizes the RMSE between the calculated water level and the measured water level, as shown inHere, and are the observed and calculated water level at station at time , respectively; is the total number of time steps; and is the number of water level measurement stations. The main goal in drawing up the objective function is to match the modeled water level hydrograph with measured hydrograph, as closely as possible.
3.4.3. Constraints
This optimization method has two constraints, namely, the basin-wide river system unsteady-flow model and the upper and lower limits of the river resistance coefficient. The constraints are described below.
(1) The Basin-Wide River System Unsteady-Flow Model. This study uses an unsteady-flow simulation model for basin-wide river systems [29–32] to carry out hydraulic calculations. In this method, mathematical models based on Saint Venant’s equations are transformed into characteristic equations. Then, the second multimode characteristic method (MMOC-II) gives the numerical solution.
(2) Upper and Lower Limits during Parameter Validation. Given that this study investigates the calibration of real river resistance coefficients, it is impossible to know the exact range of the parameter values in advance. If the range is too large, the feasible solution space will be too big, which will affect convergence time. If the range is too small, it might be impossible to find a suitable solution. However, Cowan [5] offered an equation to estimate the river resistance coefficients using the effective resistance of several factors, including channel bottom material, the windiness of the channel, and the distribution of plant life and man-made structures. Considering that the Dahan Stream, Sindian Stream, and Tamsui River channels within the research zone are wide and have mostly been dredged, their range of river resistance coefficients was set between 0.020 and 0.040. For the Keelung River, resistance coefficient ranges for each reach were set between 0.020 and 0.065 in accordance with the research study conducted by Lai and Tsay [35].
3.5. Automatic Parameter Calibration Procedure
In this study, automatic parameter calibration is mainly accomplished by using the DDS algorithm to repeatedly run the unsteady-flow simulation model for basin-wide river systems. First, an initial solution set for the parameters to be determined is created from within the range of feasible solutions. This solution is substituted into the unsteady-flow simulation model for basin-wide river systems to obtain the simulated water level. Then, the RMSE between the simulated water level and the measured water level is calculated. Finally, the simulated result is inputted into the DDS algorithm to produce a new candidate solution. This process is repeated until convergence conditions are satisfied or the preset maximum number of searches is reached. The objective is to minimize the value of the objective function.
4. River Resistance Coefficient Automatic Calibration and Validation Results
This section is divided into three parts, one for test case, one for calibration, and the other for validation. The corresponding data observed during 2012 and 2013 are collected in this study and are also divided into three groups: one for the test case, one for calibration, and one for the validation.
4.1. Dynamically Dimensioned Search Algorithm: Test Case
This study utilized actual water level measurements to carry out parameter optimization; therefore, the following factors may influence the optimization results.
(1) Water Level Measurement Error. Any errors in the measured data will affect the value of the objective function and thereby alter the results of the optimization process.
(2) Number and Location of Reaches. Considering the complicated nature of the current research zone comprising natural rivers, the parameters are set up individually for each reach. In this way, fewer parameters can be used to express the special characteristics of each reach, and, therefore, the optimization process can be greatly simplified. However, different methods of distinguishing between reaches will lead to discrepancies in the optimization results.
In order to prove that the DDS algorithm is capable of solving this research problem and the convergence property, a test simulation is first performed. This section of the paper describes the test simulation in which all of the potential uncertainty factors are removed from the to-be-determined parameter optimization process, with the objective that a basin-wide optimal solution can easily be found.
4.1.1. Test Case Setup
In this test case, the process can be divided into two stages. In the first step, since the actual river resistance coefficients in the Tamsui River basin cannot be known, a set of default values for (low-water level resistance coefficients) and another set for (high-water level resistance coefficients) are created. These default sets are then assumed to be the set of actual river resistance coefficients and are substituted into the unsteady-flow simulation model for basin-wide river systems to complete a water level simulation. In step , the model parameters used in step are assumed to be unknown, and the DDS algorithm is used in combination with the model’s water level output from step to carry out parameter calibration. The objective is to evaluate whether or not the results of the automatic parameter calibration process converge to values close to the original sets of assumed river resistance coefficients from step .
The water levels measured at the four boundary water level stations, as well as the water levels calculated for each reach using the default river resistance coefficients in the model, are substituted into the river hydraulic model for simulation. There are resistance coefficients and that need to be calibrated. The data points of Soulik typhoon event from 7/11/2013 15:00 to 7/14/2013 12:00 were included in the test case, and the time step was 20 minutes. Therefore, each water level station had 208 water level data points. After that, the DDS algorithm was used to produce a candidate solution at every iteration. Stop condition for calculation was a maximum iteration step of 40,000.
4.1.2. Test Results
Test results show that the objective function begins to converge after 29,367 model iterations. A 2.29 GHz CPU with 4 GB of RAM and the program Matlab 7.10 with an execution file of Fortran are used. The calculation time was approximately 17 seconds per iteration, so convergence occurred after 139 hours. The minimum convergence error value was 0.05 m. Figure 3 (low-water level) and Figure 4 (high-water level) show the comparison of the default parameter values with the calibrated parameter values. In addition, the default parameters and calibrated parameters were separately substituted into the unsteady-flow model to calculate the water level at each hypothetical water station. Figure 5 shows plots comparing the calculated water level hydrographs at the hypothetical water level stations in the model. These plots show that the method resulted in a very good water level fit at each station.
In summary, it is clear that default resistance coefficient values and calibrated resistance coefficient values from the test case were not exactly the same; however, they had very similar trends and distributions. The modeled water level outputs also supported the positive results of this test. One possible reason that the test case did not result in the same parameter values was that the huge size of the feasible solution space lengthened the time required during latter stages of the search. Nevertheless, this did not affect the high efficiency of the early stages of the search. This section shows that if the size of the feasible solution space can be appropriately reduced, and sources of uncertainty and error can be effectively controlled in future research, then the DDS algorithm is sufficiently robust to solve the problem posed in this study. It is not unreasonable to expect that this method will lead to an optimal basin-wide solution.
4.2. Results of Parameter Calibration Using Actual Measured Water Levels
In the previous section, an independently designed test case was set up. The purpose was to determine whether or not the methods proposed in this study had the ability to solve the river resistance coefficient optimization problems. In this section, water measurements from an actual flood event are used to calibrate the river resistance coefficients with the DDS algorithm. The distribution of river reaches and water measurement stations in the research zone is shown in Figure 2. River resistance coefficients and are calibrated using water level data from Typhoon Soala (7/31/2012~8/3/2012). Besides, in order to demonstrate the advantages of the DDS algorithm in terms of convergence results and computation time, the same parameter calibration will also be carried out using a Simulated Annealing Real-value-coding Niche Genetic Algorithm (SARvcNGA) conducted by Huang [23] with the same computer hardware.
Typhoon Soala was the 9th typhoon announced by Taiwan’s Central Weather Bureau in 2012. It was announced at around 8 p.m., Taipei time. Its path was similar to that of Typhoon Morakot in 2009. The Central Weather Bureau issued the first maritime typhoon alert for eastern sea area of Taiwan at 20:30 on July 30th. At 20:30 on July 31st, during the 8th typhoon alert, both land and sea alerts were issued. The Central Weather Bureau rescinded both the land and maritime portion of the typhoon alert at 14:30 on August 3rd. A total of 31 alerts were issued during the course of the storm. At its strongest, Typhoon Soala had a minimum atmospheric pressure of 960 kPa at its center, and a level-7 storm wind radius of 220 km as measured by the Central Weather Bureau. According to Central Emergency Operation Center, at least five people were killed and two were missing in Taiwan, in addition to 16 injured. Agricultural losses across the island were estimated at NT$ 812.51 million (US$ 27.14 million) by August 6th.
The time block used to calibrate hydraulic parameters was 7/31/2012 15:00–8/3/2012 12:00. The time step was 20 minutes; therefore, there were 208 data points from every water level measurement station. The water level measurement stations in the Dahan Stream watershed and along the main course of the Tamsui River were Sinhai Bridge, Gueiyang, Liuguan, Taipei Bridge, Dihua, Tudigungbi, and river mouth stations. In the Sindian Stream watershed, the water level stations were Zhongzheng Bridge and Shuangyuan. In the Keelung River watershed, the water level stations were Dahua Bridge, Wudu, Jiangbei Bridge, Shehou Bridge, Nanhu Bridge, Chengmei Bridge, Nanjing, Yangguang Bridge, Beian, Dajhih Bridge, Jianguo, Chengde Bridge, Jiantan, Shezihdao, and Jhoumei. Among these, the Sinhai Bridge, river mouth, Zhongzheng Bridge, and Dahua Bridge stations were boundary-point water level stations in the model. Data from the 20 stations were included in the objective function calculations.
During each iteration, the objective function is calculated and recorded. The convergence process of the two algorithms, DDS and SARvcNGA, during automatic calibration is shown in Figure 6. For DDS algorithm, stop condition for calculation was a maximum iteration step of 20,000. The objective function began to converge after approximately 10,223 iterations. Minimum convergence error was 0.237 m. The time required for each iteration was about 17 seconds. The model ran for around 2 days using a 2.29 GHz CPU with 4 GB of RAM. The resistance coefficient values calibrated for each reach are arranged in Table 1. Figure 7 compares the measured water levels with water levels calculated from calibrated river resistance coefficients at some water level station.
And for SARvcNGA, the objective function began to converge after the 43rd generation. The model ran 11,673 times. The minimum convergence error was 0.271 m. The model ran for approximately 7 days. By comparing the convergence results, it is clear that DDS algorithm used in this paper is both more efficient and more accurate than the previously proposed method, SARvcNGA.
4.3. Verifying Calibration Results with Measured Water Level Data
After calibrating the river resistance coefficients, it is necessary to verify that the calibrated parameters can accurately reflect the true roughness values of the channels. In this section, the calibrated parameters with the DDS algorithm (and SARvcNGA) are substituted into the basin-wide unsteady river flow model again to calculate simulated water levels during different typhoon floods. Below are the results of the two floods used for validation.
4.3.1. Typhoon Trami
Data from the 8/20/2013 15:00–8/23/2013 12:00 time block were entered into the objective function. The time step was 20 minutes; therefore, each water measurement station had 208 data points. The water levels calculated by the model were compared to measured water levels. For DDS algorithm, the calculated RMSE was 0.286 m. (For SARvcNGA, the calculated RMSE was 0.376 m.) Figure 8 contain plots comparing measured water levels with water levels calculated using the automatically calibrated river resistance coefficients.
4.3.2. Typhoon Usagi
Data from 9/20/2013 15:00–9/23/2013 12:00 were substituted into the objective function. The time step was 20 minutes; therefore, there were 208 data points from each water measurement station. Water levels calculated by the model were then compared to measured water levels. For DDS algorithm, the calculated RMSE was 0.284 m, which is very similar to the previous validation result. This shows that the calibration results can be applied to different typhoon flood events with excellent and robust performance. (The calculated RMSE was 0.349 m for SARvcNGA.) Figure 9 contain plots comparing measured water levels with water levels calculated using the calibrated resistance coefficients at some water measurement stations.
4.4. Discussion
Comparing the results of the test case with actual calibration results, we found that although actual calibration results fit reasonably well, the test case results fit even better. It is evident that, under real conditions, parameter optimization is affected not only by measurement data, but also by some unexpected factors. One could not determine the bed resistant coefficient by looking merely at the vegetation and bottom grain size. The unexpected factors cause the difficulty of determining the bed resistant coefficient.
Table 2 lists the Froude number of every junction during peak of Typhoon Saola. We have to admit that weirs or bridges across the river may affect the water levels [38–41] and the model we use in this study cannot deal with supercritical flow. But since there is no weir in this study area and no supercritical flow in the domain, too detailed supercritical flows around local points become insignificant. As we mentioned previously, what we are calibrating in this study is the “resistance coefficients” that include bridged and other uncertain effects, and it is one of the most difficult parts; we need a better methodology.
By using the DDS algorithm, the resistance coefficients are adjusted by letting the computed water levels match the observed ones as close as possible. Therefore, the flow-resistance coefficients used in this study have been considered as a combination of several factors influencing the resistance and could represent the characteristics of the reaches. This roughly satisfied the whole basin-wide area.
Since flow in the reach area is confined to a relatively narrow valley, we suppose that using 1D modeling is enough. The results presented here should, however, be treated with caution, when extending of the methodology to other reaches and flood events. According to the previous references, they may result in several aspects, shown below:(1)A 1D approach has been used for flood propagation in this study. The 1D modeling shows several limitations and drawbacks due to its intrinsic inadequacy for the description of rivers with wide floodplains and high sinuosity [1, 42, 43]. Some chosen cross sections are shown in Figure 10. If the study area is with a wide floodplain of complex morphology, it may therefore require a 2D approach.(2)Momentum transfer between main channel and floodplain was not taken into account in particular. Several improvements have been recently introduced in the 1D modeling both in steady [44, 45] and unsteady [39, 46] situations to take into account momentum transfer between main channel and floodplain. Since the 1D model we used shows some limitations, momentum flux is difficult to approximate by simply using 1D variables such as water level and area of a cross section. We will try another model which has considered more previous published works in the future. Two-dimensional modeling may also be more appropriate if other hydraulic processes, such as turbulent momentum and exchange between channel and floodplain waters, have to be taken into consideration.(3)The value in this study is programmed to vary with depth, , besides its ordinary variation along the reach. This variation in depth mainly emerges from the change in the cross section, and the change in vegetation and bank texture along the depth. So we have and , which is large depending on different bed texture in a composite cross section.
(a)
(b)
(c)
5. Conclusions
Significant breakthroughs have already been made using the DDS algorithm to automatically quickly optimize river resistance coefficients without sacrificing accuracy. Below are the conclusions of this study:(1)The DDS algorithm is superior to other heuristic algorithms such as SARvcNGA because it does not have control parameters that require calibration, such as initial temperature and temperature reduction factors of the simulated annealing algorithm, or the mutation probability and the number of ethnic groups of the genetic algorithm. The DDS algorithm has simple steps and features that are conducive for easy operation. Moreover, it can minimize the effect that different users have on the final solutions.(2)This study set up a parameter optimization test case that imitated the real research zone. The goal was to eliminate measurement uncertainties. The test results proved that the DDS algorithm could actually seek out acceptable hypothetical river resistance coefficient values in a very large feasible solution space, and that the simulated water level histories fit well with the hypothetical water level histories.(3)Based on the assumptions and limitations shown in the above section, a simplified 1D model has been used although the effect of vegetation, momentum exchange, and sediment transport was not considered in particular to limit the work. Difficulties in obtaining the true resistance coefficient values can be expected. However, through careful analysis of the research zone and simplification of the parameters to be optimized, it is still possible to reasonably and effectively model the actual flow conditions. It is hoped that terms representing the momentum exchange between the main channel and the floodplain added to 1D unsteady-flow modeling will be used in the future to get more accurate result. Vegetation, man-made structures, and even sediment transport considered in hydraulic model also may bring more accurate result. Moreover, two-dimensional models, capable of providing accurate simulations of the hydraulic processes occurring in the floodplains, also may be employed to improve the results although they could be computationally expensive and require topographic data which were difficult to gather until recently.(4)This research utilized a 2.29 GHz CPU with 4 GB of RAM. It was capable of completing the parameter optimization in approximately 2 days. In comparison, the method proposed by Huang [23] required 7 days to converge using the same hardware. This demonstrates that the DDS algorithm is not only a powerful optimization tool, but also an efficient one.
Competing Interests
The authors declare that they have no competing interests.
Acknowledgments
The authors are grateful for the financial support from the Ministry of Science and Technology, Taiwan (Grant no. 105-2221-E-415-008). The authors would also like to thank the 10th River Management Office, Water Resources Agency, and the Hydraulic Engineering Office, Public Works Department, Taipei City Government, Taiwan, for providing data.