Analysis of Flood Storage Area Operations in Huai River Using 1D and 2D River Simulation Models Coupled with Global Optimization Algorithms

Abstract

This article addresses the issue of flood management using four flood storage areas in the middle section of Huai River in China which protect the important downstream city of Bengbu. The same areas are also used by the local population as residential and agricultural zones. An optimization problem is therefore posed, with two objectives of simultaneously minimizing the downstream flood risk in Bengbu city and the storage areas’ economic damages. The methodology involved development of river flood models using HEC-RAS, with varying complexity, such as 1-dimensional (1D) model with storage areas represented as lumped conceptual reservoirs, and 2-dimensional (2D) models with detailed representation of the terrain, land-use and hydrodynamics in the storage areas. Experiments of coupling these models with global optimization algorithms (NSGA-II, PESA-II and SPEA-II) were performed (using the HEC-RAS Controller), in which the two objective functions were minimized, while using stage differences between the river and the storage areas as decision variables for controlling the opening/closing of the gates at the lateral structures that link the river with the storage areas. The comparative analysis of the results indicate that more refined optimal operational strategies that spread the damages across all storage areas can be obtained only with the detailed flood simulation models, regardless of the optimization algorithm used.

1. Introduction

Flood Storage Areas (FSAs) have been used in many countries for flood protection of downstream areas during incoming floods. These are dedicated upstream areas where flood waters are diverted, resulting in reduction of the flood peaks downstream. Design and operation of FSAs have traditionally been focused only on the main objective of reducing the flood peak downstream [1,2,3]. Recently, the conditions within the FSAs themselves have attracted the interest of researchers and practitioners. As floods are rare events, the FSAs are often used as agricultural and residential areas, with increased vulnerabilities and risks that need to be taken into account during their operation. Furthermore, they have been identified as zones of environmental protection in integrated river basin management, acting as one of nature-based solutions for flood risk management [4,5]. Therefore, analyzing the conditions within the FSAs during their operation needs to be taken into account together with the main objective of reducing downstream flood risks, aiming at optimal operations that minimize the risks both within the FSAs and downstream. These kinds of problems have already been studied using simulation models [6]. However, there are no reports in existing literature on using simulation coupled with optimization approaches for design or operations of FSAs, which is the main topic of this article.

The approach of coupling simulation models with external optimization algorithms (sometimes called model-based optimization) has been used for diverse water-related applications [7,8]. In this work we present how this approach can be applied to the particular problem of optimal operations of multiple FSAs used for flood protection of a downstream city. The problem has some similarities with optimal reservoir operations, for which numerous approaches have been proposed, including model-based optimization [9], but it also has some important differences (parallel consideration of damages in several FSAs, analysis of both filling and emptying the FSAs if they are controlled by corresponding gates, inclusion of such gate operations as decision variables in the optimization, etc.). Another similar problem to which model-based optimization has been previously applied is the design of flood storage reservoirs as integral part of urban drainage networks [10], but here again the conditions inside the flood storage reservoirs were not considered, only the design parameters for their filling in order to reduce the flood hazard.

The case study for the work presented in this article is the middle section of Huai River in China (also known as Huai He river in some literature), where four FSAs are considered for usage when incoming floods threaten the downstream city of Bengbu. As the same FSAs are used as agricultural and residential areas by the local population, the risk of economic damage during their usage is quite high. The methodology for addressing this problem with model-based optimization approach has been incremental, starting from initial, comparatively simpler setup, towards more complex problem formulation and analysis. It consisted of three distinct parts. First, the initial setup of the simulation-optimization (already briefly presented in [11]), was with a 1-D HEC-RAS simulation model containing FSAs representation as lumped conceptual reservoirs, coupled with the widely used multi-objective optimization algorithm Non-dominated Sorting Genetic Algorithm II (NSGA-II, [12]). In the second part of the study, the same simulation model was coupled with three different multi-objective optimization algorithms in order to test whether such different optimization algorithms would yield different results. In the third part, the analysis was extended with modification of the simulation model, primarily in the representation of the FSAs by using their detailed terrain and land use information, and, by modelling the 2-D flow during their flooding. Following Section 2 in which the study area is presented in more details, the methodological setup with the above mentioned three parts is described in Section 3. Section 4 presents the obtained results from the three different methodological approaches. The discussion of these results is presented in Section 5, followed by conclusions in Section 6, where we indicate some limitations of the study that point towards future research.

2. Study Area

China has several large rivers causing significant flood risks that need to be managed. Most well-known are Yellow River in the north and the Yangtze River in the south, which have also brought largest and most devastating floods throughout history [13,14]. Huai River, located in eastern central China, is also very flood prone, especially in its middle and downstream reaches [15]. The catchment area of the river is about 270,000 km², and the average annual precipitation is about 880 mm [16]. Historically, the river had significant drainage capacity with outflow into the East China Sea, until the diversion of part of Yellow River in 12th century AD, which brought about continuous sedimentation of its downstream area [17]. Further water infrastructure interventions in 20th century have resulted in the current conditions with significant reduction of slope in the downstream sections and complete loss of the river’s connection to the sea. The current outflow point of the river is Hongze Lake [18,19]. These conditions have been important contributing factor to increased food risks in the middle and downstream sections of the river. Additionally, the precipitation distribution in the Huai River basin is highly uneven in both space and time. Spatially, the precipitation is much higher in the south, with annual average of 1400–1600 mm (reducing from mountainous regions towards the low-lying areas, and from the coastal zone inlands), and temporally, most of the annual rainfall (50–80%) comes during the wet period of June-October [20,21]. The uneven distribution of precipitation causes both floods and droughts in the basin, although floods are more prominent and more frequent disasters. The combination of these hydrological hazards with the recent increase in socio-economic development, especially rapid urbanization and more developed and high-value agricultural zones, has led to higher impact and increased risk. Major recent flood disasters have occurred in 1991, 2003 and 2007 (for an overview of floods in earlier periods, see [22]). Huai River Basin Commission, as the responsible agency for flood risk management has developed and implemented many different structural and non-structural measures for reducing and mitigating flood risks [16,23]. Structural measures include numerous flood retention reservoirs in the upstream part of the basin, dikes and diversions, channel modifications for increasing drainage capacity in the middle and downstream part, and flood protection storage areas. Non-structural measures include administrative and legal arrangements that regulate the shared responsibilities among many involved stakeholders, flood risk assessment and mapping for developing flood risk management plans, flood forecasting and warning systems, and decision support systems for flood control operations.

This study is focused on a section in the middle part of the Huai River, covering about 110 km of the river, upstream of Bengbu city (Figure 1). The considered section contains several flood storage areas (FSAs) that have been allocated for flood protection of Bengbu, which is an important regional center with about 3 million inhabitants. When incoming flood threatens to flood the city, these FSAs are deliberately flooded using control gates, which results in reduction of the flood hydrograph peak downstream. In the past (since the 1950-ies) there have been eight such FSAs in use, but with recent changes in flood risk management plans, some of these areas have been abandoned, either due to implementation of measures that have made them redundant (e.g., increased channel capacity), or due to their usage as target areas for allocation of evacuated population in cases of emergencies. Currently, only four FSAs are considered to be active for diverting incoming flood waters (see again Figure 1), with the following storage capacities: Shouxi (1140 × 10⁶ m³), Dongfeng (259 × 10⁶ m³) Tangyu (380 × 10⁶ m³) and Jingshan (534 × 10⁶ m³). As most of the time the FSAs are not used for storing flood waters, they have been developed as agricultural and residential areas. This means that in case of flood diversion, there is significant risk of economic damage and even loss of life. Optimal operations of the FSAs (i.e., how they are to be filled and emptied) can lead to reduction of these risks, while still ensuring that the downstream city of Bengbu remains protected. Such operations of FSAs are however far from obvious, especially when there are several FSAs that are considered simultaneously. An optimization problem can therefore be formulated, with FSAs operations as decision variables, and with two objectives: one related to the protection of Bengbu city, and, another one for minimization of economic damages in the FSAs. Our article addresses this problem by developing a simulation model for the flood propagation through the considered section, together with the FSAs operations using control gates, and then combining the simulation model with external optimization algorithms in order to determine optimal FSA operation strategies.

3. Data and Methods

3.1. Initial Model-Based Optimization: 1D-SA with NSGA-II

As depicted in Figure 2, the methodology for setting up the model-based optimization consists of two main components. The first one is development of a river simulation model using the HEC-RAS modelling system that simulates the propagation of a flood hydrograph through the considered section. Using available flood hydrograph data from 2007, this model is developed and calibrated. The model also contains the FSAs represented as lumped conceptual reservoirs. If part of the flood water volume is to be diverted to the FSAs in a controlled manner, this is modelled by introducing appropriate lateral structures (weirs), with geometry and parameters corresponding to the actual real world structures, which are equipped with control gates. Usage of the so-called ‘inlet’ gates, allows water to flow from the river into a particular FSA. If given FSA needs to be emptied, this is achieved with the so called ‘outlet’ gates, allowing water to flow from the FSA back into the river. HEC-RAS has several different methods for controlling the gate operations, one of which is the stage (water level) difference between the river and an FSA for the inlet gates (or vice-versa for the outlet gates). For a particular simulation, the user can specify values for the stage differences at which the gates need to be opened or closed. Some additional data are also required, such as speed of gate opening/closing, maximum gate openings, weir formula to be used, etc. In the approach taken here, the stage differences for controlling the inlet gates have been selected to be used as decision variables in the optimization problem.

This brings the second component in the methodology, which is the optimization problem formulation, and its solution using an external optimization algorithm. Two objectives that need to be minimized have been formulated: first one related to downstream risk at Bengbu city and a second one related to the total damages in the FSAs. Both of these objectives depend on the water levels (stages) in the river and in the FSAs, which in turn depend on the gate operations, controlled by the stage differences as decision variables. The external optimization is carried with a genetic algorithm (NSGA-II, in this case) that runs the HEC-RAS simulation through a number of generations (iterations) with particular population size (number of candidate solutions). For each model run (function evaluation) different values of stage differences are used (as ‘proposed’ by the optimization algorithm), and the values of the two objective functions are calculated from the water levels provided as output of the model. From one generation to another the algorithm progressively moves towards better solutions (with smaller values of the objective functions). At the end of the optimization, a set of optimal strategies is obtained, representing different trade-offs between the two objective functions considered (the so-called Pareto set). Note that this approach is different to the one commonly taken in reservoir operation optimization studies, where releases in discrete time periods from the reservoirs are considered as decision variables. In fact, when using stage differences as control variables they are treated as parameters to be optimized, an approach that is sometimes called optimization of myopic policies [25].

The connection between the two parts, HEC-RAS simulation model and the NSGA-II optimization algorithm is realized using MATLAB code that critically relies on the HEC-RAS Controller API (Application Programming Interface), which has recently been made available [26]. This API allows external programs to access variables and parameters of HEC-RAS models, so that they can be used for different purposes (automatic calibration, optimization, sensitivity and uncertainty analyses, etc.). In this case, the API was used at each model run, first to set the values of stage differences ‘proposed’ by NSGA-II as inputs for the next HEC-RAS simulation, and, secondly, to extract the required values of water levels needed for calculating the values of the two objective functions. The following two sub-sections present the two components of the methodology in more detail.

3.1.1. HEC-RAS Simulation Model

Data for setting up the HEC-RAS simulation model of the considered section in the middle part of Huai River were provided by Huai River Basin Commission. The main purpose of the model was to simulate the propagation of a flood hydrograph from 2007, which lasted about two months (late June–late August) and had a peak value of nearly 8000 m³/s, occurring on 11 July, 2007 (Figure 3). This hydrograph was provided as upstream boundary condition at Lutaizi station. Rating curve was provided as boundary condition at the downstream end. The total length of the model was 110 km, and 75 cross sections, measured in the field, were available along this length. Data on additional lateral inflows were available at three locations as indicated in Figure 4 and they were included in the model. After several tests, computational time step of 15 min was selected for coupling purposes with the optimization algorithm, with storage of output data (results) every hour. Initially, the model was calibrated on discharges, using data available at Huanian station located in the middle of the considered section (see again Figure 1), leading to values of Manning’s n roughness coefficient in the range 0.02–0.025 along the main channel. A secondary calibration was then performed on water levels at the same Huanian station, by modifying (re-calculating) the rating curve at the downstream end, upon obtaining some additional improved cross section data at that location. The final calibration results were satisfactory, with RMSE of about 0.5 m for water levels and 25 m³/s on discharges. Given that in the overall analysis water levels are used in the optimization problem, further improvements would certainly be desirable, but these could not be attained with the available data. For demonstrating the methodological approach for dealing with this problem the calibration results were acceptable.

The model used gate-controlled lateral structures for linking the 1-D river representation with the four FSAs, which were represented as lumped conceptual reservoirs (i.e., using one storage-volume relationship). This model was therefore named ‘1D-SA’ model (the name is also used in later comparisons with other, more complex models). An important next step with the 1D-SA model was the sensitivity analysis for determining upper and lower limits of stage differences for the inlet gates that should be used in the optimization. This was initially done for each FSA, separately for the inlet and the outlet gates, and separately for stage differences for opening or closing the gates. This was followed by analysis of combinations of inlet and outlet gates, with different controlling stage differences for opening and closing, together with different combinations of FSAs.

The first conclusion from this analysis was to fix the stage differences for outlet gates at values of 0.1 m for opening 0 m for closing. Furthermore, it was concluded that the stage difference for closing an inlet gate needs to be smaller than the stage difference for opening it. This conclusion led to the selection of decision variables (stage differences) and corresponding bounds, as indicated in Figure 5. Note that Dongfeng FSA does not have an outlet gate, and flood water there stays for a longer period of time, until it infiltrates or evaporates.

3.1.2. Optimization Formulation and Solution with NSGA-II

As mentioned earlier, the optimization problem was formulated with two objective functions, first one related to the downstream risk of flooding in Bengbu, and second one related to weighted sum of total damages in the storage areas, both of which need to be minimized:

MIN F_1 = w_1 (y(X) − W_l) + w_2 (d(X))

(1)

subject to:

$$y\left(X\right)<S_i$$

$$MIN F_2={\displaystyle \sum }_{i=1}^n{w}_i{g}_i\left(X\right)$$

(2)

where: X—vector of decision variables (stage differences); F₁—downstream risk of flooding Bengbu comprised of two components: flood severity and flood duration; y(X) —river water level at Bengbu (dm); W_l—flood warning level at Bengbu (20.3 m), used for estimating flood severity; d(X) —duration of river water level above flood warning level at Bengbu (days), used for calculating flood duration; w₁ and w₂—weights for contribution of the two components to F₁ (assumed equal); S_l—safeguard water level at Bengbu (22.57 m), used as constraint that should not be violated; F₂—total weighted flood damage in all storage areas; n—number of storage areas; w_i—weight for damage calculation in each storage area, proportional to number of inhabitants; g_i(X) —total flood damage in a storage area comprising estimates of agricultural and residential damage.

The constraint that ensures non-exceedance of the safeguard level at Bengbu (S_l) is handled by introducing a penalty to the first objective function. This penalty is an additional function that introduces significant increase in the value of that objective function for solutions with water level at Bengbu that are larger than the safeguard level, so that the genetic algorithm rejects such solutions during its iterations.

Equation (1) indicates that the first objective function is formed as weighted sum of two downstream risk factors for Bengbu city: the severity of the flood (how many decimeters is the water level above the warning level) and the duration of the flood (how many days is the river water level above the warning level). The units of decimeters and days were used in order to have comparable orders of magnitudes of the two different risk factors, which are added together for calculating the total ‘composite’ risk. This formulation was chosen as both of these factors are important when considering the overall risk of flooding the city. Furthermore, the formulation allows certain level of flexibility in varying the contribution to the overall risk from the two components, which can be used in future experimentation.

The total damage in the FSAs is the 2nd objective function, given in Equation (2). For any FSA it consists of estimated agricultural and residential damage, which depend on the water level in the FSA. For damage calculations, one critical water level value is the so-called ‘Zhuangtai’ level, above which all housing should be located (normally about 3 m above the bottom of the FSA). Using this, the estimates of damages for each storage area are calculated as follows:

• If maximum inundation level > Zhuangtai level

residential:

$$P=p\left[s_1\left(I_{max}-H\right)+s_2{I}_{max}\right]$$

agricultural:

$$C=ab{d}_1{I}_{end}$$

(3)

• If maximum inundation level < Zhuangtai level

residential:

$$P=p{s}_2{I}_{max}$$

agricultural:

$$C=ab{d}_2{I}_{end}$$

(4)

where: P—residential damage (RMB), C—agricultural damage (RMB), and, in equation 1 − g(X) = P + C, p—population in the storage area, s₁—housing asset value (RMB/person), s₂—equipment asset value (RMB/person), I_max—maximum inundation level (m), H—height from bottom of storage area to Zhuangtai level (m), a—farming area (ha), b—crop value per unit area (RMB/ha), I_end—final inundation level, d₁—high damage rate, d₂—low damage rate. The available and assumed values of these parameters are given in

Table 1. Data for calculating flood damages in the storage areas (adapted from [11]).

FSA	p (person)	a (ha)	H (m)	s1 (RMB/person)	s2 (RMB/person)	b (RMB/ha)	d1	d2
Shouxi	86,000	11,333	3	400	100	4500	0.6	0.4
Dongfeng	19,800	3.733	3	400	100	4500	0.6	0.4
Tangyu	25,800	6000	3	500	200	4500	0.6	0.4
Jingshan	11,700	5733	3	500	200	4500	0.6	0.4

:

The calculation of the formulated objective functions was implemented in MATLAB, to be used by the NSGA-II genetic algorithm in the optimization. As mentioned earlier, using the HEC-RAS Controller API the two components of HEC-RAS and NSGA-II were connected. After several test runs, an optimization was performed with 5,000 function evaluations (HEC-RAS model runs), with population size of 50 and 100 generations, which lasted about 20 hours on a standard laptop.

3.2. Analysis with Different Optimization Algorithms: 1D-SA with NSGA-II, PESA-II, SPEA-II

This analysis aimed at investigating whether the solutions would be different with unmodified 1D-SA model, but with different global optimization algorithms. For this purpose, the same model setup was now connected with three different evolutionary algorithms. Furthermore, this analysis was carried out in the Multi Objective Evolutionary Algorithm (MOEA) framework developed in Java [27], which contains implementations of dozens of evolutionary algorithms and numerous supporting features for analyzing the progress of evolutionary optimization and comparing different algorithms. The framework has three important classes/objects named Executor, Analyzer and Instrumenter, each with a distinct role. Executor is used for setting up and running the optimization, Instrumenter assists in analyzing the behavior of the algorithm during its run, while Analyzer offers comparative analysis of different algorithms with data collected after the runs.

Using the MOEA framework three different evolutionary optimization algorithms were tested: NSGA-II (same as the one used in the initial setup with MATLAB), Strength Pareto Evolutionary Algorithm II (SPEA-II) [28,29], and Pareto Envelope-based Selection Algorithm II (PESA-II) [30,31]. All algorithms are improved versions of their earlier versions. Details of the algorithms are not of main focus here, but interested readers can check the indicated references, as well as some comparative studies [32,33]. Using the MOEA framework, for each of the three tested algorithms several metrics (parameters) describing the progress of the evolution were extracted, among which were generational distance and hyper-volume. These are in fact indicators about the convergence of the optimization algorithms towards a global optimum. Their analysis may indicate when the evolutionary algorithm reaches stable solutions, beyond which there will be no further improvements in the values of the objective functions. This in turn can give an indication about the number of required function evaluations, in this case HEC-RAS model runs. The analysis with all algorithms was carried out with 10,000 function evaluations, using 200 generations with population size of 50.

3.3. Analysis with More Complex Simulation Model: 1D-SA(Terrain-based) and 1D-2D, with NSGA-II

In this step the analysis was extended with modification of the HEC-RAS simulation model, while keeping the original NSGA-II algorithm implemented in MATLAB for optimization. Modifications of the initial HEC-RAS model were primarily by introducing more detailed representation of the FSAs, using terrain and land use information. Two model variants were developed. The first one modelled the filling of the FSAs in a similar manner to the initial 1D-SA model (from the lowest point of the FSA), just that their conceptual representation was based on actual geometry obtained from terrain data. The storage-elevation curve was still used, but updated with more detailed terrain data. Furthermore, in this variant, the damage calculations were based on damage curves developed from land use data. This model was named ‘1D-SA(Terrain-based)’. In the second variant the HEC-RAS capabilities for 1D-2D flow modelling were used, and the FSAs were represented as 2D flow areas. Damages were calculated grid by grid, depending on the actual flooding and gridded land use. It needs to be noted that because the 1D-2D connections with gate-controlled lateral structure were not fully implemented in HEC-RAS, the actual gate operations were still controlled with stage differences as in the 1D-SA model. In additional processing step, the so obtained flows through the gates were supplied as boundary conditions to the 2D flow areas representing the FSAs for modelling the 2D flow hydrodynamics. This model was named ‘1D-2D’ model. More details on the data processing steps and the setup of these two models are provided in the following sub-sections. Optimization experiments were performed with 2500 function evaluations (50 generations with population size of 50), which was earlier confirmed to be sufficient for obtaining stable optimal solutions. These were carried out for three different models: 1D-SA (initial setup, for purposes of comparison), 1D-SA(Terrain-based) and 1D-2D.

3.3.1. Terrain and Land Use Data Processing

Raw digital elevation data from SRTM (https://www2.jpl.nasa.gov/srtm/), the Shuttle Radar Topography Mission, were downloaded with resolution of 30 × 30 m (Figure 6a). The HEC-RAS tool named RAS Mapper was then used for combining river cross-section data with the terrain data.

To estimate the land use needed for calculating damages in the FSAs a land cover map was developed from a Landsat satellite image (Figure 6b) downloaded from USGS web site (https://earthexplorer.usgs.gov/). Starting with selected polygons with known land cover data, supervised classification was performed using the maximum likelihood method in ArcGIS, to obtain a final land cover map, with four classes as shown in Figure 7. Further processing was then performed in RAS mapper to match the land use information with the computational grid in the 2D flow areas representing the FSAs.

3.3.2. 1D-SA(Terrain-based) Model Setup

As mentioned earlier, the setup of 1D-SA (terrain-based), was in fact quite similar to 1D-SA. The only difference was in the use of improved storage-elevation curves, which were now developed from the actual terrain data. However, the differences in these curves were not very large, especially not in the lower ranges of elevation, which were of most interest, because the FSAs would be normally flooded to a maximum depth of about 5–6 m.

Larger differences compared to the original 1D-SA model were introduced in the damage calculations, as the damage estimates were obtained from damage curves that were developed from the land cover data. These damage curves were calculated by introducing flood depths in the FSAs at steps of 0.1 m and calculating the damages from so ‘flooded’ areas by considering the underlying land cover. The same equations for damage calculations were used (see Equations (3) and (4)), just that they were applied for each grid, with the actual land cover type (residential or agricultural). The obtained damage curves for the four FSAs are shown in Figure 8.

It should be noted that the damage calculations are not in fact carried out in HEC-RAS, but as part of the optimization with NSGA-II. Upon extraction of water levels in the FSAs from the HEC-RAS outputs, the damages are simply read from the damage curves and used in the optimization, as described earlier. This approach allows faster overall optimization, while still using detailed terrain and land use data.

3.3.3. 1D-2D Model Setup

This model used the full 2D representation of flooding in the FSAs. Such a model allowed capturing of actual flooding dynamics and calculating damages in different grids as they occur in time. To achieve this, extraction and storage of flood depths for different steps in time (every six hours) from the 2D HEC-RAS model was performed, followed by calculation of damages for each flooded grid in a given time step, and aggregating damage calculations across all time steps and all FSAs. This procedure was performed after each model run, to obtain the value of the second objective function that guides the iteration of the optimization algorithm. These intermediate processing requirements prolong the computational time needed for the optimization, in addition to the already longer computational time required to model the 2D flow in the FSAs. In order to reduce the time of optimization, coarser 2D computational grid of 100 × 100 m was adopted for optimization, which with computational step of 15 min still required nearly 10 min per model run. With 2500 function evaluations, this led to a very long optimization of about 15 days on a standard laptop.

4. Results

4.1. Results of the Initial Setup with Model 1D-SA

The result obtained from the initial setup in terms of Pareto set (objective functions space) is presented in Figure 9. Clearly, the obtained solutions are in two separate clusters, with the indicated solutions A and B having practically same downstream risk (1st objective), but very different values for the damages in the FSAs (2nd objective). Summary data for these two solutions are provided in

Table 2. Damages (weighted by population) for solutions A and B (adopted from [11]).

Solution	Total	Shouxi	Dongfeng	Tangyu	Jingshan
A	1.96 × 108	1.56 × 108	1.61 × 107	2.05 × 107	3.84 × 106
B	4.03 × 107	0	1.62 × 107	2.02 × 107	3.86 × 106

and

Table 3. Decision variables (stage differences in m) for solutions A and B (adopted from [18]).

Solution	x1	x2	x3	x4	x5	x6	x7	x8
A	11.1	4.3	8.6	3.1	6.6	0.4	4.8	0.5
B	/	/	8.8	3.1	6.6	0.5	4.8	0.5

. The results show that all solutions in the upper cluster include usage of the largest FSA—Shouxi, whereas the lower cluster of solutions does not include the usage of that FSA. However, intermediate solutions, in which Shouxi FSA is only partly used, have not been identified.

For any solution from the Pareto set, further detailed analysis of the actual gate operations can be performed using the HEC-RAS output. For demonstration purposes Figure 10 shows the actual operations of the gates in two FSAs (Dongfeng and Tangyu) with solution B.

The reduction of flood hydrograph at Bengbu station with solution B is shown in Figure 11. In terms of the risk factors included in the 1st objective function there is change from having severity of 5 dm above warning level and duration of 15 days in the original hydrograph, to severity of about 2 dm and duration of 4.5 days with solution B. Figure 11 also shows that with solution B the flood hydrograph is reduced around the peak, while later on the water level is somewhat higher than the original situation, as water is released back into the river from the FSAs.

The obtained result from the initial setup demonstrated the feasibility of the proposed approach for identifying optimal operation strategies. However, the separation of the optimal strategies in two clusters required further investigation. One possibility for obtaining more diverse set of solutions was sought in the second part of the study by employing different types of optimization algorithms.

4.2. Results of the Analysis with Different Optimization Algorithms

The first result presented is from the comparison of NSGA-II run in MATLAB and in MOEA (Figure 12a). The Pareto sets are quite similar, and the two clusters identified earlier are still present. Figure 12b presents the Pareto sets obtained from the three different algorithms. Although there are some small differences in the identified solutions, the same clustering is still present and the solutions are grouped depending on whether Shouxi FSA is used or not. Further analysis of individual solutions in terms of actual operations did not reveal significant differences in operations compared to those indicated in Figure 10 above.

Finally, the results for generational distance and hyper-volume variations are presented in Figure 13, which indicate that stable optimal solutions are in fact obtained after about 2,500 function evaluations, regardless of the algorithm used.

In terms of running time of the algorithms the results are quite similar, as indicated in

Table 4. Running times of the three tested optimization algorithms.

Algorithm	Frequency of Data Extraction (per N of Function Evaluations)	Running Time (min)
NSGA-II	10	409.3
	50	406.9
	Without data extraction	350.4
PESA-II	10	409.1
	50	406.3
	Without data extraction	349.9
SPEA-II	10	440.2
	50	435.9
	Without data extraction	420.5

which gives the values per 1000 function evaluations. The frequency of data extraction by the MOEA Instrumenter for calculating the required metrics influences the speed of execution, but in similar amount for all three algorithms. Without such data extraction the running times are shorter by about 10%.

The main conclusion from this part of the study is that different optimization algorithms produce similar results for the posed problem, with clustered solutions for optimal FSA operation strategies, depending on the usage of Shouxi FSA. Furthermore, this analysis indicated that about 2500 function evaluations are sufficient for obtaining stable optimal solutions. This information was used in the next part of the study where different optimization results were sought via the modification of the HEC-RAS simulation model.

4.3. Results of the Analysis with Different Simulation Models

First, the results from the obtained Pareto sets at the end of the optimization are presented in Figure 14.

Clearly, when using the two new models, 1D-2D (terrain-based) and 1D-2D, some new solutions are identified, located between the two clusters identified earlier. These solutions are in fact indicating operational strategies with more partial usage of Shouxi FSA. Second important conclusion from these results is that usage of the detailed terrain and land use information in the new models leads to much lower damage estimates in the FSAs, by an approximate factor of five. This is due to the more refined spatial flooding pattern produced by these two models, and due to exclusion of some areas from the damage calculations, such as natural forested areas or open water.

Three characteristic solutions are identified in Figure 14 for more detailed analysis, one for each of the models used in the optimization. These three solutions are selected because they provide very similar values for the downstream risk at Bengbu city, while having different damages in the FSAs.

Figure 15, Figure 16 and Figure 17 show the obtained flooding pattern and the gate operations, for each selected solution. Note that even though the initial 1D-SA model was run without any detailed terrain representation, the flooding in the FSAs is shown using detailed terrain, by providing the 1D-SA calculated flows into the FSAs. This is done for comparison purposes with the two new models.

Figure 15 indicates significant flooding in the FSAs, particularly in Shouxi, Dongfeng and Tangyu. This is confirmed with the gate operations for these FSAs (note again that Dongfeng FSA does not have an outlet gate). In particular, the inlet gate for Shouxi is opened once (around the time of the hydrograph peak) and stays continuously open for about 9 days. As Shouxi is the largest and most populated FSA, this brings about large damages. The gates of Dongfeng and Tangyu open slightly before the hydrograph peak, and they also stay continuously open quite long (Tangyu about 20 days and Dongfeng about 6 days). The inlet gate of Jingshan FSA on the other hand opens around the hydrograph peak, but stays open only for a couple of days. These gate operations modify noticeably the downstream stage hydrograph at Bengbu, particularly by reducing its peak, and by stage increase later-on after the flood peak, when outlet gates release water back into the river.

Figure 16 shows the results for one of the ‘compromise’ solutions obtained with the optimization using 1D-SA(Terrain-based) model. Compared to results from Figure 15, the flooding of Shouxi FSA is quite reduced and Dongfeng is practically not used (only two short inlet gate openings). The inlet gates of Shouxi are with four successive open/close operations, starting about three days before the hydrograph peak, and all together lasting about 6 days. The inlet gates of Tangyu and Jingshan FSAs are opened at similar moments in time, just before the hydrograph peak, with Tangyu gate opening slightly earlier. Both of these gates stay open for quite a long time (12–13 days), which brings significant flooding and damages in these two FSAs. The effects on the downstream stage hydrograph are similar as in the previous case, shown in Figure 15.

The results from the selected optimization solution with 1D-2D model presented in Figure 17 are quite similar to those of Figure 16 when 1D-SA (terrain-based) was used. Indeed, the main difference is that there are more individual open/close inlet gate operations for Shouxi FSA (eight), but all occurring over the same period as in Figure B. The damages in that FSA are therefore similar as in the previous case, even though the flooding pattern is somewhat different, because the 1D-2D model fills the FSA as water enters from the river via the lateral structure through the gate, while the 1D-SA (terrain-based) fills up the FSA from below, starting from its lowest elevation (see upper left parts of Figure 16 and Figure 17 with flooding patterns in Shouxi FSA). The gate operations of the other three FSAs presented in Figure 17 are very similar to those of Figure 16, with some small modifications for the inlet gate of Dongfeng FSA and somewhat increased flooding.

5. Discussion

The overall results presented in the previous section indicate that the use of more complex models in the optimization leads to identification of ‘compromise’ solutions, outside of the clusters identified earlier with the initial model 1D-SA. These solutions are with more partial use of Shouxi FSA. Also, because of the more detailed representation of the terrain and land cover, the overall damages are much smaller compared to the initial model where the damage calculations were lumped per whole FSA. In terms of flood damage calculations during optimization and identification of corresponding operational strategies—there is no significant difference between the 1D-SA (terrain-based) and 1D-2D model. This is an indication that the much faster 1D-SA (terrain-based) model can be used in optimization and still give reliable results. The flooding patterns, however, are modelled more realistically with the 1D-2D model.

The obtained results also indicate that the more complex HEC-RAS models during optimization provide more realistic damage estimations and this leads to identification of more refined operational strategies that distribute the damages across several different FSAs. However, the most complex 1D-2D model requires very long computational time and it is less suitable for coupling with external optimization. The ‘intermediate’ 1D-SA (terrain-based) model runs fast and can be used in optimization, while still giving very similar optimal operation strategies to the 1D-2D model. Therefore, a strategy can be adopted where the 1D-SA (terrain-based) model is used in optimization to determine several potentially good strategies of FSAs operations, which can then be tested and potentially refined with the 1D-2D model using simulation only. Such a strategy in a way confirms the well-known view of optimization as a methodology not necessarily used to find the best solution, but rather a methodology to eliminate a large number of poor solutions.

Considering the methodology of coupling of HEC-RAS simulation models with external optimization algorithms the important contribution of this work is the use of stage differences between the river and the FSAs from the simulation model as decision variables in the optimization. This usage of stage differences as control/decision variables is of particular significance, because it effectively reduces the number of decision variables used in the optimization. The complex system of the middle section of Huai River presented here was optimized using only eight decision variables. More traditional approaches used in reservoir optimization studies that use releases in particular periods of time as decision variables require discretization of time and would lead to larger number of decision variables and more complex optimization problem.

The methodology presented relies on developing complex HEC-RAS models, and their manipulation using the HEC-RAS Controller API. The work demonstrates the feasibility of performing model-based optimization using these tools, which is a confirmation of their value, especially of the API, which has recently been developed.

6. Conclusions

The presented methodology in this article can serve as basis for developing operational strategies to be used in actual management of the FSAs. Nevertheless real world application of this methodology can be achieved after addressing a number of additional important issues. First, the proposed methodology is primarily useful for designing operational strategies that could be pre-prepared and selected for implementation in case of emergency when incoming flood hydrographs would be forecasted. However, the incoming flood hydrograph may be with very different characteristics compared to the one used in this study. Therefore, the immediate next step is to address the uncertainty of the optimal strategies of FSAs operations with respect to different incoming flood hydrographs. Such analysis can lead to identifying different sets of operational strategies for flood hydrographs with different characteristics (multiple peaks, different timing of the peaks, different steepness of raising and falling limb, etc.). On basis of such analysis, further explorations for employing parts of this methodology for real-time FSAs operations can be pursued.

Further improvements and extensions of this work can be in other directions as well. The developed HEC-RAS models can be improved and refined with more discharge and water level data for better calibration. More accurate data for the terrain and land use in the FSAs could lead to better damage estimates, together with a more refined computational grid. This may lead to better and more detailed estimation of economic damages, as well as obtaining more accurate flooding pattern. However, such improvements are recommended to be used with the 1D-SA(Terrain-based) model during optimization, because of reduced computational time. The 1D-2D model can then be used only for simulation to test selected promising solutions resulting from optimization.

Regarding the formulation of the first optimization objective, the equal weights of the two risk factors related to severity and duration can be revised, or even a different formulation can be proposed. The hard constraint of avoiding flood levels above safeguard level is important, but it seems that the overall risk with the current formulation is somewhat overestimated and alternative formulations can be explored.

Although this study indicates that the choice of optimization algorithm does not have large influence on the obtained solutions, further analysis and possible improvements can certainly be sought regarding algorithm speed and efficiency in terms of identifying stable optimal solutions.