Printer Friendly

Automatic Calibration of an Unsteady River Flow Model by Using Dynamically Dimensioned Search Algorithm.

1. Introduction

The river resistance coefficient, also known as the n 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 d50 of the riverbed to estimate the value of n. 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: (1) 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; (2) 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: (1) setting up the initial solution; (2) selecting a candidate solution; and (3) 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 D parameters (called decision variables). The only three algorithm parameters to set in the DDS algorithm which must be initialized before are as follows: (1) the neighborhood perturbation size parameter r = 0.2; (2) the maximum number of searches m; and (3) the upper and lower limits of the parameters to be determined ([x.sup.min], [x.sup.max]). 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 F, see (1), to complete the setup.

[mathematical expression not reproducible]. (1)

(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" {N}. 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 D parameters to be determined. Among them, J parameters are selected in neighborhood, {N}, based on the probability P(i) as shown in (2), where i 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.

P(i) = 1 - ln(i)/ln(m), (2)

where m is the maximum numbers of searches.

For example, when the number of parameters is 10 (i.e., D = 10) and the maximum number of searches is 1000 (i.e., m = 1000), the probability of the first search (i = 1) can be calculated as

P(1) = 1 - ln(i)/ln(1000). (3)

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 to

[mathematical expression not reproducible], (4)

where [[sigma].sub.j] = r([x.sup.max.sub.j] - [x.sup.min]) and N(0,1) 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

[mathematical expression not reproducible]. (5)

Further, the objective function F([]) is evaluated and the best current solution is updated if necessary, as shown in

If F([])[less than or equal to][]

[] = F([]), then [] = []. (6)

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 [km.sup.2] 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 [km.sup.2] and an average gradient of 1/37. The upper reaches are the Shimen Reservoir catchment, which has an area of 759 [km.sup.2]. 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 [km.sup.2], 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 [km.sup.2]. 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 (1) simultaneously handle tidal flows, floods, and rapidly varied flows due to natural or manmade circumstances and (2) 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 At is set to 15 (s) for the stability. Because we roundly set Ax 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 n value in this study is programmed to vary with depth, z, besides its ordinary variation along the reach; that is, n = n(x, z) = n(x, z(x, t)). This variation in depth mainly emerges from the change in the cross section (e.g., the main channel and overbanksection), 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, z, 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 in

n = [n.sub.d] z < [z.sub.d]

n = [n.sub.d] + s x (z - [z.sub.d]), [z.sub.d] [less than or equal to] z [less than or equal to] [z.sub.u] (7)

n = [n.sub.u] z > [z.sub.u]

s = ([n.sub.u] - [n.sub.d])/([z.sub.u] - [z.sub.d]), (8)

where s is the average gradient of the river resistance coefficient over the range of allowable water levels; [z.sub.u] and [z.sub.d] 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 nu and [n.sub.d] correspond to water levels [z.sub.u] and [z.sub.d], respectively [29, 34]. The calibration of the Gauckler-Manning coefficient, n, is the objective of this research. In the following section, water level data from a severe typhoon flood event are used to calibrate n. 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, [n.sub.d] and [n.sub.u]. 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 in

min[1/T 1/T[[T.summation over (t=1)][I.summation over (i=1)][[[z.sup.o.sub.i,t] - [z.sup.c.sub.i,t]].sup.2]].sup.1/2]. (9)

Here, [z.sup.o.sub.i,t] and [z.sup.c.sub.i,t] are the observed and calculated water level at station i at time t, respectively; T is the total number of time steps; and I 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 (MMOCII) 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 [n.sub.d] (low-water level resistance coefficients) and another set for [n.sub.u] (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 (2), the model parameters used in step (1) are assumed to be unknown, and the DDS algorithm is used in combination with the model's water level output from step (1) 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 (1).

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 [n.sub.d] and [n.sub.u] 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 [n.sub.d] and [n.sub.u] 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/201215:00-8/3/201212: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 n value in this study is programmed to vary with depth, z, 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 [n.sub.d] and [n.sub.u], which is large depending on different bed texture in a composite cross section.

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.


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.


[1] M. S. Horritt and P. D. Bates, "Evaluation of 1D and 2D numerical models for predicting river flood inundation," Journal of Hydrology, vol. 268, no. 1, pp. 87-99, 2002.

[2] H. Apel, G. T. Aronica, H. Kreibich, and A. H. Thieken, "Flood risk analyses--how detailed do we need to be?" Natural Hazards, vol. 49, no. 1, pp. 79-98, 2009.

[3] N. M. Hunter, P. D. Bates, M. S. Horritt, and M. D. Wilson, "Simple spatially-distributed models for predicting flood inundation: a review," Geomorphology, vol. 90, no. 3-4, pp. 208-225, 2007.

[4] A. Strickler, Beitraege zur Frage der Geschwindigkeitsformel und der Rauhigkeitszahlen fuer Stroeme, Kanaele und geschlossene Leitungen, vol. 16, Mitteilungen des Eidgenoessischer Amtes fuer Wasserwirtschaft, Bern, Switzerland, 1923 (German).

[5] W. L. Cowan, "Estimating hydraulic roughness coefficients," Agricultural Engineering, vol. 37, no. 7, pp. 473-475, 1956.

[6] V. T. Chow, Open-Channel Hydraulics, Mcgraw-Hill, New York, NY, USA, 1959.

[7] F. M. Henderson, Open-Channel Flow, The Macmillan, New York, NY, USA, 1966.

[8] B. C. Yen, "Hydraulic resistance in open channels," in Channel Flow: Centennial of Manning's Formula, B. Yen, Ed., Water Resource Publications, Littleton, Colo, USA, 1991.

[9] C. A. M. E. Wilson, "Flow resistance models for flexible submerged vegetation," Journal of Hydrology, vol. 342, no. 3-4, pp. 213-222, 2007.

[10] F. Huthoff, D. C. M. Augustijn, and S. J. M. H. Hulscher, "Analytical solution of the depth-averaged flow velocity in case of submerged rigid cylindrical vegetation," Water Resources Research, vol. 43, no. 6, Article ID W06413, 2007

[11] J. C. Curran and W. C. Hession, "Vegetative impacts on hydraulics and sediment processes across the fluvial system," Journal of Hydrology, vol. 505, pp. 364-376, 2013.

[12] C. Wang, S.-S. Zheng, P.-F. Wang, and J. Hou, "Interactions between vegetation, water flow and sediment transport: a review," Journal of Hydrodynamics Series B, vol. 27, no. 1, pp. 24-37, 2015.

[13] L. Becker and W. W. Yeh, "Identification of parameters in unsteady open channel flows," Water Resources Research, vol. 8, no. 4, pp. 956-965, 1972.

[14] L. Becker and W. W. Yeh, "Identification of multiple reach channel parameters," Water Resources Research, vol. 9, no. 2, pp. 326-335, 1973.

[15] A. M. W. Lal, "Calibration of riverbed roughness," Journal of Hydraulic Engineering, vol. 121, no. 9, pp. 664-671, 1995.

[16] T. H. Yang, Calibration of Manning's n in an open channel with trapezoid sections [M.S. thesis], National Taiwan University, Taipei, Taiwan, 2001 (Chinese).

[17] D. F. Jones, S. K. Mirrazavi, and M. Tamiz, "Multi-objective meta-heuristics: an overview of the current state-of-the-art," European Journal of Operational Research, vol. 137, no. 1, pp. 1-9, 2002.

[18] B. Wan and W. James, "SWMM calibration using genetic algorithms," in Proceedings of the 9th International Conference on Urban Drainage, pp. 1-14, 2002.

[19] W. Fan and R. B. Machemehl, "Using a simulated annealing algorithm to solve the transit route network design problem," Journal of Transportation Engineering, vol. 132, no. 2, pp. 122-132, 2006.

[20] H.-J. Shieh and R. C. Peralta, "Optimal in situ bioremediation design by hybrid genetic algorithm-simulated annealing," Journal of Water Resources Planning and Management, vol. 131, no. 1, pp. 67-78, 2005.

[21] F. Zhao and X. Zeng, "Simulated annealing-genetic algorithm for transit network optimization," Journal of Computing in Civil Engineering, vol. 20, no. 1, pp. 57-68, 2006.

[22] M. H. Chan, Application of real-valued-coding niche genetic algorithm to automatic calibration of an unsteady river flow model [M.S. thesis], National Taiwan University, Taipei, Taiwan, 2004 (Chinese).

[23] Y. C. Huang, Simulated annealing real-valued-coding niche genetic algorithm for automatic river simulation model calibration [M.S. thesis], National Taiwan University, Taipei, Taiwan, 2006 (Chinese).

[24] B. A. Tolson and C. A. Shoemaker, "Dynamically dimensioned search algorithm for computationally efficient watershed model calibration," Water Resources Research, vol. 43, no. 1, Article ID W01413, 2007.

[25] A. Behrangi, J. Vrugt, B. Khakbaz, and S. Sorooshian, "Comment on 'dynamically dimensioned search algorithm for computationally efficient watershed model calibration' by Bryan A. Tolson and Christine A. Shoemaker," Water Resources Research, vol. 44, no. 12, 2008.

[26] E. D. White, Z. M. Easton, D. R. Fuka et al., "Development and application of a physically based landscape water balance in the SWAT model," Hydrological Processes, vol. 25, no. 6, pp. 915-925, 2011.

[27] H. Yen, J. Jeong, W.-H. Tseng, M.-K. Kim, R. M. Records, and M. Arabi, "Computational procedure for evaluating sampling techniques on watershed model calibration," Journal of Hydrologic Engineering, vol. 20, no. 7, 2015.

[28] C. H. Tu, Application of dynamically dimensioned search algorithm to automatic calibration of an unsteady river flow model [M.S. thesis], National Taiwan University, Taipei, Taiwan, 2009 (Chinese).

[29] I. L. Wu, T. K. Tsay, and C. Lai, "An unsteady-flow simulation model for basin-wide river systems," in Proceedings of the 2nd International Conference on Urban Disaster Reduction, Taipei, Taiwan, November 2007.

[30] T. K. Tsay and I. L. Wu, "Development of basin-wide river numerical model and dynamic simulation for flood scenarios," in Proceedings of the 6th Taipei International Digital Earth Symposium, Taipei, Taiwan, May 2008.

[31] I. L. Wu, P. H. Chen, and T. K. Tsay, "Using a basin-wide river numerical model to evaluate flood control improvement measures-a case study on zhonggang main drainage at Xinzhuang District, New Taipei City," Journal of the Chinese Institute of Civil and Hydraulic Engineering, vol. 28, no. 1, pp. 45-55, 2016.

[32] C. Lai, Simulation of Unsteady Flows in a River System (Operational Manual), Hydrotech Research Institute, National Taiwan University, Taipei, Taiwan, 2002.

[33] C. Lai, T.-K. Tsay, C.-H. Chien, and I.-L. Wu, "Real-time flood forecasting," American Scientist, vol. 97, no. 2, pp. 119-125, 2009.

[34] I. L. Wu, T. K. Tsay, and C. Lai, "Numerical modelling of basin wide flood flow using variable flow-resistance coefficients," Journal of Flood Engineering, vol. 2, no. 1-2, pp. 99-120, 2011.

[35] C. Lai and T. K. Tsay, "Model development of keelung-river flood forecast system," Research Reports W10D90C008, Yen Ching-Ling Industrial Research Institution, National Taiwan University, Taipei, Taiwan, 2002 (Chinese).

[36] C. Lai, "Comprehensive method of characteristics models for flow simulation," Journal of Hydraulic Engineering, vol. 114, no. 9, pp. 1074-1097, 1988.

[37] R. W. Schaffranek and C. Lai, "Friction-term response to boundary-condition type in flow models," Journal of Hydraulic Engineering, vol. 122, no. 2, pp. 73-81, 1996.

[38] L. Brandimarte and M. K. Woldeyes, "Uncertainty in the estimation of backwater effects at bridge crossings," Hydrological Processes, vol. 27, no. 9, pp. 1292-1300, 2013.

[39] P. Costabile, F. Macchione, L. Natale, and G. Petaccia, "Comparison of scenarios with and without bridges and analysis of backwater effect in 1-D and 2-D river flood modeling," Computer Modeling in Engineering and Sciences, vol. 109, no. 2, pp. 81-103, 2015.

[40] T. J. Fewtrell, J. C. Neal, P D. Bates, and P J. Harrison, "Geometric and structural river channel complexity and the prediction of urban inundation," Hydrological Processes, vol. 25, no. 20, pp. 3173-3186, 2011.

[41] T. Picek, A. Havlik, D. Mattas, and K. Mares, "Hydraulic calculation of bridges at high water stages," Journal of Hydraulic Research, vol. 45, no. 3, pp. 400-406, 2007.

[42] V. Tayefi, S. N. Lane, R. J. Hardy, and D. Yu, "A comparison of one- and two-dimensional approaches to modelling flood inundation over complex upland floodplains," Hydrological Processes, vol. 21, no. 23, pp. 3190-3202, 2007

[43] P Costabile, F. Macchione, L. Natale, and G. Petaccia, "Flood mapping using LIDAR DEM. Limitations of the 1-D modeling highlighted by the 2-D approach," Natural Hazards, vol. 77, no. 1, pp. 181-204, 2015.

[44] F. Huthoff, P. C. Roos, D. C. M. Augustijn, and S. J. M. H. Hulscher, "Interacting divided channel method for compound channel flow," Journal of Hydraulic Engineering, vol. 134, no. 8, pp. 1158-1165, 2008.

[45] S. Proust, D. Bousmar, N. Riviere, A. Paquier, and Y. Zech, "Nonuniform flow in compound channel: a 1-D method for assessing water level and discharge distribution," Water Resources Research, vol. 45, no. 12, 2009.

[46] Z. Cao, J. Meng, G. Pender, and S. Wallis, "Flow resistance and momentum flux in compound open channels," Journal of Hydraulic Engineering, vol. 132, no. 12, pp. 1272-1282, 2006.

Fu-Ru Lin, (1) Nan-Jing Wu, (2) Chen-Hao Tu, (3) and Ting-Kuei Tsay (1)

(1) Department of Civil Engineering, National Taiwan University, Taipei City 10617, Taiwan

(2) Department of Civil and Water Resources Engineering, National Chiayi University, Chiayi City 60004, Taiwan

(3) MWH Americas Incorporated, Taiwan Branch, Taipei City 10549, Taiwan

Correspondence should be addressed to Nan-Jing Wu;

Received 7 July 2016; Revised 19 October 2016; Accepted 18 December 2016; Published 1 February 2017

Academic Editor: Paolo Lonetti

Caption: Figure 1: Map of the Tamsui River basin.

Caption: Figure 2: Channel schematization of the Tamsui River basin.

Caption: Figure 3: Comparison of the default [n.sub.d] with the calibrated [n.sub.d] in test case.

Caption: Figure 4: Comparison of the default [n.sub.u] with the calibrated [n.sub.u] in test case.

Caption: Figure 5: Plots comparing the simulation water level hydrographs (test case).

Caption: Figure 6: Comparison of the convergence process of DDS with SARvcNGA.

Caption: Figure 7: Plots comparing the observed with the calculated water level hydrographs (calibration case).

Caption: Figure 8: Plots comparing the observed with the calculated water level hydrographs during Trami Typhoon.

Caption: Figure 9: Plots comparing the observed with the calculated water level hydrographs during Usagi Typhoon.

Caption: Figure 10: Cross sections and model calibrated variable resistance coefficients at (a) Wudu, (b) Dajhih Bridge, and (c) Taipei Bridge Station.
Table 1: The resistance coefficients nd and nu calibrated by DDS
algorithm for each reach.

Distance from     Reach    Cross            [n.sub.d]   [n.sub.d]
river mouth (m)   number   section number

23,870-21,215     1        T36.A-T32        0.035       0.039
21,215-16,775     2        T32-T24.A        0.025       0.034
16,775-13,600     3        T24.A-T20        0.03        0.036
13,600-9,830      4        T20-T14          0.024       0.028
9,830-8,920       5        T14-T13          0.034       0.023
8,920-6,075       6        T13-T09          0.024       0.027
6,075-0           7        T09-T00          0.03        0.026
27,495-21,215     8        H10.A-H01        0.034       0.022
49,497-43,647     9        K94-K80          0.064       0.04
43,647-38,737     10       K80-K68          0.047       0.035
38,737-36,748     11       K68~K61          0.02        0.038
36,748-32,832     12       K61~K50          0.036       0.032
32,832-29,522     13       K50-K43          0.026       0.025
29,522-25,879     14       K43-K35          0.045       0.028
25,879-19,417     15       K35-K19.A        0.025       0.033
19,417-18,267     16       K19.A-K17        0.043       0.025
18,267-16,692     17       K17-K14          0.03        0.02
16,692-14,545     18       K14-K10          0.053       0.037
14,545-11,445     19       K10-K05          0.03        0.049
11,445-8920       20       K05-K01          0.043       0.035

Table 2: Froude number during peak of Typhoon Saola
(at 08/02/2012 0:00).

Junction   Velocity (m/s)   Area ([m.sup.2])   Width (m)

1               1.87            1797.56         960.18
2               1.37            3985.03         2902.00
3               2.43            2648.25         1091.97
4               2.27            2838.12         1248.46
5               1.60            3676.42         2296.04
6               1.31            2744.94         2095.21
7               1.51            2577.26         1708.95
8               0.96            5396.25         5641.67
9               1.06            1198.21         1125.71
10              1.88             156.20          83.27
11              2.12             197.79          93.40
12              1.87             213.89         114.53
13              1.51             280.07         186.09
14              1.23             364.52         296.02
15              0.97             486.54         499.22
16              0.95             500.48         525.77
17              1.02             625.52         613.74
18              0.78             860.27         1097.29
19              0.92             717.89         784.58
20              0.84             896.98         1067.20
21              1.16             794.94          68743

Junction    Discharge      Froude number
number     ([m.sup.3]/s)

1             3365.21          0.21
2             5472.24          0.16
3             6422.53          0.23
4             6451.91          0.26
5             5886.68          0.14
6             3596.14          0.12
7             3886.77          0.18
8             5161.52          0.10
9             1275.37          0.11
10            293.00           0.28
11            418.83           0.27
12            399.43           0.28
13            421.51           0.22
14            448.87           0.16
15            474.19           0.11
16            476.41           0.12
17            637.53           0.14
18            674.45           0.10
19            656.87           0.10
20            753.91           0.09
21            919.27           0.13
COPYRIGHT 2017 Hindawi Limited
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2017 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Title Annotation:Research Article
Author:Lin, Fu-Ru; Wu, Nan-Jing; Tu, Chen-Hao; Tsay, Ting-Kuei
Publication:Mathematical Problems in Engineering
Article Type:Report
Date:Jan 1, 2017
Previous Article:An Optimization Method of Passenger Assignment for Customized Bus.
Next Article:A New Approach of Waveform Interpretation Applied in Nondestructive Testing of Defects in Rock Bolts Based on Mode Identification.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters