Accounting for the uncertain effects of hydraulic interactions in optimising embankments heights: Proof of principle for the IJssel River

Most alluvial plains in the world are protected by flood defences, for example, embankments, whose primary aim is to reduce the probability of flooding of the protected areas. At the same time, however, the presence of embankments at one area influences hydraulic conditions of downstream areas located on the same river. These hydraulic interactions are often neglected in current flood risk management. The aim of this study is to explicitly acknowledge hydraulic interactions and investigate their impact on establishing optimal embankment heights along a stretch of the IJssel River. We find that the current approach leads to a single solution, while taking into account hydraulic interactions substantially expands the number of promising solutions. Furthermore, under a reference scenario, the current approach is in fact suboptimal with respect to both downstream locations and the system as a whole. Under uncertainty, it performs adequately from a system viewpoint, but poorly for individual locations, mostly due to risk overestimation downstream. Overall, the current approach proves to be too short‐sighted, because spatial trade‐offs among locations are neglected and alternative solutions remain hidden. Acknowledging the effect of hydraulic interactions provides policy makers with a broader and more comprehensive spectrum of flood risk management strategies.


| INTRODUCTION
The alluvial plains along most large lowland rivers around the world (e.g., the Rhine, the Po, and the Elbe River) are protected against flooding by embankments or other flood defences. Embankments have the primary aim of reducing the probability of flooding of the protected area and they have historically been the most commonly adopted flood risk reduction measure. Embankments sometimes substantially influence the way the protected alluvial plains develop, both economically and demographically (White, 1945). The Netherlands represents an emblematic case: structural defences, such as embankments, floodwalls, and dams have been built over the years and they currently amount to a total length of about 3,500 km, only accounting for the so-called primary defences (Kind, 2014).
Although embankments represent a successful flood risk management measure, their adoption is recognised to alter the hydrological regime of rivers. For example, Di Baldassarre, Castellarin, and Brath (2010) demonstrate the increase in the flood peaks experienced at a downstream location of the Po River as a consequence of the progressive enhancement of flood defences over time. Conversely, Van Mierlo et al. (2007) illustrate how potential breaches upstream lead to a reduction of flood load downstream. These examples illustrate the existence of complex, and yet understudied, hydraulic interconnections between planned interventions (e.g., raising embankment height) upstream and the associated unintended consequences (e.g., higher water levels or increased flood damage) downstream. This highlights the importance of considering what in the present paper is referred to as 'hydraulic system behaviour', i.e., the change in hydraulic loads at one location as a consequence of the state of the embankment system at other locations (Van Mierlo et al., 2007;Vorogushyn, Lindenschmidt, Kreibich, Apel, & Merz, 2012).
Several studies investigated the effect of hydraulic system behaviour on flood hazard and risk. For instance, Apel, Merz, and Thieken (2009) built a dynamic-probabilistic model to assess the effect of hydraulic system behaviour on flood frequency in contrast to traditional flood frequency analysis. They found that their model was able to provide a much more realistic flood frequency curve for downstream locations, especially for high discharge events, where embankment breaching mechanisms become relevant. Courage, Vrouwenvelder, van Mierlo, and Schweckendiek (2013) performed a flood risk analysis by comparing results with and without considering hydraulic system behaviour. They found that hydraulic system behaviour leads to different individual as well as economic risks. De Bruijn, Diermanse, and Beckers (2014) came to similar conclusions. They investigated the effect of hydraulic system behaviour on societal flood risk in the Netherlands and found that the number of fatalities N occurring with a given frequency F (N) is significantly lower when hydraulic system behaviour is considered. This applies especially to extreme flood events, which are those of more concern from a societal viewpoint. These studies show the relevance of taking into account the effects of hydraulic system behaviour in flood risk analysis and suggest that its inclusion may lead to different decisions and alternative investment schemes. However, a due consideration of hydraulic system behaviour in the design and planning of flood risk management measures is still lacking: current plans are usually based on flood risk analyses which assume hydraulic loads at each embankment as being independent from those nearby or upstream (De Bruijn, Diermanse, Van Der Doef, & Klijn, 2016;Vorogushyn et al., 2017). For instance, flood protection standards in the Netherlands are based on a well-known and successfully applied embankment height optimization model which was first introduced by van Dantzig (1956) and then developed further by, for example, Brekelmans and den Hertog (2012) and Eijgenraam, Brekelmans, Den Hertog, and Roos (2017). Although the value of these models is undeniable, flood probabilities of the protected areas are considered to be independent of one another. They thus ignore the change in the hydraulic load along the river stretch as a consequence of the state (e.g., failure, increase in safety) of embankments elsewhere.
The neglect of hydraulic system behaviour in flood risk management is because of two main reasons (De Bruijn et al., 2016). First, considering hydraulic system behaviour requires dealing with multiple uncertainties, such as the location of breaching ('which embankments will fail and in what order'?), the moment of breaching, the breach growth rate and the final breach width (Vorogushyn, Merz, Lindenschmidt, & Apel, 2010). Second, considering hydraulic system behaviour would further complicate the decisionmaking process in flood risk management (De Bruijn et al., 2016;Van Mierlo et al., 2007). Deciding upon flood risk management measures at one location would require taking into account the interests of communities elsewhere, both upstream and downstream of that location. In fact, the European Flood Directive (Directive 2007/60/EC, 2007 prescribes that this should be performed, but it is seldom applied in practice. The aim of this study is to investigate what taking into account hydraulic system behaviour might mean for the choice of optimal embankment heights along a stretch of the IJssel River in the Netherlands. The analysis is carried out applying the Many Objective Robust Decision-Making (MORDM) framework (Kasprzyk, Nataraj, Reed, & Lempert, 2013).
The paper is structured as follows: we (a) introduce the case study and the optimization problem, explain (b) the Many Objective Robust Decision Making framework, (c) the simulation model, (d) discuss the results, and (e) provide conclusions and recommendations for further research.

| THE CASE STUDY AND THE OPTIMIZATION PROBLEM
The IJssel River is a branch of the Rhine River in the Netherlands flowing north for about 125 km before discharging into the IJsselmeer. The present work focuses on a stretch of this river between the cities of Doesburg and Deventer (see Figure 1). Five locations of interest are identified, each representative of a different embankment stretch.
Typically, the problem of finding optimal embankment heights for each stretch is approached by searching for the least total costs solution (Kind, 2014). In other words, the optimal embankment height is considered the one for which the sum of embankment raising costs and expected annual damage is the lowest.
Embankment raising costs are simulated as in Eijgenraam et al. (2017): where u is the degree of embankment heightening; parameters c and b are fixed and variable costs, respectively, λ is a scale parameter and W is the cumulative embankment heightening over the entire planning period. Parameters c, b, and λ are assigned per stretch of the embankment system. However, some of the stretches considered in the present study are only a portion of those for which cost functions are defined. New cost functions are estimated by multiplying the original cost functions by the ratio of the lengths of the original embankment stretch and the one considered in the present study. Values of parameters c, b, and λ and the lengths of each stretch are shown in Table 1. The cost function is meant to apply to a sequence of optimal embankment height decisions over time. Therefore, W represents the increase in embankment heightening costs as a consequence of the previous heightening. However, the scope of the present work is to find a single optimal embankment height and to study how the consideration of hydraulic system behaviour affects this choice. Thus, W is assumed to be equal to zero. The expected annual damage (EAD) is computed per location as: where L is the flood damage (€); H is the water level in the river (m + m.s.l.), with H min being the lowest water level causing flood damage and H max is the water level corresponding to the 12,500-year return period event (maximum conceivable event in this study); u is the embankment height; p(H) is the exceedance probability of a given water level H. The EAD represents the average annual flood damages that communities would expect over the entire planning period. However, not all the expected damages have to be valued equally. Typically, the farther in the future losses are, the less important they are considered by policy makers. For this reason, EAD is discounted as follows: where EAD d is the discounted EAD over the planning period T and r is the discount rate. Typically, cost-benefit analysis studies in the Netherlands apply a discount rate of 4.5%. However, our analysis neglects the effect of future economic growth on damage development and, as a consequence; a lower discount rate is needed to make investments in the increase of embankments height a viable solution. We thus, assume a fixed discount rate equal to 1.5% over a planning period of 50 years. Finally, the optimization problem can be formulated as follows: where u is the decision variable, that is, increase in the embankment height, which can take values ranging from 0 to 1 m with an interval of 10 cm. Equation (4) introduces an optimization problem where total costs at each location have to be minimised. If, as in current practice, hydraulic system behaviour is not considered, this optimization problem can be solved for each location separately. However, accounting for hydraulic system behaviour requires solving Equation (4) as a many-objective optimization problem. This is because of two main reasons. First, when acknowledging hydraulic system behaviour, a given set of optimal solutions no longer depends on the single decision objective of many locations, but rather on the many decision objectives of the entire embankment system. In such a system, flood risk reduction at a given location does not solely depend on measures implemented at that location, but it may also be accomplished by acting elsewhere. Second, the European Flood Directive (Directive 2007/60/ EC, 2007) prescribes, founded on the solidarity principle, that flood risk management plans 'shall not include measures which, by their extent and impact, significantly increase flood risks upstream or downstream'. This can only be achieved if trade-offs among locations are explicitly taken into account. Finally, a many-objective optimization approach allows a posteriori decision support, the aim of which is not to dictate the adoption of a given solution, which may be biased by upfront specified preferences of decision objectives or from their premature aggregation, but rather to support discussion and provide policy makers with a set of reasonable choices. In this way, a policy maker having a certain preference about the decision objectives (e.g., by assuming each location as equally important and thus looking for a least system-wide total costs solution), can still aggregate them accordingly.

| MANY OBJECTIVE ROBUST DECISION MAKING
Given the inherent uncertainties related to hydraulic system behaviour and the many-objective nature of flood risk management planning aiming to properly account for it, we solve the problem in Equation (4)  1. Policy problem formulation. This step requires determining which system elements and decision objectives are important and should be included in the simulation model.

Generating alternatives. This step employs Many
Objective Evolutionary Algorithms (MOEAs) (Coello Coello et al., 2007) to find a Pareto-approximate set of solutions, namely solutions for which it is impossible to improve a single objective without deteriorating the performance of at least one other objective, relative to a reference situation. Thus, by providing the best approximate set of Pareto optimal solutions, MOEAs allow displaying critical trade-offs emerging from alternative policies without a priori attributing preferences (or weights) to any of the decision objectives. The success of any MOEA in finding an approximation of the Pareto front is measured according to the convergence (the evolution of the Pareto front) and diversity (degree of distribution of the solutions over the entire Pareto front) of the solutions. In this study, the ε-NSGAII search algorithm is used (Kollat & Reed, 2005), which exploits adaptive population sizing to provide more diverse solutions during the search. 3. Uncertainty analysis. In this step, the previously found Pareto optimal solutions are stress-tested under uncertainty to evaluate performance across a wide range of scenarios. Ideally, in this step, all previous assumptions related to, for example, parameter values, model structure and problem formulation are relaxed. 4. Scenario Discovery. In this step, statistical clustering techniques are applied to the large dataset of model output generated in step 3. Scenario Discovery (Bryant & Lempert, 2010;Kwakkel & Jaxa-Rozen, 2016), uses the Patient Rule Induction Method (PRIM) (Friedman & Fisher, 1999) to find orthogonal subspaces in the model input space (i.e., the space spanned by the uncertain factors) for which the resulting output is substantially different from the typical model output. A subspace is described by subintervals for one or more uncertain factors. PRIM returns a series of increasingly smaller subspaces. This series presents a trade-off between coverage (percentage of the cases of interest captured by a given box) and density (number of cases of interest over the total number of cases of a given box). Users can select their preferred subspace based on this trade-off.
The analysis is carried out through the Exploratory Modelling and Analysis Workbench (EMA-Workbench) (Kwakkel, 2017), an open source toolkit developed in the Python programming language.

| UNCERTAINTIES AND THE SIMULATION MODEL
Three are the uncertain factors considered, they relate to (see Figure 2): 1. Embankment strength. This is represented by the conditional failure probability sampled from the fragility curve, that is, a curve that indicates the embankment's probability of failure given a water level. The lower the sampled conditional failure probability, the higher the water level that would cause failure (i.e., critical water level), the stronger the embankment; This study makes use of fragility curves used in policy support flood risk management studies in the Netherlands; 2. Breach growth. Growth is assumed to follow an exponential model; however, the considered growth rate changes substantially. Three models are possible, that is, where the maximum breach width is reached after 1, 3, or 6 days; 3. Maximum breach width, B max . Final breach widths can assume values between 30 and 350 m.
The above uncertainties apply to each of the five locations, resulting in a total number of 15 uncertainties. Possible values of the identified uncertainties are summarised in Table 2. Figure 3 provides a schematization of the simulation model and the required modelling steps. The following subsections describe each step in details. The model is entirely implemented in the Python programming language.

| Pre-processing
The pre-processing involves (a) calibration of the routing scheme and (b) adjustment of fragility curves to the flood protection standards in place.
F I G U R E 2 (a) Breach growth models. (b) Adjusted fragility curves per location Flood routing is modelled by applying a Muskingum method (Todini, 2007). In the Muskingum method, the flood wave is routed by solving the continuity equation between the upstream and downstream ends of a river reach while aggregating any geomorphological and hydraulic characteristic into two parameters (Todini, 2007). In our model, we apply a Muskingum scheme for each two subsequent locations in Figure 1, which thus requires four different sets of calibrated parameters. Calibration is performed against the results of the SOBEK, a 1D hydraulic model, by following the method of least-squares as in (Karahan, 2012).
Fragility curves represent the probability of failure of an embankment section given a hydraulic load. The study focuses on overtopping breaching mechanisms triggered by high water levels. Indeed, neglecting the effect of hydraulic loads, such as water level duration and failure mechanisms, such as piping is a simplification. This implies breaches are only possible on the ascending limb of the hydrograph leading to an overestimation of the flood attenuation effect after a breach. Fragility curves are adjusted in a way as to comply with the actual flood protection standards (in this case, 1/1250 for each stretch) by applying Crude Monte Carlo as F I G U R E 3 Schematic representation of the flood risk system simulation model  Diermanse, De Bruijn, Beckers, and Kramer (2015). For each stretch, given N realisations of water levels and critical water levels, the probability of failure can be defined as follows: where Z is the limit state function, defined as the difference between critical water levels (strength) and water levels (load), I is one when Z is negative, and zero otherwise. The process of adjusting the fragility curves entails iteratively shifting them (thus, changing critical water levels) until pf equals the target failure probability.

| Event generation
An event is defined by the flood hydrograph, the conditions of the embankments system (embankment strength and the breach growth dynamic) and the adopted interventions (embankment height). Flood waves are generated by associating a sampled maximum discharge to a normalised hydrograph. Maximum discharges are generated following a Generalised Extreme Value distribution Type I, that is, a Gumbel distribution, as in De Bruijn et al. (2014) and Diermanse et al. (2015): where coefficients a and b are equal to 1,316.45 m 3 s −1 and 6,612.5 m 3 s −1 , respectively. Such parameters are found by fitting the distribution to high discharges at Lobith, that is, where the Rhine River enters the Netherlands. The Rhine then bifurcates into three branches, one of which is the IJssel River. High discharges at Lobith and the corresponding flood levels at an upstream location on the IJssel River have been compared using SOBEK. On average, SOBEK simulates the IIssel River as discharging about 15% of the water that enters at Lobith. We used this figure to adapt Equation (5) to find a distribution of maximum discharges for the IJssel River. Expected annual damages are calculated based on 100 upstream high discharges. A flood hydrograph is then generated by multiplying the sampled maximum discharge with one of the plausible normalised hydrographs calculated for Lobith in the GRADE project (Generator of Rainfall and Discharge Extremes) (Hegnauer, Beersma, van den Boogaard, Buishand, & Passchier, 2014). The sampled 100 upstream discharges and the normalised hydrograph are reported in Figure 4. Finally, values of the exogenous uncertainties (embankment strength, breach growth rate and maximum breach width) and intervention (embankment height increase) are sampled.

| Event simulation
As the flood hydrograph propagates through the river channel, failure probabilities and impacts on the hydrograph are evaluated at every location of interest. At each location, (a) the coming discharge is translated into water levels through a stage-discharge relationship and (b) the water level is compared with the critical water level. For each model run, water levels are compared with sampled critical water levels at each location and for each time step. If the water level exceeds the critical water level, a breach is simulated using a weir formula and water will flow from the river into the protected area. Downstream discharges will then become lower. When considering hydraulic system behaviour, this discharge is hence subtracted from that in the F I G U R E 4 (a) Sampled upstream high discharges and associated return period. (b) Normalised hydrograph main channel. Conversely, when hydraulic system behaviour is ignored, the discharge in the main channel is held constant.

| Damage estimation
Losses are estimated by comparing the modelled water levels with water level-damage functions. These functions are based on results of the VNK project (in Dutch: Veiligheid Nederland in Kaart, in English FLORIS: Flood Risks and Safety in the Netherlands) (Jongejan & Maaskant, 2015). VNK is a major flood risk analysis project which provides, for the areas considered in the present study, damage estimates for the 1:125, 1:1250, and 1:12,500 per year floods.

| Generating alternatives assuming a reference scenario
As shown in Table 2, in the reference scenario embankment failure is supposed to occur as soon as the water level equals the one corresponding to the 0.5 failure probability and, when embankments fail, the breach will grow up to 175 m after about 3 days.
Under the reference scenario, the optimization problem as defined in Equation (4) is solved twice, namely by accounting for hydraulic system behaviour and neglecting it. Results are presented as a parallel plot in Figure 5 where each single line represents a solution.
For a fairer comparison of the outcomes, all optimal solutions in Figure 5, regardless whether they account for or ignore hydraulic system behaviour, are re-evaluated accounting for hydraulic system behaviour. This is because, ultimately, reality is such that upstream breaches do cause some degree of flood attenuation and the performance of each solution must therefore be evaluated accordingly. The question is whether ignoring this phenomenon during the design phase will lead to differences in the identified solution(s). When neglecting hydraulic system behaviour, there is a single unique optimal solution that minimises the total costs across all locations (blue line in Figure 5). Instead, when hydraulic system behaviour is taken into account, a set of 17 Pareto optimal solutions is identified (light orange lines in Figure 5). Thus, neglecting hydraulic system behaviour may lead decision makers to a solution based on cognitive myopia (Hogarth, 1981), that is, a situation in which the problem formulation is too narrow and possible alternative courses of actions remain hidden. The wider set of solutions illustrates that trade-offs exist between locations when deciding upon optimal increases of embankment height. Some solutions may lead to very low total costs upstream, by for example, raising upstream embankments while neglecting the downstream hydraulic load (thus risk) transfer; F I G U R E 5 Parallel plot of optimal solutions. The first five vertival axes indicate decision variables (i.e., degree of embankment rasing), the subsequent five axes indicate decision objectives scores (i.e., total costs) and the last axis indicates the system-wide total costs (i.e., sum of each decision objective score). The solution in blue represents the optimal solution found when hydraulic system behaviour is neglected. solutions in light orange represent the approximate Pareto set of solutions found when hydraulic system behaviour is taken into account. Of these latter solutions, solutions D and S are depicted in bold as a dotted and continuous line, respectively while other solutions may lead to very low downstream costs, by for example, keeping upstream embankments low.
Overall, the single optimal solution found when ignoring hydraulic system behaviour represents a conservative choice of optimal embankment heights because it requires the highest embankment heights among all solutions. A policy maker could still choose to be conservative while acknowledging hydraulic system behaviour by, for instance, opting for solution D in Figure 5. Interestingly, these two conservative solutions are the same for all stretches except for the most downstream one, A.5. For this stretch, solution D can suffice with a lower embankment height and, also, incurs lower total costs than the optimal solution that is found when ignoring hydraulic system behaviour. Consequently, the optimal solution found neglecting hydraulic system behaviour is Pareto dominated by solution D.
Solution S represents the overall least cost solution, where embankments are raised at stretches A.1 and A.3 only. This leads to a situation in which location A.2 bears the greatest burden in terms of total costs and, by so doing, lower total system costs can be achieved. Thus, the optimal solution identified while neglecting hydraulic system behaviour is in fact Pareto suboptimal with respect to solution D and suboptimal from a system perspective with respect to solution S. Figure 6 helps investigating this further. Figure 6a shows the expected annual damage, investment costs, and total costs at stretch A.5 as a function of the embankment height. The embankment heights of the other stretches are held constant at the level required by both solution D and the optimal solution found when neglecting hydraulic system behaviour. When hydraulic system behaviour is taken into account, a lower expected annual damage (blue dotted line) results from the hydraulic load reduction because of upstream flooding. This reduced estimated expected annual damage causes a reduction of total costs (green dotted line) as well as a shift to the left of the minimum of the total cost curve, which moves from an optimal embankment raising of 0.5-0.3 m. Figure 6b shows the system total costs (i.e., the sum of all total costs) as a function of the degree of the total embankment height increase (i.e., the sum of all embankment height increase). In the figure, dots represent all solutions in terms of expected annual damage, investment costs and total costs. A line has been fitted through the points to obtain an approximated system-wide function through these three cost measures. Solution S is by far the system-wide least total costs solution and may be regarded as an outlier. However, even when neglecting solution S, the overall optimum is reached at a total degree of embankment height increase of approximately 1.1 m, which is lower than the 1.5 m required by the optimal solution found when hydraulic system behaviour is neglected.
F I G U R E 6 Change in expected annual damage (blue), investment costs (red) and Total costs (green) for increasing embankments heights for location a.5 (a) and the system as a whole (b)

| Uncertainty analysis
In this section, all solutions previously found are evaluated under multiple scenarios. In this context, solutions are set of embankment heightening, while scenarios analyse the uncertainty of the solutions' performance by varying and combining parameters in the ranges given in Table 2. Range of values for each type of uncertainty and average values used in the reference scenario. The overall aim of this section is to test the robustness of the previous findings, namely the conditions of sub-optimality of the solution found neglecting hydraulic system behaviour. This is accomplished by re-evaluating the decision objectives of each solution under uncertainties about the probability of failure, the final breach width and the breach growth rate. For the same reason stated in the previous section, all solutions are evaluated accounting for hydraulic system behaviour. A Latin Hypercube Sampling (LHS) technique is adopted to generate 2000 scenarios, from the value ranges specified in Table 2.
The aim of the analysis is to study the robustness of each solution (a) in retaining Pareto optimality and (b) with respect to the system-wide performance. In doing so, two different robustness metrics are adopted. For a review of robustness metrics, the reader is referred to McPhail et al. (2018).
Robustness in retaining Pareto optimality is measured with a satisfying robustness metric, that is, metrics that evaluate the range of scenarios having an acceptable performance. In this case, the acceptable performance relates to the capability of a solution to be among the set of Pareto solutions over all scenarios. Thus, for each solution, the likelihood of being among the Pareto set, that is, the ratio between the number of scenarios in which a given solution is Pareto dominant and the total number of scenarios, is calculated. F I G U R E 7 (a) Likelihood to retain Pareto optimality (y-axis) and mean regret relative to system-wide performances (x-axis) of all solutions.
(b) Boxplots of the regret relative to system-wide performances of the trade-off solutions T A B L E 3 Statistics of the regret to system-wide best performance of the trade-off solutions in Figure 7 Likelihood to retain Pareto optimality Robustness with respect to the least system total costs solution is measured with a regret-based robustness metric, which is defined as follows: Regret P j, s = P best,s −P j,s j = 1,…, n f g,s = 1,…,m f g ð6Þ where n is the number of solutions and m is the number of scenarios. In other words, for a given scenario s and a solution j, regret is defined as the difference in performance P (i.e., system total costs) between the solution j and the best performing solution, that is, the one resulting in the system's least total costs. Results are shown in Figure 7. In Figure 7a, all solutions are plotted in terms of their likelihood to retain Pareto optimality and the mean regret to the best system-wide performance. An ideal solution would be found in the bottom-left side of the plot, that is, having high capability of retaining Pareto optimality and low regret. Interestingly, none of the solutions perform in this way. Rather, a set of trade-off solutions can be identified, where improvement in the likelihood to retain Pareto optimality comes at the expense of regret to system-wide best performance, and vice versa. Hence, the uncertainty analysis allows excluding from the set of plausible final F I G U R E 8 Scenario discovery results for cases when the solution found neglecting hydraulic system behaviour does not (CASE A) and does (CASE B) retain Pareto optimality. Each row shows two graphs relative to (1) the density-coverage trade-off trajectory resulting from PRIM's boxslicing process (left) and (2) which uncertain factors and the associated range of values best describe the cases of interest (right). Boxes with a coverage of at least 0.7 and a density of at least 0.8 are selected whilst non-significant values (i.e., values whose quasi P-value is higher than 0.05) are dropped F I G U R E 9 Mean regret relative to system-wide performances for increasing degree of the total embankment height increase. Results for solutions where hydraulic system behaviour is accounted for or neglected are depicted in orange and blue, respectively. When more than one solution have the same total degree of embankment height increase, a black vertical line is shown representing the error bars with a 95% confidence interval solutions those outside the trade-off curve. This would not have been possible relying solely on the analysis under the reference scenario. Interestingly, solutions D and S are not among the plausible set of solutions. Moreover, the optimal solution found while neglecting hydraulic system behaviour, which is Pareto dominated under the reference scenario, does belong to this trade-off set of solutions. It represents the one that is least capable of retaining Pareto optimality and also the one that yields the least regret with respect to system-wide best performance. In that sense, it qualifies as a 'better safe than sorry' solution (Klijn, Asselman, De Kruif, Bloemen, & Haasnoot, 2016). Figure 7b shows the distribution of regret to system-wide best performance of the trade-off solutions. The related statistics are shown in Table 3.

| Scenario discovery and trade-off analysis
Scenario Discovery is used to investigate two distinct types of outcome: 1. The solution found neglecting hydraulic system behaviour does not retain Pareto optimality; 2. The solution found neglecting hydraulic system behaviour does retain Pareto optimality.
The results of this investigation are shown in Figure 8. Case a is explained by scenarios where stretch A.5 is relatively strong. The interpretation is rather straightforward: when downstream stretches do not fail easily, a solution that overdesigns such embankments, as the one found neglecting hydraulic system behaviour, proves sub-optimal.
Case b is complementary to a. It occurs when stretch A.5 is weak whereas upstream stretches, for example, A.4 and A.1, are relatively strong. Thus, when subject to flood events higher than the design flood, the non-failure of upstream stretches implies a transfer of loads downstream, where the embankments then fail instead. An overdesign of the downstream embankment stretches allows for better coping with the exacerbation of hydraulic load coming from upstream, thus preserving the Pareto optimality of such a solution.
F I G U R E 1 0 Likelihood to retain Pareto optimality for increasing degree of the embankment height increase per location. Results for solutions where hydraulic system behaviour is accounted for or neglected are depicted in orange and blue, respectively. The black vertical lines are error bars representing the 95% confidence interval The trade-off between Pareto dominance and regret to system-wide best performances is explored by looking at the relationship of these two measures with the degree of embankment height increase. In particular, Figure 9 shows that the mean regret decreases when increasing total degree of embankment height. This explains why conservative solutions, like the optimal solution found while neglecting hydraulic system behaviour, are among the ones leading to the least mean regret.
Conversely, Figure 10 shows that the likelihood of retaining Pareto dominance decreases with an increasing degree of embankment height at each location. This makes conservative solutions, like the optimal solution found while neglecting hydraulic system behaviour, to be among the worst performers in retaining Pareto optimality.
Under uncertainties on how, where and to what extent flood attenuation effects will take place, an approach that finds optimal embankment heights while neglecting them guarantees a low regret to system-wide total cost. At the same time, however, because of this neglect, the solution is not capable of guarantying optimality at each stretch, where lower investments would have been justified had flood attenuation been taken into account. Ultimately, the choice of a final plan is part of the decision-making process and it depends upon the preferences of decision makers with respect to the two decision robustness criteria. However, one could hypothesize what a reasonable solution would be under different contexts.
For example, under the assumption of centralised funding for flood risk management, policy makers may opt for the solution that neglects hydraulic system behaviour if they are confident that overspending downstream will not outweigh the saved total costs of the whole system. If funding is instead not centralised and a negotiated solution must be agreed upon, one that better retains Pareto optimality is more likely to bring consensus among parties.
Finally, one must be aware that considering economic risk only is too limited. Decisions should be made including a wider set of risk measures, for example, societal and individual risk. De Bruijn et al. (2014) found that societal risk is overestimated when hydraulic system behaviour is not accounted for. Our results are in line with that, yet, because of the associated uncertainty, the question of 'how safe is safe enough' still holds, especially when human lives are involved.

| CONCLUSIONS
In the present paper, we investigated the effect of ignoring hydraulic system behaviour (i.e., the change in the hydraulic loads at one location as a consequence of embankment breaching at other locations) on making decisions about optimal embankment heights, exemplified by a case analysis of the IJssel River in the Netherlands. For the analysis, we applied the Many-Objective Robust Decision-Making framework.
Current practice in flood risk management often ignores hydraulic system behaviour in establishing optimal embankment heights; this leads to a single optimal solution, that is, a unique set of optimal embankment heights. In contrast, taking into account hydraulic system behaviour in the optimization problem widens the solution space substantially. Instead of finding one single optimal solution, our case revealed a Pareto set of 17 solutions. Current practice thus leads decision makers to a solution based on cognitive myopia, that is, a situation in which the problem formulation is too narrow and possible alternative courses of actions remain hidden. Furthermore, under conditions of perfect knowledge about the behaviour of the embankment system, the solution that would qualify as optimal according to current practice proves in fact sub-optimal with respect to both downstream locations (i.e., it is not Pareto dominant) and with respect to the system as a whole (i.e., it is not the system-wide least total costs solution). These two conditions of sub-optimality were investigated further under uncertainty.
The uncertainty analysis revealed a trade-off between the ability to retain Pareto optimality under uncertainty on the one hand and the regret with respect to the least total system costs on the other hand. The optimal solution found following current practice is among the trade-off solutions. It shows the best system-wide performance, and also the least capability of retaining Pareto optimality. This lack of Pareto optimality is attributable to its poor performance on the most downstream embankment stretch, where money is spent unwisely, which is in line with what we found under the reference scenario. The good system-wide performance is instead attributable to the conservative nature of such a solution, where embankments are raised following a worst-case approach, that is, neglecting the flood attenuation effects of possible upstream breaching. In other words, it rigorously applies the precautionary principle of 'better safe than sorry', at extra costs.
The modelling results suggest that policy makers willing to pursue flood risk management of large-scale systems, while considering risk transfers that take place within such system, must account for hydraulic system behaviour. This would make explicit decision conflicts arising among the parties involved, thus allowing, as demanded by the EU Flood Directive, a due consideration of fairness in the decision-making process. However, the proposed approach becomes computationally unfeasible if either applied to a substantially higher number of stretches or a more detailed and computationally demanding simulation model is employed. The proposed modelling framework can thus be used for identifying interesting solutions to be explored and tested further with more detailed models.
Further research will focus on developing a flood risk management plan while including a broader set of possible flood risk management strategies, societal and economic risk as well as other failure mechanisms. Furthermore, considerations of equity between stretches will be given thorough attention. A risk transfer decision objective will be formalised with the aim of studying its effects on the attractiveness of the solutions.