A novel simulation–optimization strategy for stochastic‐based designing of flood control dam: A case study of Jamishan dam

This study presents a novel stochastic simulation–optimization approach for optimum designing of flood control dam through incorporation of various sources of uncertainties. The optimization problem is formulated based on two objective functions, namely, annual cost of dam implementation and dam overtopping probability, as those are the two major concerns in designing flood control dams. The nondominated solutions are obtained through a multi‐objective particle swarm optimization (MOPSO) approach. Results indicate that stochastic sources have a significant impact on Pareto front solutions. The distance index (DI) reveals the rainfall depth (DI = 0.41) as the most significant factor affecting the Pareto front and the hydraulic parameters (DI = 0.02) as the least. The dam overtopping probability is found to have a higher sensitivity to the variability of stochastic sources compared to annual cost of dam implementation. The values of interquartile range (IQR) indicate that the dam overtopping probability is least uncertain when all stochastic sources are considered (IQR = 0.25%). The minimum annual cost of dam implementation (2.79 M$) is also achieved when all stochastic sources are considered in optimization process. The results indicate the potential of the proposed method to be used for better designing of flood control dam through incorporation of all sources of uncertainty.


| INTRODUCTION
Flood is a natural hazard which can cause significant environmental and social losses. Therefore, flood control is an important issue for the mitigation of flood damages across the world (Sharafati & Zahabiyoun, 2014). Construction of floodwater storage facilities such as reservoir dam is the most common strategy of flood control which can also provide other benefits, including water supply during the dry season and power generation (Ghosh, 1997).
Dam overtopping is a major concern of flood control dams. Destruction of about 35% of dams globally occurs due to overtopping (Sun, Chang, Miao, & Zhong, 2012). Therefore, estimation of required flood control volume and occurrence probability of overtopping are the most important components of designing and operation of flood control dams, especially for earth/rock-fill dams. Generally, the flood control volume and the probability of overtopping are estimated deterministically based on the design flood return period and flood routing (Paik, 2008). The main disadvantage of this deterministic methods is their inability to address the uncertainties arise due to variability in hydraulic and hydrological properties (Kwon & Moon, 2006).
Several analytical methods have been proposed to address the uncertainties in the estimation of flood control volume and the probability of dam overtopping (Afshar & Mariño, 1990;Gui, Zhang, & Wu, 1998;Kwon & Moon, 2006). Besides, some approximation approaches such as First Order Variance Estimation (Lee & Mays, 1986;Paik, 2008;Yenigun & Erkek, 2007), Rosenblueth Probabilistic Point Estimation (Kuo, Yen, Hsu, & Lin, 2007), Monte Carlo Simulation (Ardeshirtanha & Sharafati, 2020;Sharafati, Yaseen, & Pezeshki, 2020;Yuan et al., 2017) and Latin hypercube sampling (Hsu, Tung, & Kuo, 2011) are widely utilized to assess dam failure probability by considering the uncertainties associated to hydrological and hydraulic parameters. Though the consequences of dam failure are often devastating, their probability of occurrence is usually very low. Estimation of risk of such rare events using analytical methods needs joint probability function to get an accurate estimation. Therefore, estimation of dam failure probability using analytical methods is often not precise adequately. To improve risk estimation efficiency, different stochastic methods have been proposed in recent years. Though the estimation of dam overtopping probability has become an important topic in stochastic hydrology, use of stochastic simulation-optimization for the optimal stochastic-based design of dam projects is still very limited (Afshar, Masoumi, & Solis, 2015;Kuo, Hsu, Tung, Yeh, & Wu, 2008;Mousavi, Anzab, Asl-Rousta, & Kim, 2017).
In any construction project, an optimal stochasticbased design is decided based on a compromise between the risk of project failure during the lifetime and the initial cost of the project implementation (Afshar, Rasekh, & Afshar, 2009). Therefore, the optimal design of flood control dam project is a multi-objective problem where the objectives are to minimize the annual cost of dam implementation and dam overtopping probability. The classical optimization approaches such as linear, nonlinear, and dynamic programming are often not adequate to provide an optimal solution of multi-objective problems because these approaches have a single solution (Reddy & Nagesh Kumar, 2007). Meta-heuristic optimization approaches such as MOPSO is often suggested which can generate neuromas Pareto front solutions in each run. The MOPSO method provides an opportunity to select a favorable solution from several Pareto front solutions.
A novel multi-objective stochastic simulationoptimization approach is proposed in this study for optimal designing of flood control dam. A stochastic rainfall pattern generator, a rainfall-runoff model and a flood routing model are combined to develop the stochastic simulation model. The MOPSO optimization approach is linked with the simulation model to determine the optimum sizes of the flood control dam. The proposed method is utilized for designing flood control facility of Jamishan dam in Iran. It is expected that the stochastic optimization method developed in this study through the incorporation of all sources of uncertainty can be used for the estimation of dam overtopping probability with least uncertainty and implementation of dam with least cost.

| DAM OVERTOPPING PROBABILITY
Dam overtopping is a common failure mode which is defined as the loading exceeding the resistance. The resistance term is defined as the height of dam or elevation of dam crest (E Dam ) which is assumed to be constant. Besides, the loading term is adopted base on the maximum reservoir water level (E Reservoir Max ) which comprises the maximum routed flooding level (E Flood Max ) and the wave height (H Wave ) as follows: The wave height can be determined as follows: where H S and H r are the wind set-up and wave run-up, respectively. The maximum routed flooding level (E Flood Max ) can be determined using reservoir routing equation as follow: where f(.) is a function, I ave and Q ave are the average inflow and outflow discharges respectively, Δt is the time interval, and Δs is the reservoir volume changes. The inflow discharge can be simulated using rainfall-runoff modelling, and the outflow discharge can be computed using the spillway rating curve. Equation (3) indicates that the maximum routed flooding level (E Flood Max ) in a reservoir depends on inflow (I ave ), outflow (Q ave ) and initial reservoir water level. The volume of a long duration flood is very high while the peak inflow (I ave ) is usually less or equal to maximum outflow capacity, which means the change in the reservoir water level is less or equal to 0 ΔS Δt ≤ 0 À Á . Hence, the routed flooding level is less or equal to the initial reservoir water level. In contrast, a short-period flood with I ave > Q ave can cause a routed flooding level more than the initial reservoir water level ΔS Δt > 0 À Á . The overtopping can occur for a flood with large peak and volume. Besides, the highest peak discharge can occur when rainfall duration is close to the time of concentration. Hence, the deterministic value of rainfall duration is considered close to the time of concentration to estimate the high peak discharge, while the high flood volume is obtained through the simulations of high rainfall depth events under saturated soil condition.
The probability of dam overtopping (P f ) can be computed as follows (Tung, Yen, & Melching, 2006): where P (Á) denotes the probability. The Equations (1)-(4) reveal that uncertainty associated with dam overtopping probability is related to the variability of maximum reservoir water level. Hence, the possibility of estimation of dam overtopping probability is explored in this study by employing meteorological variables, hydrological parameters, hydraulic parameters and operational variables as the random variables.
A Monte Carlo (MC) simulation is employed to address the considered uncertainties and simulation of a large number of system outputs.

| Uncertainties in meteorological variables
The uncertainties in meteorological variables like rainfall and wind can have significant effect on the variability of dam overtopping. Therefore, the uncertainties in rainfall and wind are incorporated in the estimation of dam overtopping probability.

| Quantification of rainfall uncertainty
The Rain Pattern Generator (RPG) model, which is a stochastic model, is utilized to generate random rainfall events by considering uncertainty in rainfall features including rainfall depth (DE), rainfall duration (DU), and rainfall pattern (TY). The RPG extracts the empirical conditional probability distribution of each rainfall variable based on observed rainfall events and then generates rainfall events. For this purpose, all the observed rainfall events are preprocessed through a process consists of four phases, namely, filtering, classification, extracting the dimensionless matrix and derivation of probability distribution. In the filtering phase, all employed events (2,658 number of events) are filtered based on the minimum values of the rainfall depth (2 mm) and duration (100 min). The remaining events (1,724 of events) are used for classification, where the events are divided into three main groups: (a) the depth group which is further divided into six subgroups (DE 1 , DE 2 , …, DE 6 ), (b) the duration group, consists of five subgroups (DU 1 , DU 2 , …, D U5 ), and (c) the pattern group comprising four subgroups (TY 1 , TY 2 , …, TY 4 ). The events are made nondimensional through the third phase. For each subgroup, a nondimensional matrix is extracted using its nondimensional events.
The last phase is the derivation of probability distribution. In this phase, the observed events are divided based on the depth subgroups (e.g., DE 1 , DE 2 , …, DE 6 ). Then, the events stored in depth subgroups are classified based on their duration and pattern. For instance, in DE 1 subgroup, there are other combinations such as (DE 1 − DU 1 , …, DE 1 − DU 5 ) and (DE 1 − TY 1 , …, DE 1 − TY 4 ) where, DE 1 − DU 1 comprises the events of DE 1 and DU 1 subgroups. Ultimately, the relative frequency of each combination is computed.
The sequence used for RPG is given below: 1. The rainfall events are divided into the subgroups of depth, duration and pattern. 2. For each subgroup, the occurrence probability is computed. 3. Using the findings provided in the steps 1 and 2, the probability distribution of rainfall features [e.g., P (DU i j DE j ) or P(TY k j DE j )] are obtained. 4. The frequency analysis is used to generate a rainfall depth. 5. Using the rainfall depth generated in the step 4 and the P(DU i j DE j ), a rainfall duration is generated. 6. Using the rainfall depth generated in the step 4 and the P(TY k j DE j ), a rainfall pattern is generated.

| Quantification of uncertainty in wind velocity
Wind-generated waves have a noticeable effect on dam overtopping. Thus, wind velocity is an important factor of dam overtopping. The generated waves by wind comprise two components, the wind set-up and wave run-up. Wind set-up is the rise in the reservoir level due to wind stresses on the surface of the water which can be calculated as follow (Hsu et al., 2011): where V (mph) is the average wind speed over the reservoir, F (miles) is the fetch distance, and D (ft) is the mean water depth along the fetch line. Wave run-up is the vertical extent above the still water level to which an incident wave run-up toward the upstream face of the dam. It can be calculated as follow (Hsu et al., 2011): where a and b are the coefficients whose values depend on the embankment geometry (Cheng, Yen, & Tang, 1982). In this study, the values of a and b are considered 1.16 and 8.88, respectively. Equations (2), (5), and (6) indicate that the main source of wave height variability is wind velocity. The uncertainty of wave height can be quantified using the probability distribution function of wind velocity. The extreme value distributions are widely used in frequency analysis of wind velocity (Hsu et al., 2011). The Gumbel distribution is used in this study for the estimation of uncertainty of wave height. For this purpose, the random wave height is computed using Equations (2), (5), and (6) and the random wind velocity is generated using MC.

| Hydrological uncertainty
The Hydrologic Modelling System (HEC-HMS) is used in the present study to simulate the flood events. The HEC-HMS, which can simulate complex hydrologic processes for a wide range of basins, was introduced by the U.S. Army Corps of Engineers (USCE). The HEC-HMS model comprises several modules including infiltration, runoff transform, channel routing, base flow, snowmelt, evapotranspiration, soil moisture, and automatic calibration.
To compute rainfall excess, HEC-HMS employs land surface infiltration models such as Soil Conservation Service (SCS) curve number, initial and constant-rate and Green and Ampt model. As transform functions, HEC-HMS utilizes SCS, Clark's, Snyder's, Kinematic wave, and ModClark unit hydrographs to compute direct runoff from rainfall excess. The baseflow is calculated using exponential recession and linear reservoir methods (USCE, 2000).
Herein, Green and Ampt infiltration loss model, SCS unit hydrograph and exponential recession baseflow are considered for simulating flood events. Furthermore, Generalized Likelihood Uncertainty Estimation (GLUE) method is used to quantify the uncertainty of HEC-HMS parameters including the Green and Ampt loss parameters, for example, initial loss, moisture deficit, wet suction, hydraulic conductivity and impervious area, baseflow parameters, for example, flow at the start of a storm, threshold flow and discharge recession constant, and SCS unit hydrograph parameter namely, lag time. For this purpose, a large number of model parameter sets are generated based on the prior probability distribution of model parameters (Table 1). The corresponding model outputs are compared with the available observation using root mean square error (RMSE) (Wang, Frankenberger, & Kladivko, 2006).
The behavioral and nonbehavioral sets are obtained based on a behavioral threshold value. For the ith behavioral parameter set, θ i , the likelihood measure, l(θ i ) and re-scaled likelihood weights, l w (θ i ) can be computed as follows (Sharafati, Yasa, & Azamathulla, 2018): where RMSE min is the lowest calculated RMSE value among K behavioral solutions. Assuming the prior parameter probability distribution as the uniform distribution, Equation (7) can be used for extracting the posterior parameter probability distribution (Freni, Mannina, & Viviani, 2009). The extracted posterior distributions can be used to generate HEC-HMS parameters and corresponding flood event in each MC simulation.
In this study, the behavioral parameter sets are identified using 10,000 LHS (Latin Hypercube Sampling) sampled parameter sets using a behavioral threshold value of 1% (Her & Chaubey, 2015). The posterior parameter probability distributions of HEC-HMS parameters are presented in Figure 1. Figure 1 shows that the probabilistic characteristic of each posterior parameter distribution is not uniform, which means the uniform distribution is not adequate for quantification of uncertainty of HEC-HMS parameters.

| Hydraulic uncertainty
The outflow of the ogee spillway can be computed using a rating curve as follow (Mays, 1999): where l e is the effective length of spillway crest, C d is the discharge coefficient, C 0 is the ogee crest coefficient, K sp and K He H d are the ratio of upstream face slope of the spillway and crest head to design head, respectively, and H is the head on spillway crest. The effective length of spillway crest (l e ) would be calculated as follows (Mays, 1999): where l n is the net length of crest, N P is the number of piers, K P is the coefficient of pier contraction, and K a is the coefficient of abutment contraction. In this study, the variability of C 0 , K sp , K He H d , K P , and K a are considered as the main sources of hydraulic uncertainty. The random characteristic of these variables is quantified using uniform and rectangular distributions, as presented in Table 2.

| Variability of initial reservoir water level
Based on the dam operation policy, the reservoir water level can fluctuate between the minimum and maximum operating levels. It means there can have a wide range of values of the initial reservoir water level. The initial reservoir water level is a crucial variable in the evaluation of dam overtopping probability. Therefore, it is considered as a random variable and quantified by triangular probability distribution with a (upper limit), b (peak location), and c (lower limit) parameter values of 1,587, 1,584.5, and 1,564 m, respectively. These values are obtained from the time series of reservoir water level.

| Development of overtopping probability assessment approach
A new approach is proposed in the present study for the assessment of dam overtopping probability. The procedure followed for the estimation of dam overtopping probability is outlined below.
1. Define the value of dam crest elevation (E Dam ) as resistance term. 2. Generate random rainfall events using RPG model. 3. Generate random parameters of the HEC-HMS model.
4. Simulate the HEC-HMS model using the results obtained in steps 2 and 3, and then compute corresponding hydrographs. 5. Define the length of spillway (l n ), number of piers (N P ) and the width of piers (W P ) to generate a random rating curve of ogee spillway. 6. Define the elevation of spillway crest value (E Sp ). 7. Generate random initial reservoir water level (IRWL). 8. Route the inflow hydrograph based on the results obtained in steps 4-7, and then compute the maximum routed flooding level (E Flood Max ). 9. Generate random wind speed to calculate the corresponding random wave height (H Wave ).
where n f is the number of failure simulations and n T is the total number of MC simulations.

| THE COST FUNCTION OF DAM
The annual cost of dam implementation (ACDI) which comprises of the initial construction (C c ) and annual maintenance costs (C m ) of dam is used to quantify the cost function of dam as follows (Yazdi, Golian, & Roohi, 2017): where α is the maintenance cost ratio which is assumed 0.05, CRF is the capital recovery factor which can be computed as follow (Peterson, 2005): where r and t are the interest rate and the lifetime period, which are assumed 0.1 and 50 years, respectively, in the present study.
The C c has a direct relation with the dimension of dam body and spillway size. In this study, C c is computed from dam crest elevation (E Dam ), spillway crest elevation (E Sp ), length of spillway (l n ), number of piers (N P ) and the width of piers (W P ). Based on the quantitative survey and estimation analysis performed in Jamishan dam, the following equation is obtained to compute the C c (million dollars): Equations (11)-(13) indicate that ACDI is a deterministic function which is only related to the dam body and spillway dimensions. Furthermore, it is assumed that the Jamishan dam would be destroyed entirely after overtopping. Hence the cost associated with the overtopping event is same as ACDI. The remediation cost of dam structure or the cost incurred due to downstream consequences are ignored in this study due to unavailability of the required data.

| Problem formulation
Two objective functions, namely the annual cost of dam implementation (ACDI) and the dam overtopping probability (P f ) are used to optimize the size of flood control facilities in dam project. Hence, the optimization problem can be formulated as follows: which are subjected to the following constraints: l Min n ≤ l n ≤ l Max

| Optimization approach
The MOPSO is used as an optimization approach in this study to solve the optimal stochastic-based design problem. Coello, Pulido, and Lechuga (2004) proposed MOPSO, which combines the Pareto front and the technique of grid making as a multi-objective optimization approach. In MOPSO, a new population set is initiated where each of the particles in population is evaluated. The positions of nondominated solutions are then stored in the repository as the first memory of particles. Next, the velocity of ith particle (V t i ) is updated as follow: where t is the number of iteration, w t is the weight of inertia in tth iteration which can be damped in the next iteration by multiplying damping ratio (r d ) as w t + 1 = w t Á r d , c 1 , and c 2 are the cognitive and social coefficients, respectively, r 1 and r 2 are the random values between 0 and 1, P t i,Best is the best position of ith particle in tth iteration, P t g,Best is the global best position in tth iteration which is taken from the repository, and X t i is the position of tth particle in tth iteration. The next position of the tth particle can be determined as follows: The new position of each particle is assessed, and the nondominated ones are added to the repository. Next, the extra member of the repository is randomly removed based on repository size. The process is repeated until it reaches a termination criterion such as maximum iterations. The parameters of MOPSO used in the present study are given in Table 3.

| Simulation-optimization framework
A new simulation-optimization approach is proposed in this study to determine the optimal stochastic-based design of flood control dam. The proposed simulation-optimization approach is implemented using the following steps: 1. Go to optimization section 2. Generate initial population of decision variables (e.g., dam crest elevation, length of the spillway, number of piers, the width of piers, and elevation of spillway crest) 3. Compute the spillway rating curve using generated spillway specifications in step 2. 4. Determine dam crest level using generated dam specification in step 2. 5. Go to the simulation section 6. Compute the dam overtopping probability using defined uncertainty sources and the parameters estimated in steps 3 and 4. 7. Get back to optimization section 8. Calculate the annual cost of dam implementation using specifications of dam and spillway generated in step 2. 9. Consider the calculated value of the annual cost of dam implementation (step 8) as the first objective function (F 1 ). 10. Consider the calculated value of dam overtopping probability (step 6) as the second objective function (F 2 ). 11. Repeat steps 3-10 for all members of the initial population. 12. Store nondominated solutions in the repository. 13. Based on the termination criterion (e.g., maximum iterations), apply leader selection and mutation operators on the initial member of the repository. 14. Store the final member of the repository as Pareto fronts points.
After performing simulation-optimization, the optimal values of dam crest elevation, length of the spillway, number of piers, the width of piers and elevation of spillway crest can be extracted from Pareto fronts points. The framework proposed for simulation-optimization is presented in Figure 2. After completion of the simulation, the value of dam overtopping probability is considered as the second objective function (F 2 ) in the optimization section ( Figure 2). The optimized specifications of dam and spillway are generated in the optimization section from the inputs considered in the simulation section.

| Distance index
To assess the effect of stochastic sources on obtained Pareto fronts, the defined cases are compared with the Base Case (defined in Table 4). The distance between Pareto front functions of Base Case and Case i (ith Case) is measured using distance index (DI) as follow (Iwayama, Hirata, & Aihara, 2017): where x min and x max are the minimum and maximum values of the feasible interval between the Pareto fronts of the Base Case and Case i , f Base Case and f Case i are the Pareto front function of Base Case and Case i , respectively.

| Total weight factor index
The selection of the best Pareto front solution is a major concern in MOPSO. In this study, the best Pareto front solution for each defined case is selected using a pseudo weight factor, W j i . This factor can be computed using Deb (Deb & Algorithms, 2001) formula as follow: where F Max i and F Min i are the maximum and minimum values of i th objective function, F j i is the jth value of the ith objective function and l is the number of repository members.
In a multi-objective optimization problem, each solution can have two or more pseudo weight factors based on the number of objective functions. Therefore, the total weight factor index, W j T is defined to select the best Pareto solution as follows: where M is the number of objective functions. The best Pareto solution has the highest value of the total weight factor.

| DESCRIPTION OF THE CASE STUDY AREA
The proposed optimization approach is applied to Jamishan dam, located in the west of Iran. Iran is dominated by the arid and semi-arid climate where the water is highly scarce. Besides, water existence is subjected to various constraints of watershed properties. The water shortages and limitations emphasis the need for engineering solution of the problems related to water management and sustainability. The catchment of Jamishan dam is an upstream sub-basin of the Seymareh river. The area of the catchment is 524 km 2 and the maximum elevation is 3,000 m. The average annual air temperature and precipitation of Jamishan catchment are 10 C and 441 mm, respectively. According to the basic design information, the height of Jamishan dam, total volume of the reservoir, length of spillway, dam and spillway crest levels are 53 m, 62.8 mcm, 30 m, 1,590 m, and 1,587 m, respectively. The observed meteorological data (e.g., precipitation and wind speed) are extracted from Songhor station, which is located near to the center of the catchment. A total of 2,658 observed rainfall events measured in 10-min time interval are used in this study for the extraction of the probability distributions of rainfall characteristic (e.g., depth, duration and pattern). 7 | RESULTS AND DISCUSSION

| Description of defined cases
In this study, the effect of stochastic sources on optimization results (e.g., Pareto front and best solution) is investigated using 13 defined cases (e.g., Base Case, Cases 1 -12). In Base Case, all employed variables (e.g., meteorological, hydrological, hydraulic, and operation variables) are considered as stochastic while in other cases (e.g., Cases 1 -12), at least one stochastic source is deemed to be deterministic. For instance, in Case 1, hydrological parameters are considered deterministic while other variables are considered stochastic. Details of the defined cases are presented in Table 4. Table 4 shows the effect of rainfall duration and depth on dam overtopping investigated using Cases 5 -8 and Cases 9 -12, respectively. It should be noted that the correlation between rainfall depth and duration is considered using Cases 5 -12. For instance, the rainfall depths are considered deterministic in Cases 9 -12 (rainfall depth for return period = 500, 1,000, 5,000, and 10,000 years) while rainfall durations are generated randomly by empirical conditional probability distributions (see Section 2.1.1).
The number of MC simulations is a crucial issue in the simulation of risk assessment problems. The appropriate number of MC simulations (N MC ) is computed using the following equation: where Zα = 2 is the is the (1 − α = 2 ) percentile of the standard normal distribution, σ is the standard division of the simulations, and ε is a defined acceptable difference between each simulation and the mean value of previous simulations (Ata, 2007). The overtopping probability, σ generated in Base Case is 0.21%. Hence, the adequate number of MC simulation is 6,776 considering α = 5% and ε = 0.005%. This study considers 10,000 simulations for assessing the convergence of MC results. For instance, values of dam overtopping probability against the number of simulations for a member of Base Case Pareto Fronts is shown in Figure 4. Figure 4 indicates the value of dam overtopping probability is converged to a fixed value after 7,000 simulations. Hence, it is considered that the overtopping probability assessment with 10,000 MC simulations can provide reliable results for the case study.

| Effect of stochastic sources on obtained Pareto fronts
The effect of the uncertainties in hydrological parameters, initial reservoir water level, hydraulic parameters, rainfall duration and depth, and wave height on optimization results is investigated through Cases 1 -4, 7, and 10, respectively. In each case, the corresponding stochastic source is considered deterministic. The deterministic values of parameters/variables are determined based on the following concepts: F I G U R E 3 Map of Jamishan dam catchment and the locations of the observation stations 1. The calibrated parameters as presented in Table 1 are considered as the deterministic values of hydrological parameters (Case 1) 2. The designed hydraulic parameters of Jamishan dam (Table 2) are considered as the deterministic values of hydraulic parameters (Case 2).
3. The initial water level of designed flood routing of Jamishan dam (conservative initial water level) is considered as the deterministic value of initial water level (Case 3) which is 1587 (m. a. s. l.) 4. The 100-year return period wave height (1.01 m) of Jamishan dam is considered as the deterministic value of wave height (Case 4) 5. The designed rainfall duration (12 hr) which is close to the time of concentration of Jamishan dam catchment is considered as the deterministic value of rainfall duration (Case 7) 6. The 1,000-year return period rainfall depth, which is also the design rainfall depth of Jamishan dam is considered as the deterministic value of rainfall depth (Case 10).
The deterministic values of rainfall depth (T = 1,000 years) and duration (close to the time of concentration) are considered based on the conventional design criteria. The obtained Pareto fronts of all the foregoing cases are compared with Base Case, as shown in Figure 5a.  Figure 5a reveals that among all the cases, the obtained Pareto front of Case 10 (rainfall depth case) had the highest distance from the Base Case (black curve). Therefore, it can be remarked that rainfall depth has the highest effect on optimization results among all the stochastic variables considered. The DI values obtained for all the cases (the distance between Base Case and other cases) are presented in Figure 5b. The polynomial equations are fitted to obtain the Pareto fronts of the defined cases and to extract the mathematical function of each Pareto front, which are given in Table 5. Figure 5b shows that rainfall depth (DI = 0.41), hydrological parameters (DI = 0.36) and initial reservoir water level (DI = 0.15) are the first, second, and third effective stochastic sources while the hydraulic parameters (DI = 0.02) has the least effect.
To assess the effect of stochastic sources on the uncertainty of Pareto fronts solutions, the boxplots of the foregoing cases are presented in Figure 6. Figure 6a shows that the minimum and maximum values of 50% percentile (Q 50% ) of overtopping probability belong to the Cases 10 (Q 50% = 0.65%) and 1 (Q 50% = 1.87%). This means the relative difference of Q 50% (RD Q 50% ) of overtopping probability is about 188% ( 1:87− 0:65 ð Þ 0:65 × 100). Figure 6b shows that the minimum and maximum values of Q 50% of the annual cost of dam implementation are 2.96 M$ (Case 2) and 3.03 M$ (Case 4), respectively, and the RD Q 50% value is 2.4% ( 3:03− 2:96 ð Þ 2:96 × 100 ). The results indicate that the stochastic sources have higher impacts on dam overtopping probability compared to dam construction cost. The comparison of the values of the interquartile range (IQR) of dam overtopping probability (Figure 6a) and the annual cost of dam implementation (Figure 6b) reveals that the minimum and maximum values of IQR belong to the Base Case (IQR = 0.25%) and Case 4 (IQR = 0.34%), respectively. Therefore, it can be remarked that Pareto front of Base Case has the least uncertainty of overtopping probability compared to other cases.

| Effect of rainfall uncertainty on obtained Pareto fronts
The effect of rainfall duration and depth on obtained Pareto fronts are investigated using Cases 5 -8 (rainfall duration = 6, 9, 12, and 15 hr) and Cases 9 -12 (rainfall depth for return period = 500, 1,000, 5,000 and 10,000 years). As demonstrated in Table 4, the rainfall durations are considered deterministic, and other uncertainty sources are generated randomly for Cases 5 -8. The similar condition exists for Cases 9 -12, where the rainfall depths are considered fixed, and other variables (e.g., rainfall duration, hydrological parameters, hydraulic parameters, initial water level, wave height) are generated based on their PDFs. For instance, all stochastic variables except the rainfall depth are generated randomly for Case 12 while all the stochastic variables generated randomly used for Base Case. Hence, the effect of rainfall depth uncertainty on obtained Pareto fronts can be captured by making a comparison between the results obtained in Base Case and Cases 9 -12.
The Pareto front curves and the boxplots of Pareto front solutions for the cases are presented in Figures 7 and 8. Figure 7b shows that the mean value of the interquartile range (IQR mean ) (0.23 + 0.39 + 0.33 + 0.32/4) for rainfall duration cases (Cases 5 -8) is 0.32 while it is 0.26 for rainfall depth case (Figure 8b) rainfall duration has a higher effect on the uncertainty of Pareto fronts solutions compared to rainfall depth. Figures 7b,c and 8b,c show that RD Q 50% values of overtopping probability in both rainfall duration and depth cases are 675% and 417%, respectively, while those for the annual cost of dam implementation are 6.8% and 2%, respectively. It means that RD Q 50% values of the annual cost of dam implementation are significantly smaller than overtopping probability, or the annual cost of dam implementation for all Paetro fronts are less sensitive to rainfall duration and depth compared to dam overtopping probability.

| Relation between stochastic sources and best solution
To obtain the best solution for the defined cases, the values of total weight factor index are computed using Equation (24) and presented in Table 6. The values of the extracted optimized objective functions after considering the effect of stochastic sources on optimized cost is presented in Figure 9. Figure 9 shows that the minimum annual cost of dam implementation belongs to Base case (

| Impact of rainfall uncertainty on the best solution
The results obtained in previous sections revealed that uncertainties in rainfall duration and depth have the most significant effect on the obtained Prato front. Therefore, the impact of the variations of rainfall duration and depth on the best solution is assessed using Cases 5 -12. The annual cost of dam implementation and the overtopping probability of the best solutions for the mentioned cases are presented in Figure 10. Figure 10a,c reveals that the optimal cost of flood control facilities decreases with the increase of rainfall duration and increases with the increase of rainfall depth. At the same time, dam overtopping probability has vice versa relationships (Figure 10b,d).

| SUMMARY AND CONCLUSIONS
A novel multi-objective stochastic simulationoptimization approach by combining RPG, HEC-HMS, and flood routing models is presented in this study for the optimal stochastic-based design of flood control dam. The MOPSO optimization approach is linked with the runoff simulation model to determine the optimum sizes of the dam flood control system such as dam crest elevation, flood control height, length of the spillway, number and width of spillway piers. The annual cost of dam implementation and dam overtopping probability are utilized as two objective functions. The combination of various sources of uncertainty such as meteorological, hydrological, hydraulic, and dam operation is considered using 13 unique cases.
The results reveal that the obtained Pareto front strongly depends on the stochastic variables used. The impact of the sources of uncertainty on obtained Pareto front, quantified using distance index (DI) revealed that the uncertainty of rainfall depth (DI = 0.41) has the highest impact on the obtained Pareto front followed by hydrological parameters (DI = 0.36) and initial reservoir water level (DI = 0.15), while hydraulic parameter (DI = 0.02) has the least effect. The RD Q 50% values reveal that dam overtopping probability is highly sensitive to stochastic sources compared to the annual cost of dam implementation. The solutions of Pareto fronts show the least uncertainty in dam overtopping probability (IQR = 0.25%) when all stochastic sources are considered (Base Case). It is also found that the rainfall duration (IQR mean = 0.32) has a higher effect on the uncertainty of Pareto fronts compared to the rainfall depth (IQR mean = 0.26). The total weight factor index used to select the best Pareto front solution indicates an increase in the dam overtopping probability and optimum cost of flood control facilities when stochastics sources are ignored. Overall, it can be concluded that stochastic variables should be considered in a simulation-optimization model to remove the effect of random variables on the optimum design of a flood control dam.
The spillway discharge simulation method proposed by Mays (1999) (Equations (8) and (9)) is adopted in this study. In future, other complex approaches based on computational fluid dynamics (CFD) can be used for better assessment of overtopping probability. Besides, other costs, including dam remediation cost or the cost incurred due to downstream consequences, can be considered into optimization analysis in future.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon justifiable request.