Evaluation of flood risk reduction strategies through combinations of interventions

Large, complex coastal regions often require a combination of interventions to lower the risk of flooding to an acceptable level. In practice, a limited number of strategies are considered and interdependencies between interventions are often simplified. This paper presents the Multiple Lines of Defence Optimization System (MODOS)‐model. This quick, probabilistic model simulates and evaluates the impact of many flood risk reduction strategies while accounting for interdependencies amongst measures. The simulation includes hydraulic calculations, damage calculations, and the effects of measures for various return periods. The application and potential of this model is shown with a conceptual and simplified case study, based on the Houston‐Galveston Bay area. The analyses demonstrate how the MODOS‐model identifies trade‐offs within the system and shows how flood risk, cost, and impact respond to flood management decisions. This improved understanding of the impact of design and planning choices can benefit the discussions in finding the optimal flood risk reduction strategy for coastal regions.

A storm surge can heavily affect a coastal region, causing both direct and indirect damage to people, structures, and the environment. The impact of storm surges can be reduced by implementing flood risk reduction measures, in the form of structural interventions like levees and storm surge barriers, or non-structural measures such as wetlands and oyster reefs. However, coastal zones that combine multiple coastal morphologies, like dunes, barriers islands, bays, and estuaries, often require a combination of flood risk reduction measures, which complicates the search for the best strategy.
For flood protection systems purely consisting of structural elements, several studies have used probabilistic risk analysis to elaborate on the interdependence-how one measure affects the performance of the other-between structures (Courage, Vrouwenvelder, van Mierlo, & Schweckendiek, 2013;De Bruijn, Diermanse, & Beckers, 2014;Tsimopoulou, 2015). These studies have shown that a flood risk reduction measure can strongly influence the effectiveness of other elements of the flood defence system. This is true for both defences placed adjacent, protecting the same area (resembling a series system, which is as strong as the weakest link), and for defences placed in multiple lines of defence, like a parallel system. This interdependency amongst measures complicates the assessment of the risk reducing abilities of the overall flood protection system. Besides flood defence structures, non-structural measures like flood risk zoning and flood proofing of local structures can also contribute significantly to reduce flood risk (Aerts et al., 2014;Dawson et al., 2011;Kreibich, Thieken, Petrow, Müller, & Merz, 2005). Additionally, Nature-based Solutions can be considered, which provide both risk reduction and ecological value. It is important to take into account these different interventions (and their combinations) in the development of flood risk reduction strategies for flood-prone areas. A key question is which combination of these interventions is the most effective.

| Literature and knowledge gaps
Although the optimisation of risk reduction provided by single structures has been investigated extensively in the past (Eijgenraam et al., 2014;Kind, 2011;Tung, 2005), only specific combinations of two interventions (levees, hard structures, natural foreshores) have been analysed by hand (Courage et al., 2013;Vuik, Jonkman, Borsje, & Suzuki, 2016). Common techniques for mathematical optimisation (i.e. mixed integer or dynamic programming) fail for these cases, as the flood defences cannot be viewed separately. Therefore, studies on this subject often assume total independence (Kind, 2014) or focus on small-scale decisions and processes, making it harder to represent complex systems (Dupuits, Diermanse, & Kok, 2017).
The use of multiple lines of flood defence in flood management strategy was introduced under various names (Tsimopoulou, 2015). Examples are "hierarchical flood protection system" (Custer, 2015) or "multiple lines of defence strategy" (Lopez, 2009), amongst others. Lopez (2009) focussed on the inclusion of both structural and non-structural measures, but offers a qualitative model, where results are harder to reproduce or compare. Instead, Custer (2015) used hierarchical probabilistic modelling to quantify the impact of multiple structure in a parallel system. Although this does result in quantitative risk analysis, flood risk management is often part of a broader context-i.e. inclusion of non-structural measures-which cannot be captured with this method.
Risk analysis on the scale of regional floodplain systems has been described in earlier research (Aerts et al., 2014;Gouldby, Sayers, Mulet-Marti, Hassan, & Benwell, 2008;Moser, 1996;Woodward, Gouldby, Kapelan, & Hames, 2013). Similar to this paper, Woodward et al. (2013) attempted to combine simplified hydraulic modelling with common optimisation techniques for comparison of flood management strategies. Although they managed to gain satisfactory results for systems with multiple defences placed in a series system, their approach is not applicable to systems with multiple lines of defence. Current approaches do not capture the complexity of implementing and combining different interventions and their effects on risk and other dimensions. To address these shortcomings, this paper presents an approach combining quantitative risk assessment of interdependent structural and non-structural interventions (e.g. Nature-based Solutions, spatial planning, and disaster management) with multi-layer flood protection.

| Objectives and scope
Decision making in flood risk management would greatly benefit from a clearer understanding of the effects of a portfolio of interventions on flood risk throughout a coastal region. However, probabilistic flood risk assessment for combinations of flood defence structures can be computationally intensive and time consuming. This paper presents the Multiple Lines of Defence Optimization System (MODOSmodel) as a viable model for the evaluation of coastal defence strategies. MODOS is a fast risk-based model that simulates and evaluates the impact of many alternative flood risk reduction strategies, each consisting of a combination of measures. The computation time is significantly reduced, as basic (hydraulic) formulas are used for the flood and damage simulation, which allows for many strategies to be compared.
The main characteristics of the model include the ability to (a) consider both structural and non-structural measures, (b) allow for economic a`nd non-economic performance indicators, and (c) have a generic setup, which is easily adapted to other flood-prone regions around the world. In this paper, the model will be demonstrated with a conceptual and simplified case study in the Houston-Galveston Bay area in Texas.
The main purpose of the model is to support decision making in the early phases of design. Especially when large, complex flood-prone regions are concerned, many planning decisions are required in these early design stages, while only limited information is available. This method therefore focusses on the simulation and evaluation during initial screening or the conceptual design phase. The assessment of flood risk reduction strategies is based on a multiple-objective approach. The economic assessment compares the construction costs of a combination of interventions with its ability to reduce flood risk. Here, risk is defined as a combination of scenarios, each of which has a probability of occurrence and a potential negative consequence (Kaplan & Garrick, 1981). The consequence includes damage to the region itself as well as damage to the flood protection structures. Besides the economic assessment, the flood risk reduction strategies will also be evaluated separately based on non-economic impact criteria such as environmental impact. These impacts will be ranked using simplified indicators. The model evaluates the flood risk reduction measures, which allows the model to compare strategies in terms of costs, risk reduction, and other impacts. The MODOS-model can simulate and evaluate the impact of numerous flood risk reduction strategies in a large, complex region. MODOS is implemented in Python and consists of two parts: the simulation model and the evaluation model ( Figure 1).
The simulation model combines three layers of information for its simulation: data on the spatial layout of the region, the flood risk reduction strategy (which is a chosen combination of interventions), and the incoming storm (represented as hydraulic boundary conditions). This is shown in Figure 2(left).
In the simulation model, the study area itself is schematised as a combination of lines of defence and protected areas (see Figure 2, right). The lines of defence are the locations where structural flood defences can be placed to retain water. Here, calculations on failure probability, construction cost and discharge are made. The protected areas consist of basins between lines of defence, for which water levels and damages are calculated. In the schematisation of Figure 2, defence A1 (land barrier) and A2 (storm surge barrier) are located adjacent. Damage to one will hardly change the load on the other. Barriers A and B are in different lines of defence. Damage to the barriers A1 or A2 will greatly affect the loads on B1 and B2. The flood risk can be grossly underestimated when this effect is not considered.
The simulation model includes two types of calculations: the hydraulic calculations and the damage calculations. The hydraulic calculations result in inundation levels across the region, which the damage calculation uses to calculate the amount of damage caused by the flood. To speed up the calculations, important characteristics of the measures and the region are determined beforehand. For example, the strength of a flood protection structure is defined by a fragility curve, which is a graphical representation of how the failure probability of a structure depends on its main load (US Army Corps of Engineers [USACE], 1996). Likewise, the inflicted damage in the affected basins is linked to the local water level in the form of a damage curve, which will be derived based on local damage estimates and land use data.
The evaluation model uses the simulations to compare numerous different strategies. First, it repeats the simulation for different storm intensities to build a "risk curve" (Figure 3). This shows the potential inflicted damage for an event and the probability of that event. By comparing this risk curve with the risk curve of the initial local situation, the risk reduction can be computed.
Subsequently, this risk curve is built for numerous strategies. This is possible because of the short computation time (i.e. 10 s on a consumer machine for a single storm). The number of different strategies required to sufficiently explore the design space depends on the regional complexity and the  amount of potential measures. The resulting data is stored and used for further analysis and evaluation. Currently, the model is presented for present conditions that is, for a static situation. In future versions, the lifetime of measures (e.g. degradation) and potential future scenarios (e.g. SLR and subsidence, economic growth) will be included as factors. Using the current model, the effect of future changes can be assessed by means of a sensitivity analysis.

| Schematisation
As the flood progresses into the project area, it encouters lines of defence and protected areas, one by one (Figure 4). When the flood encouters a line of defence, MODOS calculates the failure probability and the discharge into the next FIGURE 3 An often-used version of the risk curve is the frequencydamage curve. This shows the probability of exceedance of damage. The dashed line represents the risk curve in the current situation (without additional risk reducing interventions). Implementing flood risk reduction measures aims to lower the probability of experiencing damage or to lower the amount of damage resulting from a flood. As a result, the curve shifts towards the lower-left and the total risk decreases FIGURE 4 Schematisation of how different outcome scenarios are taken into account. In the example region (up), there are two lines of defence with each two flood defence types. The chart (down) shows how the expected value of damage can be calculated from the different possible outcome scenarios area. The protected areas in between are divided into several basins, based on watersheds. Here, the water levels in the basins and the hydraulic conditions going forward are calculated. The expected impact of one storm depends on whether the implemented flood risk reduction interventions will hold or fail. Given a few flood risk reduction measures, scenario outcomes can differ significantly within the same strategy. This is considered by simulating all scenarios. Figure 4 shows a schematised example, where each line of defence consists of two flood defences, each of which can fail or hold. This system can lead to 16 potential scenario outcomes. After the simulation, the expected values of damage for each scenario are weighted according to their probability and summed to result in the expected damage for that flood risk reduction strategy for one storm.

| Flood risk reduction measures
A flood risk reduction strategy consists of a combination of potential measures. It is a choice between several locations, effects and (if applicable) barrier heights. The location and effect varies based on the type of measure. Therefore, they are divided into four different categories as shown in Table 1.
The effect of flood risk reduction measures differs between structural and non-structural measures. Structural measures are mostly protective measures such as barriers and levees and reduce the likelihood of flooding. Nonstructural measures could affect the risk in several ways. Measures like Nature-based Solutions are often not sufficiently effective by themselves. Economic benefit and safety provided by these measures are hard to quantify, while they often offer co-benefits like ecological value or a better living-environment (Vuik et al., 2016). These types of measures are mostly used to reduce waves for damage reduction (Mazda, Magi, Ikeda, Kurokawa, & Asano, 2006) or to reduce wave impact on nearby flood defences as a part of a hybrid solution (Vuik et al., 2016). Non-structural measures can thus also be incorporated by modelling their effect on different components of risk.

| Simulation model 2.2.1 | Hydraulic calculations
The simulation of the storm revolves around the water flow through the region. The storm surge encounters either a line of defence or a protected area. The hydraulic calculations consist of three relations that are used throughout: Within the lines of defence, the model calculates the failure probabilities of the flood protection structures and the discharge flow into the next layer. Within the protected areas, water levels in the basins are calculated. Appendix shows pseudo-code for how these choices are made within the model. It also shows the common hydraulic formulas used in the model.
The failure probability is calculated by comparing the incoming storm surge water level with the fragility curve of the flood defence structure. A fragility curve is a graphical representation of the conditional probability distribution (Schultz, Gouldby, Simm, & Wibowo, 2010;Van der Meer, Ter Horst, & Van Velzen, 2009). It represents the strength of the structure by showing how the probability of failure depends on one type of loading (see Figure 5). To limit the model's computation time, the strength related calculations for the fragility curve are performed beforehand. Without additional knowledge of the structure, a fragility curve can be used in the form of a cumulative normal distribution. Because a fragility curve only takes one load into account, other loads (e.g. waves) must be included when the fragility curve is made. In this model, additional loading will affect the μ and σ in (1): where Pf is the failure probability; h is the water level (m); μ is the mean value. Based on flood defence type and incoming wave height (m + MSL); σ is the standard deviation. Based on flood defence type and incoming wave height (m). Surge flow to the next layer is calculated using simple hydraulic relations. In this model, there are three forms of flow: Overtopping/Overflow, Broad Weir and Open channel flow, calculated with formulas from EurOtop (2016), Brunner (1995), and Stoeten (2013), respectively. Table 2 shows which situations are possible and how they are considered.  The water level on the inside of the line of defence changes based on the discharge. In the model, the storm surge is represented by a time series. For every time step, the inside water level can be found with where h p (i), h p (i + 1) is the water level in the area protected by the line of defence at time step i, i + 1 (m); Q(i) is the sum of all discharges across the flood protection (m 3 ); Δt is the length of the time step in seconds (s); A p (i) is the surface area protected by the line of defence (m 2 ). The maximum water levels in the watershed basins can be found by focussing on the protected areas between the lines of defence. The discharge calculated at the line of defence fills the area, leading to a time series of the water level at the inside. Here, the maximum water levels in all basins are calculated, taking wind set-up into account (see Figure 6). In earlier research on the Galveston Bay, this method showed to be effective in predicting water levels resulting from past storms (Dupuits, Kothuis, & Kok, 2017;Stoeten, 2013). The wind set-up peaks in the direction of the wind. The area is divided into eight wind directions. The primary wind direction is included as a time series. With the centre of the inundated area as reference point, using the deviation of every basin from the main wind direction, the water level in each basin is calculated: where h b is the maximum water level in the watershed basin (m + MSL);h p is the maximum average water level in the protected area (m + MSL); S is the wind set-up (m); C b is the wind direction factor (between − 0.5 and 0.5 ). This calculation also plays a role in the transition to the next layer. Hydraulic conditions at the following line of defence depend on its relative position to the protected area. For example, wind set-up will lead to different water levels High inside water levels Broad weir/overtopping Backflow for flood defences placed in different wind directions as seen from the source of the flood. The surge at this next line of defence is built by scaling the time series of the water level in the protected area. The maximum significant wave height is highly influenced by local conditions and therefore hard to predict. For conceptual design, a generally used rule-of thumb is to assume the significant wave height to be half of the water depth (Molenaar & Voorendt, 2018).

| Damage calculations
Damages are calculated by finding the inundation level, as a function of the difference between the maximum water level and the land elevation. To increase accuracy, every basin is divided into height contours, which consists of several land use types (see Figure 7). Flood height(h f ) is inserted in the damage curve to find the portion of property damaged for one land use type. This is divided in damage to the structure and damage to the content. This is summed for every land use type {1, 2, …, u} to find the damage for one contour (D T ). In turn, this is summed across all contours {1, 2, …, n}, all watersheds {1, 2, …, m} and weighted for the probability of each outcome scenario {1, 2, …, s} to result in the expected damage for one flood risk reduction strategy(D s ): where D s is the expected damage for one flood risk reduction strategy ($); D T is the expected damage for one land use Process of calculating damages for a given flood level. Every basin is divided into height contours. Subtracting the elevation from the water level gives the maximum inundation. This is input for the damage curves, which estimates damage for one land use type ($); P s is the probability of outcome scenario; V st , V ct = aggregate value of one land use type in terms of structures and content, respectively ($); p st (h f ), p ct (h f ) is the estimated portion of value damaged of structures and content, respectively; h f is the flood inundation in contour (m); h b is the maximum water level in watershed (m + MSL), and h cr is the elevation of contour (m + MSL).

| Construction and repair cost calculation
Construction and repair cost are important performance metrics to compare strategies. For every simulation, the construction cost of the measures is estimated, based on local research or reference projects. The cost of construction depends on the type of measure and its length (which is assumed constant). For some measures, height can also be adjusted. This creates additional costs, which are considered as "variable costs." Repair cost of a flood risk reduction measure is included as a percentage of the construction cost and is only added in situations where that measure fails.

| Environmental score schematisation
Large scale coastal interventions can have a significant impact on the environment. For example, the construction of dams in the Netherlands has changed part of the environment in these estuaries (Smits, Nienhuis, & Saeijs, 2006). On the other hand, Nature-based interventions, like nourishments and wetland recreation, may even improve environmental conditions. It is therefore important to include an assessment of these impacts in decision-making. However, no uniformly accepted method is available to evaluate environmental impacts of coastal interventions, leading to previous studies struggling to quantify environmental impact (Kind, 2014;Lopez, 2009). To illustrate how environmental impacts can be included, an environmental score is included to compare between strategies. The score takes two processes into account: • Structural measures negatively impact their direct surroundings. This impact increases with construction height, as a higher structure will have a relatively higher impact on natural habitats. • Inundation of ecologically important areas. Flooding in the watersheds that house nature reserves causes damage and should be avoided. Higher water levels in these watersheds result in a more negative environmental score.
The inclusion of an environmental score allows for comparison between strategies. It is emphasised that this score is highly simplified, not verified and is included to illustrate how environmental factors can be included. A verified approach can be added in later research, which will benefit the broader applicability of the model.

| Evaluation model
The MODOS-model allows us to discover how the risk profile of the region reacts to different flood risk reduction strategies and identify interesting trade-offs. These analyses are done with the use of the "Exploratory Modelling and Analysis (EMA)-workbench" (Kwakkel, 2017a(Kwakkel, , 2017b. It has been used for a variety of research topics in the past (Halim, Kwakkel, & Tavasszy, 2016;Kwakkel & Cunningham, 2016;Kwakkel & Jaxa-Rozen, 2016). In this context, it generates potential flood risk reduction strategies, which are then evaluated using the MODOS model. By systematically sampling potential strategies spanning the entire design space, we can investigate and quantify the impact of design choices.
The EMA-workbench includes a variety of analysis and optimisation tools. Currently, the used techniques are Feature Scoring, Pair Wise Plotting, and Scenario Discovery (Bryant & Lempert, 2010;Kwakkel & Jaxa-Rozen, 2016). Feature Scoring is a form of sensitivity analysis which ranks the input parameters based on their impact on the output parameters. Pairwise plotting visualises the relation between different output parameters. Scenario discovery is a relatively recent approach for identifying subspaces within the model input space that have a high concentration of results that are of interest. The dominant machine learning algorithm for Scenario Discovery is the Patent Rule Induction Method (PRIM, see Friedman & Fisher, 1999). This analysis finds input ranges which meet a set of performance thresholds, which can be defined by the modeller. With these three analysis techniques, it is possible to explore the design-space and investigate the impact of design choices.

| Introduction
The Houston-Galveston Bay is prone to a variety of waterrelated hazards. It is characterised by large variety in land use. The eastern side is mainly occupied by marshlands, while the western side is highly populated with several million inhabitants. Three other notable locations are the low-lying barrier islands, Galveston Island and Bolivar Peninsula; and the Port of Houston. The Port of Houston is one of the most important petrochemical ports in the world. 2 In 2008, Hurricane Ike caused an estimated 30 billion USD in damage, mostly due to coastal storm surge. Moreover, studies have shown that the damage could have been several times higher, should the hurricane have kept its original course 50 km south (Bedient & Blackburn, 2012; Gulf Coast Community Protection and Recovery District (GCCPRD), 2015). In that case, the hurricane would have directly hit Houston and the Port of Houston, which would result in an even more devastating disaster. The MODOS-model compares different strategies for reducing flood risk in the Houston-Galveston Bay area. These strategies consist of combinations of measures shown in Figure 8. Here, the coastal flood protection measures are schematised as the first potential line of defence. The second line of defence is assumed to separate the urban areas of Houston and Texas City from the Galveston Bay. All measures shown are included in the model; although some are combined because of simplifications in the model (see Table 3).

| Model conceptualisation 3.2.1 | Region layout
The layout of the region is based on several spatial datasets, which were used for two primary functions: quantifying losses of a given storm, and estimating additional surface area flooded during each time-step of the hydraulic model. The region is divided into watershed basins, according to the hydrologic unit code (HUC) 10 scale. Each of these watersheds is divided into height contours of 0.3048 m, according to a digital elevation model (DEM), which has a 10-m spatial resolution (see Figure 9).
Both the HUC-10 and the DEM are based on data from the U.S. Geological Survey (Bassler, 2018). The damages are calculated by combining the land use data with damage curves. Type and value of structures on a parcel level are based on data from Harris County Appraisal District (2017). The damage curves are derived from Federal Emergency Management Agency (FEMA) (2010).
Little research has been done on the environmental impact of floods and flood risk reduction measures in the Houston-Galveston Bay. Therefore, this assessment is done on a qualitative basis with an "environmental score" (see Section 2.2.4), where a lower score signals a more negative influence on the local ecology. There are several parts of the region where nature reserves are present and where flooding would negatively impact nature and wildlife (see Figure 10). Table 3 lists the flood risk reduction measures as they are used in the MODOS-model. They are mostly based on reports focussing on the Houston-Galveston Bay area. From the categories mentioned in Table 1, only flood defences and Naturebased Solutions are used. Each measure has been mentioned in earlier research to be useful for reducing flood risk. The MODOS model will create strategies by randomly combining measures and varying their heights (if applicable).

| Hydraulic boundary conditions
The hydraulic boundary conditions are based on incoming storms. The storm enters the model as a time series of water level and wave height at the Gulf of Mexico, as well as a time series of wind velocity and wind direction in the bay (see Table 4). The storm characteristics (i.e. storm track, hydrograph) are based on observations and characteristics of historical hurricanes, analysed by Sebastian et al. (2014) and Ebersole et al. (2017). Because all storms are hurricanes in this case, the wind direction in the bay depends greatly on the hurricane track. In all simulations, a path in north-western direction is chosen, with landfall near the western tip of Galveston Island. This is generally noted as the most destructive hurricane path (Bedient, Blackburn, Anderson, Brody, & Colbert, 2015;Sebastian et al., 2014). Other paths could lead to worse conditions locally, but because this path heavily affects (the port of ) Houston, it is presumed leading for flood risk in the region.
The event frequencies-hydraulic conditions for different return periods-used in the evaluation model are based on  (2015), which includes extreme value analysis on FEMA datasets. Two additional sets of hydraulic boundary conditions are used, purely for verification. These are based on analyses conducted by Jackson State University (Ebersole et al., 2017) and simulate the conditions during Hurricane Ike and "Storm 36," which is a synthetically generated worst-case scenario. Both events have been modelled by Ebersole et al. (2017), providing flood maps for verification. For the simulation of Hurricane Ike, historical data on the hurricane track is used instead of the hurricane track mentioned above.

| Simulation model validation
The simulation model consists of hydraulic calculations and the damage calculations ( Figure 1). Therefore, the validation will compare these steps with the more advanced, federally accepted models ADCIRC and HAZUS. ADCIRC is a model that allows for interactions between waves and circulation, which allows the model to predict storm surge. The ADCIRC model has been used to accurately hindcast Hurricane Ike (Hope et al., 2013). HAZUS is a multi-hazard impact model, which can estimate the human, property, and financial impacts (FEMA, 2010). The HAZUS model benchmark was assessed previously (Atoba, Brody, Highfield, & Merrel, 2018). Because damage estimates of HAZUS also depend on water levels derived from ADCIRC, it is not possible to directly validate with damage estimates by MODOS. Therefore, a Baseline model is introduced (Figure 11), in which water levels from ADCIRC (by Ebersole et al., 2015) are combined with damage curves and property value data from MODOS. When compared with actual HAZUS results, this shows the accuracy of MODOS damage calculations.   Note. Cost estimates in reference literature can vary and often not divides between constant and variable costs. a Nature-based Solutions: can be separately added to any strategy. Connecting dredge spoils inside the bay and building oyster reefs can mitigate wave impact on the (north-) western part of the Galveston Bay shore. Wetlands can be placed in front of the Texas City Levee to mitigate wave impact. Sand nourishments can mitigate wave impact and erosion at the coast. The cost of sand nourishments is roughly based on reference projects and can only be implemented along the entire barrier island. b A levee on Galveston Island and Bolivar Peninsula combined with a storm surge barrier at Bolivar Roads is called the Ike Dike (Coastal Spine), a strategy proposed by Texas A&M University Galveston. c Combination of using dredge spoils and an in-bay storm surge barrier creates the Mid Bay barrier. In this case, the storm surge gate will not be place at the river inlet, but in the middle of the bay.
The flood water level and estimated damage was calculated for each watershed. This was done for the two storms-Hurricane Ike and Storm 36-and two flood protection configurations: current situation and Coastal Spine. Figure 12 shows the verification of water levels in the watersheds, which shows how well different parts are represented.
The validation (Table 5) shows a good fit for the hydraulic model (<20% error) compared to ADCIRC results, where the main differences are caused by local hydraulic processes and no systematic bias is expected. The difference in estimated damage can be the result of a difference in modelling or a difference in property inventory. HAZUS provides a better industrial and commercial inventory, which was not available for this study (FEMA, 2016). Residential losses are significantly overestimated when the flood levels are low. An issue that may explain these errors is the limited accounting of elevated structures. At these storm intensities, coastal flood protection can limit the flood level down to the point where an elevated house can make a significant difference (Jongman et al., 2012).Because elevated houses have not been taken into account, the damage for limited water levels is overestimated in comparison to actual flood events. However, this effect is equally present across all simulations and should be accounted for when formulating recommendations.
Because of the significant lack of property inventory and accurate damage curves available for the MODOS-model, most damage estimates are underestimating the impact of flooding. Nonetheless, because the estimates are consistent between simulations, it does allow for analyses in comparison to each other. Therefore, the risk reduction will not be quantified in dollars, but in the strategy's ability to reduce risk compared to the nullscenario (expressed in percentages). The damage calculations in MODOS have also been simplified to reduce computation time. Better accuracy can be achieved by decreasing the contour resolution, which is at 0.31 m (1 ft), or by increasing the amount of land use types taken into account.

| Evaluation of flood risk reduction strategies
In order to thoroughly explore the full design space, numerous strategies have to be considered. However, this also increases the computation time. Given the amount of measures, 500 different strategies were deemed sufficient. With this number, most combinations of measures will be represented with several different levee heights. Construction costs are determined using Table 3 and the environmental score is calculated using a formula listed in Appendix. Three output parameters were chosen: Risk reduction, construction cost and environmental score.
When the risk reduction is compared with the construction cost (Figure 13), there are clearly two different trends, based on the placement of a storm surge barrier at Bolivar Roads. This gate only results in higher risk reduction when the total construction costs are relatively high. If no coastal barrier is built on Galveston Island or Bolivar Peninsula, or when either is not strong enough to withstand extreme storms, significant damage will still occur around the Galveston Bay. The storm surge barrier at Bolivar Roads will only be effective when combined with large investments in other defences. Taking a closer look at the individual strategies, two strategies that seek to optimise risk reduction are highlighted in Figure 13. Both include strong coastal defences, but the main difference is revealed inside the bay. Strategy A includes a relatively low levee along the north-western edges of the bay, while strategy B includes a gate in the river, oyster reefs in the bay and wetlands to lower wave attack. The use of Nature-based Solutions and limiting damage to nature reserves clearly has a positive effect on the environmental score of strategy B. The small difference in risk reduction also shows how dominant the coastal defences are in risk reduction. Similar comparisons show that the impact of inner-bay defences greatly depends on the placement of coastal defences. Without coastal defences, the high water levels in the bay combined with wind set-up will quickly overpower inner-bay defences, greatly reducing their effectiveness.
Environmental impact is hardly correlated with the two other output parameters (see Figure 13). This is due to the Nature-based Solutions. These measures have a large effect on environmental impact, but are relatively cheap and have a small effect on risk reduction. They can therefore be implemented within every strategy, lifting the environmental score, while hardly affecting the cost and risk reduction of the total strategy.
The data can also be analysed using Random Forest Feature Scoring (Figure 14, based on Breiman (2001) and Pedregosa et al. (2011)). A higher number implies a higher impact of the measure (left) on the output parameter (down). The values are mainly useful in comparison to each other and enable the measures to be ranked. For example, risk reduction is primarily impacted by the choice for the coastal barrier on Galveston Island and Bolivar Peninsula. In agreement with what was presented above, Nature-based Solutions greatly impact the environmental score, but have almost no impact on the construction costs or risk reduction. Modelled is the situation where Hurricane Ike would hit the region, given that the "coastal spine"-strategy is implemented. The coastal spine is a locally proposed combination of intervention mostly focussing on coastal protection  Figure 15 is an example of a PRIM-analysis, which can find strategies that comply with specific demands. For example, we can search for a combination of low construction costs (<5 billion USD) with high risk reduction (>75%). The results show the importance of the coastal barrier on land. Moreover, it shows that the costs of storm surge barriers are too high in this case to find a system that includes the storm surge barrier and still meets the thresholds, even though the impact on risk reduction is second highest (see Figure 14).
Because of the computational time of the simulation and the possibility of comparing a large number of flood risk reduction strategies, it is possible to discover how the risk profile in the region changes for different design choices, stakeholder preferences and other developments. These insights can support the decision-making process early on, when impactful decisions are required based on limited information, and lead to a better optimised strategy for reducing flood risk in the Houston-Galveston Bay area.

| DISCUSSION
The model presented in this paper is built for use during initial exploration of flood risk reduction strategies and has been optimised to be fast and broadly applicable. Therefore, there are several limitations that should be considered. First, the verification showed a good fit for the hydraulic model to more advanced ADCIRC models (<20% error for water levels). Damage estimates had a higher uncertainty (<50% error for damages around the bay), mostly due to the higher uncertainty in damage curves and lack of commercial data. Although the errors in damage estimates are relatively high, the errors are relatively consistent across strategies and storms, which allows for the strategies to be compared on their flood risk reducing abilities. Some types of measures (e.g. Nature-based Solutions) still have a large uncertainty in their effect on the system (Vuik et al., 2016). The current assessment of environmental impact in the form of a score needs further development to be as reliable as the other parameters. Other measures like buyouts, changes in (insurance) policy, and raising houses have not been evaluated yet.
Second, because the region is divided into lines of defence, it is hard to accurately include measures that work on a much smaller scale. This problem can be addressed by including more lines of defence. However, adding more lines of defence will exponentially increase the amount of possible outcome scenarios and therefore the computation  The MODOS-model allows inclusion of several other features that have not been shown in this paper. The simulation only considered coastal storm surge, because of its destructive nature. However, most coastal regions like the Houston-Galveston Bay also suffer from inland flooding from rain or river runoff. In the future, the hydraulic module of MODOS can be extended to incorporate these types of events.
The connection with the EMA-workbench allows for the exploration of different future scenarios in terms of climate change or future land use. This can be used to find a robust design, able to provide sufficient safety in an uncertain future. Eventually, multiple-objective optimisation techniques can also be used to identify the Pareto approximate set of optimal strategies. Here, the modeller can choose to include additional output variables, with practically no extra computational load. Finally, the model is specifically built to be applicable in different regions in the world, which can be established through additional case studies. Especially regions with large variation in land use and topography, where many flood risk reduction measures are possible, the model can provide valuable insights to the design discussion.

| CONCLUSION
This paper presents the MODOS-model and shows its application and potential for the Houston-Galveston Bay area. The MODOS-model can (a) simulate the effects of a flood risk reduction strategy on a complex, flood-prone region and (b) evaluate the results of different strategies, consisting of multiple interventions, based on economic and non-economic performance indicators. With the MODOS model, it is possible to probabilistically assess the impact of both structural and non-structural measures on flood risk reduction. These interventions can range from levees and storm surge barriers to Nature-based Solutions like wetlands and oyster reefs. By using simplified hydraulics, the model is kept computationally workable, which allows for many different strategies to be compared. The model setup is kept generic, which enables the model to be transferred to other locations relatively easily.
A simplified case study based on the Houston-Galveston Bay Area in Texas was used to demonstrate the model. First runs showed that the flood risk reducing capabilities of the chosen strategy mainly depend on the choice of coastal land barrier, because most strategies that seek to maximise both risk reduction and cost efficiency included either a coastal levee or floodwall. Also, Nature-based Solutions are shown to have a small impact on the construction costs and the risk reduction of the total flood risk reduction system. However, they are expected to have a clear, but local positive effect on environmental impact of the constructed system.
The use of the EMA-workbench in the evaluation model enables the modeller to include analysis of uncertainties in design, and future scenarios like future land use scenarios or climate scenarios. This is not included in this research. The MODOS model is designed to support decision making during design discussions. Especially during the conceptual design phase, when design choices are impactful and information is scarce, it can provide valuable insights and save time and resources in identifying promising strategies for reducing risks for vulnerable flood-prone areas.
NOTES 1 Portion of structure/content damaged at a given inundation depends on the land use type. This model uses seven land use types, which each have their own damage curve. 2 https://comptroller.texas.gov/economy/economic-data/ports/houston.php Table A1 provides a description of the code used for the hydraulic calculations. Listed below are the equations used in that process. In Table A1, the numbers between brackets refer to the applied equation.

TABLE A1
Hydraulic simulation used for a line of defence in the MODOS-model. Numbers between brackets refer to equations Build inside water level series per outcome scenario (Equation (2)) Find maximum water levels in the protected area per outcome scenario  (5) and (6)) Calculate Expected damages (Equation (4)) Coding used for calculating the maximum inland basin water levels in the Simulation Model for every protected area Calculate maximum basin water level (Equation (3)) Hydraulic conditions for next line of defence for each scenario: Calculate wind setup at barrier location (Equation (A6)) Calculate maximum water level at barrier location (Equation (3)) Scale time series of average bay water level if (water level higher than land elevation): wave height on land = minimum(incoming wave height, 0.5 × local water depth ) wave height on water = minimum(incoming wave height, 0.5 local water depth)