Flood hazard reduction from automatically applied landscaping measures in RiverScape, a Python package coupled to a two-dimensional flow model

Abstract River managers of alluvial rivers often need to reconcile conflicting objectives, but stakeholder processes are prone to subjectivity, time consuming and therefore limited in scope. Here we present RiverScape, a modeling tool for numerical creation, positioning and implementation of seven common flood hazard reduction measures at any intensity in a 2D hydrodynamic model for a river with embanked floodplains. It evaluates the measures for (1) hydrodynamic effects with the 2D flow model Delft3D Flexible Mesh, and (2) the required landscaping work expressed as the displaced volume of material. The most effective flood hazard reduction in terms of transported material is vegetation roughness smoothing, followed by main embankment raising, groyne lowering, minor embankment lowering, side channel construction, floodplain lowering and relocating the main embankment. Implementation of this tool may speed up decision making considerably. Applications elsewhere could weigh in adverse downstream effects, degradation of the ecology and overly expensive choices.


Introduction
Flood risk reduction ranked high on the political agenda over the last two decades, which is warranted given the high and increasing societal cost of flooding, the anticipated ongoing climate change, and economic developments in fluvial and deltaic areas (Hirabayashi et al., 2013). Here, flood risk is defined as the inundation probability times the inundation effect. The European Flood Directive (European Commission, 2007) states that it is feasible and desirable to reduce the risk of adverse consequences associated with floods, and obliges member states to create flood hazard and risk maps, and a flood risk management plan for the implementation. Flood risk management can be summarized by (1) strategy, i.e. protection against floods, living with floods, and retreat to flood-safe areas, and (2) timing of the action relative to the flood event, i.e. pre-flood preparedness, operational flood management and post-flood response (Kundzewicz and Takeuchi, 1999). Consequently, river managers are confronted with large challenges in the planning of measures in and around floodplains of embanked alluvial rivers, not only due to the number of stakeholders involved, but also due to the long lasting effect on the landscape, economic development and riparian ecosystems (Pinter, 2005).
Flood hazard management at the river basin scale consists of storing water in the headwater of the basin, retaining water instream in the middle parts and discharging the water in the downstream reaches (Hooijer et al., 2004). This is because the propagation of a flood wave, or flood wave celerity, increases with the flow velocity of the water and with the fraction of the discharge conveyed by the main channel (Jansen et al., 1979). For example, the narrowing of the floodplains by embankments and decreasing the flow resistance of the floodplain vegetation increases the flood wave celerity, which adversely affects the flood hazard downstream (Clilverd et al., 2016). Here we present a flexible tool for quantifying effects and effectiveness of common measures to lower the flood risk with the aim to support stakeholder discussions with evidence-based facts and figures. We develop and apply the tool to a specific case of a lowland deltaic floodplain at the downstream end of the river Rhine, which is a medium-sized river draining part of North-West Europe.
Typical measures at the scale of a floodplain section ( Fig. 1) have in common that they increase the water storage, and the conveyance capacity during floods. Two types of measures are considered here to lower the flood hazard, more specifically, the probability of flooding the embanked areas. The first type lowers the flood stage during peak discharges (measure type 1 to 6, Fig. 1) by creating more space for the river within the embankments. The second type comprises raising the main embankment, which enables higher water levels. The flood hazard reductions of these measures have been reported previously (Baptist et al., 2004;Remo et al., 2012), and are routinely evaluated in operational river management. The typical workflow comprises a geodatabase with spatial information that is converted to input data for a hydrodynamic model. Experts, together with stakeholders, choose what measure will be implemented, and manual adjustments are made to the geodatabase and the derived hydrodynamic model. Expert judgment drives this process, which is limited by the amount of manual work required to update the hydrodynamic model with a realistic bathymetry and land cover at the spatial extent of the measure. These processes can take years for simple measures, and more than a decade for complicated projects due to the complex and iterative nature of joint decision making. Decision support systems (DSS) for these long term planning projects in the preparedness phase are scarce, contrary to DSSs for operational flood management.
The options for flood hazard management for the lower reaches of the River Rhine in the Netherlands (Silva et al., 2004) were modelled for individual measures, and the water level lowering at the river axis were made available in a graphical user interface (WLjDelft-Hydraulics, 2008). Interactive planning of some measures was possible using geospatial software (Van der Werff ten Bosch, 2009). Application at the river-reach scale with realistic measures, however, is tedious and impractical, showing a need for automated procedures to generate these measures in larger areas. Measures can be applied with different gradations and spatial extents, to which we will refer to as 'intensities of application'. The units of this intensity vary, e.g. small and large side channels, or relocation of embankments over short or large distances. Nonetheless, each measure lowers the flood hazard and their implementation requires material displacement. Our main objectives were to (1) develop a tool to automatically position and parameterize seven flood hazard reduction measures and (2) evaluate these measures on hydrodynamic effects plus the required volume of displaced material. These aims are limited to the physical domain; evaluation on costs was outside the scope of this study, even though it is closely related to transported material. We developed the RiverScape package in Python and applied it to the main distributary of the River Rhine. The results are followed by discussion of the applicability to other alluvial rivers and future perspectives to incorporate values other than material displacement.

Materials and methods
We developed RiverScape, a Python package, which uses map algebra functions from PCRaster (Schmitz et al., 2013). RiverScape can position and parameterize landscaping measures and update the input data for the two-dimensional (2D) flow model Delft3D Flexible Mesh (DFM), which is also open source. It requires input on hydrodynamic boundary conditions, a geodatabase with layers of river attributes, and settings to determine the intensity of application for each measure (Fig. 2). Once the measures are known, we updated the 2D flow model's input in order to determine the flood hazard reduction and the flow velocities. Here, we present the methods implemented.

Study area and available data
The case study area is located in the Rhine delta, which consists of three distributaries: the Rivers Waal, Nederrijn and IJssel. We selected the River Waal, which is the main distributary of the River Rhine in the Netherlands (Fig. 3). The three main concerns here are flood risk in view of global change, navigability and ecosystem functioning. The study area spans an 94-km-long river reach with an average water surface gradient of 0.10 m/km. The total area of the embanked floodplains amounts to 132 km 2 . The main channel is around 250 m wide and fixed by groynes. The cross-sectional width between the primary embankments varies between 0.5 and 2.6 km. Meadows dominate the land cover, but recent nature rehabilitation programs led to increased areas with herbaceous vegetation, shrubs and forest. The design discharge for the River Waal is now set to 10,165 m 3 s -1 , which has an average return period of 1250 years. Such a discharge is expected to give a 3.99 m water level above ordnance datum (þOD) at the downstream end of the study area. The main channel functions as the primary shipping route between the port of Rotterdam and major industrial areas in Germany. The main channel position is fixed in place by groynes, which were partly lowered during the 'Room for the River' project (Van Stokkom et al., 2005). In the future, the design discharge will be combined with a risk-based approach that takes the potential damage and casualties within the protected areas into account (Van Alphen, 2016).
The spatial data describing the major rivers in The Netherlands are stored in an ArcGIS file geodatabase according to the Baseline data protocol, version 5 (Scholten and Stout, 2013). This protocol, specific for the Netherlands, describes the layers in the geodatabase and specifies the required attributes for each of the layers in terms of names, and properties. Baseline schematizations include layers with (1) land cover as a polygon layer of ecotopes (landscape-  (Middelkoop and Van Haselen, 1999)). ecological units), (2) hydrodynamic roughness as point, line, and polygon layers, (3) minor embankments, groynes, and main embankments as 3D lines consisting of routes and events, and (4) river geometry describing the extent of main channel, groyne fields, floodplains as polygons (Fig. 4A). We adhered to the Baseline schematization in order to allow comparison between our results and existing projects, but in principle any other method of input of the aforementioned data can be used in combination with River-Scape. All spatial river attributes are represented as vector layers, except bathymetry, which is represented as a Triangular Irregular Network (TIN) for the main channel and the floodplain. The TIN represents the ground level and does not include the groynes and minor embankments. We used the 'rijn-beno14_5-v2' schematization of the Rhine branches, which describes the layout after the finalization of the measures of the 'Room for the River' program in December 2015. In the areas protected from flooding by the embankments additional data sources were required as Baseline only covers the embanked floodplains and the main channel. The national LiDAR-based Digital Terrain Model (DTM) provided terrain elevation data. This gridded DTM has a 0.5 m resolution, and a vertical error less than 5 cm (bias and random error) in open terrain (Van der Zon, 2013). Building locations were derived from the national database of addresses and buildings (BAG, 2016).
The Ministry of Infrastructure and Environment provided the computational mesh of the flow model (Baart and Scholten, unpublished data). It consisted of 24 small quadrilateral cells across the main channel and groyne field, sized around 40 by 20 m (Fig. 5). These are connected by triangular cells to large quadrilateral cells in the floodplains, sized around 80 by 80 m. No mesh refinement was implemented around the individual groynes to limit the computation time. We extended this mesh with triangular cells for the areas protected from flooding by embankments. Discharge and water level time series between 1989 and 2014 were obtained from the gauging station at Tiel (Fig. 1) (www.live.waterbase.nl).

Hydrodynamic modeling
RiverScape was coupled to a 2D hydrodynamic model, and required a calibrated model as a starting point. In this study, we used DFM, the open source hydrodynamic model that is being developed and maintained by Deltares (2016). The computational core of DFM solves the shallow water equations based on the finitevolume methods on an unstructured grid (Kernkamp et al., 2011). DFM's computational mesh and output are stored in netcdf files that follow the UGRID conventions for specifying the topology of unstructured and flexible (triangular, quadrilateral, etc.) grids (UGRID, 2016). The computational mesh of the study area consisted of 120,000 cells, of which 71,000 were active with the current location of the major embankments.
The spatial DFM input consists of five components in either netcdf or ASCII format. Firstly, the ground level of the bathymetry is derived from the TIN in the Baseline geodatabase, which is converted to netcdf format using a predefined computational mesh. Secondly, linear terrain features, such as groynes, minor embankments, and steep terrain jumps, are defined as ASCII-formatted line elements. These linear features cause additional energy loss when submerged. They are excluded from the bathymetry to limit the number of computational cells and the associated long computation time. These linear elements are called fixed weirs and contain the coordinates, the height difference on the left and right side, and the width and slope of the linear feature. Thirdly, the hydrodynamic roughness is based on trachytopes: spatially-distributed, and stagedependent roughness values. Trachytopes are based on points for single trees, on lines for hedge rows, and on polygons for land cover derived from the ecotope map. Trachytopes in the main channel are adjusted in the model calibration. For each flow cell in the computational mesh the fractions of each trachytope is given, e.g. 0.7 for trachytope X and 0.3 for trachytope Y. Ch ezy roughness is computed at runtime within DFM using the water depth dependent roughness equations developed by Klopstra et al. (1997). Fourthly, obstacles that can not be submerged, such as bridge pillars or houses, are implemented as so-called thin dams. These represent infinitely high obstacles to flow that may consist of lines, or polygons. Finally, dry areas are only represented as polygons that render the contained computational cells inactive, whereas for thin dams the cell remains active and water can flow around the obstacle.
We generated the DFM spatial input from the Baseline geodatabase using the Baseline plugin for ArcGIS (Scholten and Stout, 2014). Further, we extended the trachytope definitions with missing codes, defined the discharge time series on the upstream boundary, and compiled a rating curve at the downstream boundary. This completed the DFM model setup.

Compilation of river attributes
RiverScape works on a gridded representation of the river to increase computational speed. The basic data consists of a terrain model, land use, measured water levels, a rating curve, and a functioning hydrodynamic model. This makes application in areas that are more data scarce than the Netherlands feasible. For the study area, vector-based data were available in Baseline, which were rasterized to a 25 m raster resolution, to ensure that the cell area of the rasters (625 m 2 ) was smaller than the cell area of the computational mesh of the 2D flow model (>800 m 2 ). Some subgrid roughness information is lost in this way, because a single flow cell may contain multiple trachytopes. In DFM, this information is maintained as each flow cell may store fractional trachytope areas. Thirteen relevant Baseline layers (Table 1) were rasterized to a common map extent and resolution in PCRaster format using the Geospatial Data Abstraction Library (GDAL) (Warmerdam, 2008). The second set of attributes, needed for the automatic positioning of measures, gave additional information on the river geometry, such as channel curvature, curve direction, and separate floodplains sections (Table 1). For example, a good location for embankment relocation is an area with a sharp right turn in the river axis, a narrow floodplain on river right, and a low total value in real estate. Floodplain width calculation was challenging as it posed a one-tomany problem: many points on the main embankment in the outer bend could be connected to a limited set of points on the channel bank, and in the inner bend many channel points could be connected to a single cell on the main embankment. Here we pragmatically calculated the distance from each embankment cell to the channel bank on a line perpendicular to the channel center line (Fig. 4E). Crossing lines were redirected towards the nearest embankment, while maintaining the highest width value. The radius and turning direction of the river were derived from fitting a circle to the river axis at each axis cell. The location of the center point of the fitted circle in river left or right determines the turning direction.
The third set of river attributes represented hydrodynamic characteristics derived from a reference run with the 2D flow model ( Table 1). The discharge was increased in a stepwise manner between low flow and design discharge, and each step was maintained for 4 days to create a stationary flow. We used discharges of 698, 1481, 1713, 2157, 2935, 4966 m 3 s -1 , which are exceeded 363, 150, 100, 50, 20, and 2 days per year, respectively based on the time series of the Tiel gauging station. The discharge values were derived from percentiles derived from the exceedance percentage in days. The choice for these exceedance values was based on their ecological significance as implemented in the classification method for the Dutch ecotope map of the large water bodies (Klijn and de Haes, 1994;Van der Molen et al., 2000). This map is periodically made with a standardized method for management purposes. For example, in areas that are inundated less than two days per year the vegetation is considered unaffected by inundation. The low exceedance values are ecologically relevant for vegetated floodplains, and the high exceedance values are related to the water bodies and side channels. In addition, we ran the model with the design discharge for the River Waal of 10,165 m 3 s -1 , which provided data on water depth, flow velocity and hydrodynamic roughness ( Fig. 4BCD), amongst others. Cell IDs were required for fast nearest neighbor interpolation of the data stored in an unstructured mesh to regular rasters.

Positioning and parameterization of seven flood hazard reduction measures
We developed automated procedures to determine the flood hazard reduction potential of seven landscaping measures by adjusting the input of the 2D flow model (Table 2). Six flood stage lowering measures in the floodplain and groyne field were determined plus embankment raising. Measure positioning was required to determine suitable locations for roughness smoothing, side channel construction, floodplain lowering, and embankment relocation. Table 3 gives a summary of these methods and settings. Groyne lowering, minor embankment lowering, and main embankment raising do not need positioning as their location is predefined in Baseline. Each of the measures was applied with six intensities of application. We omitted deepening of the main channel as a flood hazard adaptation option, because the River Rhine is already deepening due to reduced sediment input from the basin and the narrowing of the main channel with groynes (Frings et al., 2014a). Bed degradation between Tiel and the downstream model boundary was approximately 5 mm/y based on (Frings et al. (2009) , Fig. 5). Bed erosion negatively affects shipping at nonerodable outcrops, infrastructure and ecology due to lower water levels (Frings et al., 2014b) and exposes entrenched telecom cables and pipelines. Technically, the implementation would be similar to floodplain lowering.
Each of the measures requires the transport of material (soil, stones, vegetation, or high quality clay for the main embankments) when implemented in the field to create more space for the river. The volume of moved material does not necessarily equal the added volume available for water. For floodplain smoothing, the volume of moved material is larger than the water volume as the emergent vegetation also needs to be removed. Contrarily, in the case of raising the main embankments, the volume of water is much larger than the volume of soil required for raising. To compare the flood stage reduction effect of the different measures, we calculated material volume that needs to be moved, and water volume created for each of the measures. Separate volumes were calculated for (1) vegetation based on stem densities and stem diameters per roughness class (Van Velzen et al., 2003) and mean volume per vegetation class (Schelhaas et al., 2014), (2) material in groynes and minor embankments derived from the 3D attributes, and (3) soil transport, based on the differences between the current and new DTM. We limited the evaluation to the physical domain, i.e. the conveyance capacity using the lowered flood levels and the storage capacity by calculating the increased water volume. The logical extension of the evaluation on transported material would be the cost of the measures, but this was outside of the scope for this study. Flood wave celerity was not considered, since the study area is close to the river mouth and flood waves are long relative to the study reach length.

Side channels
Side channels were created in a two-step approach, which comprised positioning of the channel center line, followed by the parameterization of the cross section shape and hydrodynamic roughness. Firstly, we positioned side channels only in wide floodplain sections (Fig. 4E) without side channels currently present. Over each section, the start and end point of a side channel were positioned on the river axis alongside the upstream and downstream end of the section. The centerline of the side channel was determined by the path of the least resistance between start and end point. A high resistance value was assigned to the main channel and groyne field to force the centerline into the floodplain. A low resistance was given to existing floodplain backwaters, and a resistance based on distance to the main channel and main embankment was given to the remaining areas (Fig. 4F). The upstream end of the side channel was disconnected from the main channel to prevent large morphological changes in the main channel. Secondly, we parameterized the side channels with a trapezoidal shape of which width, depth, and cross-sectional slope can be set with user-specified values. The new side channel is only defined where its depth is below the current bathymetry, leaving existing lakes largely untouched. For the largest side channel, we set the depth as an offset of 2.5 m below the water level at the river axis exceeded 363 days per year, the width was set to 75 m, and the bank slope to 1:3 (Fig. 4G). In this study, we implemented a series of six side channels in each of the suitable floodplain sections. The six intensities of application were defined by the depth and width values scaled to 10, 20, 40, 60, 80, and 100 percent of the value used for the largest side channel.

Vegetation roughness smoothing
A low vegetation roughness increases the conveyance capacity of the floodplain area lowering the overall water levels (Baptist et al., 2007;Aberle and J€ arvel€ a, 2013). There is no standard procedure for choosing where to lower the roughness, and in practice it is based on the judgement of the river manager and contested by nature developers. We developed a method that optimizes roughness smoothing by selecting the areas where lowering the floodplain roughness is most effective in terms of water level lowering. This is the case where a high specific discharge (q, (m 2 s À1 )) coincides with a high vegetation roughness, expressed as the Nikuradse equivalent roughness length (k, (m); Fig. 4D). For example, a dense forest at the outflow point of a floodplain section would be a big obstruction to flow. We calculated a, the product of two fields q and k, and determined its cumulative frequency distribution (cfd a ). The score of cfd a at a specific percentile of the distribution was used as a threshold for positioning the roughness smoothing. Areas where a exceeded the percentile score were selected for roughness smoothing. The percentile was calculated as 100 minus a userspecified percentage of the terrestrial floodplain area (P sm ). For example, floodplain smoothing over 10% of the floodplain area is positioned where the score at the 90 th percentile of cfd a is exceeded. The vegetation type at the selected areas was changed into production meadow, the vegetation with the lowest roughness. Intensities of application were set to floodplain smoothing over 1, 5, 10, 25, 50, and 99% of the terrestrial floodplain area (Fig. 4H). The increasing increment was chosen, because of the decreasing effectiveness of this measure as the current land cover also includes production meadows.

Floodplain lowering
Floodplain lowering was positioned using a similar method as for roughness lowering. It is most effective where a high flow velocity coincides with a low water depth under peak discharge, for example in case of flow over a natural levee deposit at the upstream end of a floodplain section. We calculated the product of two fields, denoted as b, (1) the water depth and (2) flow velocity field subtracted from the maximum flow velocity at design discharge. The inverse of the flow velocity was chosen to prevent equifinality in the selection. Floodplain lowering was positioned where b exceeded the score at percentile of cfd b , where the percentile equals 100 minus a user-specified percentage. The new terrain elevation was set to the height corresponding to the 50 days per year (d/y) flood duration. We chose the roughness code for production meadows as the new roughness. This smooth land cover adds to the flood level lowering, but RiverScape is flexible in assigning new codes. Like floodplain smoothing, we increased the intensity of application by applying floodplain lowering over 1, 5, 10, 25, 50, and 99 percent of the terrestrial parts of the floodplain (Fig. 4I). Within the schematization of the flow model, bathymetry, roughness, and fixed weirs were updated (Table 1).

Embankment relocation
Embankment relocation can be implemented as a buffer around the current main embankment, but it is more efficient when the embankment is straightened locally, especially with a tortuous embankment shape in top view. Therefore, we relocated the embankment using an alpha shape derived from the embanked area. An alpha shape (Edelsbrunner and Muecke, 1994), also known as the concave hull, is based on a Delauney triangulation of a point set where long edges are removed based on the alpha value. The lower the alpha value, the more it follows the current embankments, while an infinitely high alpha value gives the convex hull. We increased the intensity by using alpha values of 500, 1000, 2000, 3000, 5000, and 7000 m (Fig. 4J). We took current built-up areas into account by creating ring dikes around areas with high building costs. Relocation was implemented by adjusting the dry areas. The vegetation type of the new floodplain area was set to production meadow.

Lowering of groynes and minor embankments
Minor embankments have been constructed to prevent the inundation of agricultural fields during minor floods and their crest level varies. Likewise, groynes have varying crest levels as some have been lowered within the 'Room for the River' project to lower flood water levels. Lowering of groynes and minor embankments was implemented using their current locations, as stored in the flow model's input. Groynes and minor embankments are both stored as 'fixed weirs' in DFM. Fixed weirs are three-dimensional lines that describe location and crest height (xyz) for each vertex along the line. In addition to the crest height, each vertex contains Alpha shapes of 500, 1000, 2000, 3000, 5000, 7000 m information on the cross-sectional shape of the fixed weir: the terrain height on the left and right of the crest (the toe of the minor embankment, or groyne), the cross-sectional slope, and the crest width. Lowering can be applied as a percentage of the current height, an absolute change in height, or by using an external height level such as a water level that is exceeded a fixed number of days per year. Due to the current differences in crest height, we chose to homogenize the differences by applying an external height. The new height consisted of the minimum of the current crest height and the external level, but the crest height should not be lower than the terrain height left or right of the crest. The intensity of application was increased by lowering to flood durations of 50, 100, 150, 200, 250, and 366 d/y for groynes and to 2, 20, 50, 100, 150, and 366 d/y for minor embankments. The 366 d/y flood duration involved the complete removal of the groynes and minor embankments from the fixed weir input of the 2D flow model.

Results
We aimed at the evaluation of flood stage reduction effects of seven landscaping measures with increasing intensity and its relation to displaced material. Calculation time for rasterizing the input data, calculating the derived information, positioning and parameterization of measures, and updating the flow model output required 0.5, 1, 2, and 1.75 h, respectively, on a i7-6700 3.4 GHz processor using a single thread. The calculation time for updating the flow model input includes conversion of the updated data to GIS-compatible raster and vector layers. The initial hydrodynamic calculation with the stepwise increase in discharge and the ensemble of realizations both took 24 h. We first describe the measures that resulted from the automated methods and then describe their effects on flood water level and channel flow velocity changes. Finally the results are expressed as a function of the required material displacement volume for the measures and the additional space for flood water due to the measures.

Positioning and parameterization of measures
The spatial layout of the six flood stage lowering measures is given in Fig. 6 for the downstream section, and in Fig. S1 in the supporting information on an A3-sized figure for the whole study area. Roughness lowering locations ( Fig. 6B and Fig. S1B) coincide with forested areas when applied to 1% and 5% of the terrestrial part of the floodplain area. Examples include the forest on river left at river kilometer (rkm) 872.5, on river right at rkm 878.5, and on river left at rkm 889. These forest patches are often located in areas with low specific discharge. At higher intensities of applications (50, or 99%), also herbaceous vegetation, single trees, and hedgerows are removed and converted to meadows.
Groyne lowering was applied to the 797 individual groynes along the main channel. The current height of the groynes is not the same over the river reach. Between rkm 876 and 886, 914 and 922, and 952 and 960 the groynes are around a meter above the water level exceeded 100 d/y, whereas in the remaining stretches the groynes are slightly below this level ( Fig. 6C and Fig. S1C). The section between rkm 911 and 928 does not contain groynes in the inner curves as they were converted to longitudinal training dams, which are treated as minor embankments in the hydrodynamic input (Fig. S1D). Not all groynes were equally affected by groyne lowering, due to the spatial differentiation of the current groyne height. When groynes are lowered to the water level exceeded 250 days per year, the median height difference at the endpoint of the groynes between the crest and the highest groyne toe is still 2.93 m, which is marginally higher than the minimal guaranteed depth of 2.8 m (Van Vuren et al., 2015).
All of the 223 km of minor embankments were lowered to increasingly long flood durations. The maximum lowering of the crest (Fig. 6D and Fig. S1D) was limited by the maximum toe height, which represents the ground level of the floodplain. Heights of the minor embankments above ground level vary strongly and they showed a 1.25 m interquartile range (0.31e1.56 m). The highest minor embankments are found upstream from rkm 883, especially on the river right.
New side channels were planned in 16 out of the 29 wide floodplain sections (Fig. S1E). Positioning of side channels was comparatively demanding computationally (20 min computation time) as all floodplain sections were addressed sequentially. The centerlines follow the midpoints between the main embankments and the main channel in sections without water bodies present, e.g. at rkm 895.5 on river left. In curved sections, the center line is drawn towards the inner part of the curve within the floodplain section (i.e. rkm 875 on river right), and when water is present, the centerline is drawn towards existing water (i.e. rkm 880 on river right). We positioned one new side channel in the floodplain section on river left around rkm 930 (Fig. 6E). In reality, a side channel was created here as well (inset Fig. 9) at a similar position. In contrast to our modeling choices, the side channel was connected to the main channel at the upstream end as well.
Floodplain lowering at 1 and 5% of the terrestrial floodplain area mainly affected artificially raised industrial areas ( Fig. 6F and Fig. S1F), such as the shipyard at rkm 897.5 on river left. At 10e25% natural levee deposits are removed, whereas at 50e99% also the low lying sections are lowered to the water level that is exceeded 50 d/y. It should be noted that the difference between floodplain lowering and groyne lowering is defined in Baseline, which states that the maximum width of minor embankments is 10 m. Wider areas are contained in the bathymetry. This is visible in the lowering pattern as elongated lines, e.g. at rkm 922 on river right (Fig. 6F).
Embankment relocation was carried out with six increasing alpha shapes values, while existing real estate was taken into account ( Fig. 6H and Fig. S1H). The larger alpha shapes almost doubled the surface area of the embanked floodplains in the tortuous upstream part and around rkm 923 to 934, where existing villages become islands in the floodplain. The straight sections led to elongated new floodplain areas, such as on both sides of the river between rkm 888 and 898.

Hydrodynamic effects of measures with increasing application intensities
We derived water levels at the river axis and depth-averaged flow velocities from the 2D flow model. The simulation period was set to three days with a 10165 ms-1 discharge, which ensured that the flow became fully stationary. Wall clock time of a single simulation on a single core was around 4 h. Differences in water level between the reference run and runs with measures gave insight in the flood stage reduction of each measure and each intensity (Fig. 7). The downstream boundary condition, a fixed water level creating a backwater effect, led to zero change at rkm 961. Changes in flow velocities from the measures (e.g. Fig. 8 for side channels) indicate possible adverse morphological effects in the main channel that could hamper navigation. Flow velocity differences at 25 and 75% of the main channel width (Fig. 9) summarize the potential morphological effects.

Flood level reductions
Roughness smoothing increasingly lowered the water levels with larger extents of application (Fig. 7A). The effectiveness reduces at larger percentages, with the water level reduction between 1 and 5% almost equal to the reduction from 25 to 50%. The 99% intensity showed a negligible effect compared to 50%. The lowering was equally distributed over the area, with a maximum lowering of 0.2 m.
Groyne lowering showed three distinct steps in the lowered profiles at rkm 883, 920, and 958 for all lowering intensities except for the lowering '366 d/y' (Fig. 7B). The three steps coincide with the sections where the groyne height exceeds the 100 d/y exceedance level. The maximum reduction is 0.1 m. The '366 d/y' groyne lowering involves the complete removal of the groynes from the flow model. The resulting 0.25 m reduction serves as a reference for a more natural river with active meandering, which is not feasible for the Waal.
Minor embankment lowering reduced water levels when lowered to the 2, or 20 d/y flood duration, with a 0.12 m maximum (Fig. 7C). Lowering to flood durations of 50, 100, and 150 days did not lead to additional water level reduction, which indicates that the lowest ground levels around the minor embankment have an inundation duration of approximately 20 d/y. Similar to groyne lowering, we completely removed the minor embankments from the flow model input, which was labeled as '366 d/y' for consistency. This led to an additional 0.1 m reduction in predicted flood levels. At 150 d/y all minor embankments are effectively reduced to a terrain jump, because the 150 d/y level is lower than the terrain height. The additional 0.1 m reduction is due to the neglect of the energy loss from the terrain jumps that are in the current bathymetry, but that are lost in the relatively coarse bathymetry model. The real world implementation would be to alter the terrain in such a way that the downstream slope of the jumps is less than 1 to 7 to avoid flow separation.
Side channel construction led to flood water level reductions that increased with increasing cross sectional area (Fig. 7D, area indicated as percentage). Each increase in intensity did not lead to equal steps in the lowering. For example, the side channel around rkm 933 on river right (Fig. 6E) showed a 0.06 m lowering from 40 to 60% intensity, which is larger than for the other steps in intensity. This nonlinearity resulted from the upstream end of the disconnected Fig. 6. River attributes, measure locations and effects for a 30 km subset of the study area. See Fig. S1 for the full extent. A) River geometry with the main land cover classes and river kilometers at the river axis. B) Roughness smoothing locations at six percentages of application. Higher percentages include the location of lower roughness. C) Groyne height above the water level at the river axis exceeded 100 days per year. D) Minor embankment locations and their maximum lowering potential, i.e. the height above the highest one-sided ground level. E) Side channel location and bathymetry at the maximum intensity. The trapezoidal cross section shape had a 75 m width and a depth 2.5 m below the 363 days per year flood duration level at the river axis. F) Locations for floodplain lowering. Higher percentages include the areas of the lower percentages. G) Embankment relocation at six different alpha shapes. Dense built-up areas are ignored. side channel. At 40% intensity it does not affect the minor embankment here, whereas at 60% the extent is larger and the minor embankment was removed increasing the discharge capacity of the floodplain. Maximum lowering was 0.38 m; small differences in lowering were present between the 10 and 20% intensities.
Floodplain lowering and embankment relocation resulted in flood level reductions that differed an order of magnitude with the other measures ( Fig. 7E and F). Maximum reductions are 1.6 and 2.1 m for lowering and relocation, respectively. Floodplain lowering increased the flood level reduction in upstream direction, whereas relocation led to strong local reductions with their own backwater effects.

Flow velocity differences in the main channel
While the main focus of our work is flood risk, here we also studied changes in flow velocity in the channel. This is important because of the morphodynamic response: a spatial gradient in flow velocity leads to a gradient in sediment transport, and the latter gradient causes erosion and sedimentation in the shipping fairway, which would require dredging. The flow velocity along the channel shows patterns as a result of channel convergence and divergence and of exchange with the floodplain. Fig. 9AeC shows minor changes in flow velocity, but still large gradients in width averaged velocities. Here we averaged flow velocities over the main channel width and compared the velocity against the reference scenario. We assume that the flow velocity in the reference run does not cause erosion and sedimentation that require dredging for fairway maintenance, but this is not strictly true because maintenance dredging is conducted frequently.
The key result is that construction of side channels, floodplain lowering and relocation of embankments have significant effects (Fig. 9DEF). Zones of reduced flow velocity appear adjacent to the modified floodplain sections. On the other hand, floodplain smoothing, groyne lowering and removal of minor embankments in the floodplain hardly cause changes in flow velocity. The local velocity reduction is up to 0.5 ms À1 for side channels and floodplain lowering upstream from rkm 880, which is significant given that typical flow velocity in the channel is 1.75 ms À1 so considerable sedimentation is expected to result. This trend is even much stronger for embankment relocation, which implies dramatic morphological change in the river channel.
The relatively modest reduction of sediment transport in the side channel and floodplain lowering measures should also be evaluated against the sediment balance of the river Waal. Over the past decades, possibly more than a century, the river bed eroded in response to the installation of the groynes, due to dredging and due to reduced upstream sediment supply (Frings et al., 2009). Our modelled reduction of sediment transport capacity in the channel counteracts this trend, meaning that floodplain modification potentially has a positive effect. These results point at the need to adapt measures along the river such that changes in the gradients of sediment transport in the channel are minimized. This is not the same as the present strategy to minimize changes in sediment transport magnitude.

Water volume increase and flood hazard reduction from seven landscaping measures
Each of the seven measures involves the transport of one or more types of materials: (1) vegetation from roughness lowering, (2) stones and soil from the adjustments of groynes, minor embankments and main embankments, and (3) soil from the ground level (Fig. 10). The material volume varied strongly per measure type and intensity. The vegetation volume for roughness smoothing was 15% larger than the water volume (Fig. 10A) with a maximum vegetation volume of 1.3 10 5 m 3 . For the lowering of groynes and minor embankments material volume equals the water volume ( Fig. 10B and C). The material volume for main embankment raising was exceeded by a nine times larger water volume (Fig. 10D). Side channel construction and floodplain lowering involved all three material types as the vegetation and minor embankments are removed as well at the measure extents. Vegetation volume is at least an order of magnitude smaller than the volume for groynes and embankments, which is again an order of magnitude smaller than soil from the ground level. Interestingly, the vegetation volume triples when floodplain lowering is doubled from 50 to 99% due to the low lying vegetated areas. For Fig. 7. Flood level lowering for six measures and six intensities of application relative to the reference model run. A) roughness smoothing, B) groyne lowering, C) minor embankment lowering, D) side channel construction, E) floodplain lowering, F) embankment relocation. Note that the y-axes of panels A to C are scaled to À0.3 m. Panel D, E and F are scaled to À0.4, À1.4, and À2.5 m, respectively, which indicates their much larger effect on flood water levels. embankment relocation, the material volume is larger than the water volume at an alpha shape of 500 m. This confirms that relocation over small areas is inefficient, but the material volume barely increases with larger alpha shapes. Relocation with the alpha shape at 7000 m provided the largest water volume, 2.8 10 8 m 3 , which is 23% of the total water volume in the study area during current design water levels.
The relation between the flood hazard reduction and the volumetric changes per measure (Fig. 11) provided a concise overview of the effectiveness of the different river management options. Solid and dashed lines indicate the material, and water volume, respectively. A small target of 0.05 m flood stage reduction could be achieved by all different measures, but the required intensity of application differed, as well as the volumes. Roughness smoothing required least material volume displaced, followed by main embankment raising, groyne and minor embankment lowering. Embankment relocation using alpha shapes required the largest material volume at a 0.05 m flood level reduction. Conversely, an ambitious target of 0.5 m flood stage reduction could only be reached using floodplain lowering, main embankment raising, and relocating the main embankment. Note that the difference between the displaced material volume of main embankment raising and the increased water volume is a factor 10 for our study area. The factor depends on the mean width of the cross sectional area, but Fig. 8. Example of modelled flow velocity differences from the side channel measures. The inset shows the actual location of a newly dug side channel (in light blue), which is quite close to the predicted location for a side channel in this floodplain (Fig. 6E). (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) Fig. 9. Width averaged flow velocities (red lines, u (m/s)) and flow velocity differences with the reference situation during design discharge (blue lines, (m/s)) for six measure types. Flow velocities differences within the main channel, at 25% (left side) and 75% (right side) of the main channel, are given in a grey fill. All velocities are averaged over 100 m in longitudinal direction. (For interpretation of the references to color in this figure legend, the reader is referred to the Web version of this article.) this clearly shows that embankment raising is an effective method. Surprisingly, the lines of material and water volume for embankment relocation cross each other at a flood level reduction of around 0.1 m, and the material volume increases by 2.5 for larger alpha shapes. Many small embankment relocations require a large material displacement, and are ineffective for flood hazard reduction.

Discussion
With RiverScape, primary geospatial data can be used to (1) quickly update a hydrodynamic model and (2) determine the two dimensional hydrodynamic effect of landscaping measures. This modeling pipeline provides a transparent data stack for a systematically modelled set of specific measures at a range of intensities. The results can be seen as endmembers of river management options as each of the measures is assessed in isolation. The ranking of the measures in terms of their potential in mean flood level lowering over the whole study area (Fig. 11) is as follows: minor embankment lowering (0.04, or 0.11 m depending on full removal of fixed weirs), floodplain smoothing (0.13 m), groyne lowering (0.04 m, or 0.20 m for full removal), side channels (0.18 m), floodplain lowering (0.92 m), and embankment relocation (1.23 m). The RiverScape routines provide flexibility in the area of application: per river section, per floodplain section, or over the whole reach as presented in this paper. Also the positioning and parameterization settings (Table 3) can quickly be adjusted in intuitive ways to create new measures.
Our tool can be applied to all alluvial rivers, provided the input data are available. Results depend on the initial land cover, the bathymetry, and the position of the main embankments relative to the main channel. However, the major alluvial rivers in densely populated areas share many characteristics with our study area. Floodplains of the Mississippi, Yangtze, Elbe, Danube and San Joaquin Rivers all have floodplains that are (1) around five times the width of the main channel, (2) largely comprise land cover with a low vegetation roughness, and (3) are fixed in position with groynes or riprap. Therefore, we believe that the ranking between flood hazard reduction and volume of displaced material (Fig. 11) can be upscaled or downscaled with main channel width and discharge. The precise relationship in other areas would require additional modeling. Our results are expected to scale less well with freely meandering rivers, such as the Red River (USA), Paran a River (Brazil), or the Guaviare River in Colombia as the natural land cover differs strongly from the River Waal.
The applicability of the tool depends also on the availability of data and a hydrodynamic model schematization. The minimum requirement is a 2D hydrodynamic model, a land cover map, and a digital terrain model. However, for most countries where the problems of reconciling river functions are urgent, these data are likely to be available in some form. In addition, global hydrodynamic models are getting more detailed by using remote sensing data and open data, e.g. OpenStreetMap (Schellekens et al., 2014). RiverScape now uses 14 layers from a ready-made geodatabase, but most layers could be derived from a hydrodynamic model. For example river axis, river kilometer, left and right shoreline (Table 1) could be derived from the flood extent at low flow. Roughness codes could be extracted from a global land cover dataset (Chen et al., 2015) and floodplain lakes and channels can be derived from global scale permanent water body products (Pekel et al., 2016). The preprocessing for RiverScape would need to be tailored to these inputs. Similarly, we used Delft3D Flexible Mesh as our 2D flow model, but interoperability with other models could be created due to the modular setup of RiverScape. This would require a suitable conversion script for each 2D flow model (e.g. LISFLOOD-FP, Telemac, Mike21, TUFlow). With Python as the scripting language, the data preprocessing and updating the 2D flow model input is likely to be possible.
Within the framework of decision support, we focused on two physical criteria for the evaluation of all measures: transported material and the flood level lowering. The results should be interpreted as the exploration of the parameter space. More detailed measures should be designed by landscape architects in combination with engineers and stakeholders. In reality, measure selection includes a number of additional parameters, which were outside the scope of this study. These aspects are both limitations of the current study and possibilities for other applications and model extensions in the future as follows.
Firstly, in the middle and upper reaches of the river, the measures should not increase the flood wave celerity to avoid adverse downstream effects (Hooijer et al., 2004). To determine the effects on the propagation of the flood wave, the stationary design discharge in this study should be replaced with a standardized flood wave with a peak discharge equal to the design discharge. This would particularly be useful in steeper and longer river reaches, while our study reach is situated at the downstream end of a large river, meaning that flood waves are relatively low and long. Measures that increase the flow velocity and the fractional discharge through the main channel increase the flood wave celerity (Jansen et al., 1979), such as roughness lowering, and groyne lowering. In contrast, floodplain lowering, embankment relocation and side channel recreation will slow down the flood wave. Flood wave celerity and attenuation could provide additional parameters for measure evaluation outside of the downstream river reaches. This is useful to quantify in a longer reach situated more upstream in a river system. A major flood wave for the River Waal takes around two weeks, which is much longer than the travel time through the area of around 12 h. In the upper parts of the catchment the flood wave length and travel time differ less and the effects on the flood wave propagation are therefore stronger.
Secondly, the set of measures could be extended with deepening of the main channel in addition to measures in the floodplain and in the groyne field. Main channel deepening is technically similar to floodplain lowering: create a mask for the area to be lowered and apply a change in bathymetry over the masked area, either as a fixed value, or as spatially distributed values. Contrary to floodplain lowering, channel deepening would require a relative change in bathymetry rather lowering to an absolute external level based on exceedance levels. For the Waal, also the permanent layers should be taken into account that were created to reduce deep scour in sharp bends.
Thirdly, the improvement of the fluvial ecology provides a secondary objective of many flood hazard measures (Buijse et al., 2002;Bernhardt et al., 2005), which is also required by the European Water Framework Directive and the American Clean Water Act (Hering et al., 2010). The changes in ecotope composition can be evaluated beforehand on the potential biodiversity for different taxonomic groups (Lenders et al., 2001;De Nooij et al., 2004;Straatsma et al., 2017). For example, in the parameterization of floodplain lowering, we assigned 'production meadow' as the new ecotope, but its biodiversity potential is rather low. Including potential biodiversity scores would provide a more complete evaluation.
Fifthly, investment costs (Eijgenraam et al., 2017), depreciation costs from higher flood frequencies for agriculture in the floodplains, and maintenance costs would provide insight in the financial feasibility. The investment costs include cost for earthwork, treatment or storage of polluted soil, dike raising, groyne lowering, and acquisition and/or demolition cost of buildings and land parcels.
Sixtly, under natural land management vegetation succession leads to a shift in vegetation from meadows or agricultural fields to herbaceous vegetation, shrubs and forest over a period of decades. The associated increase in hydrodynamic vegetation roughness lowers the conveyance capacity of the floodplains and increases the water levels (Makaske et al., 2011). The model could be extended with a vegetation succession model. The resulting time series of vegetation distributions require the conversion to trachytopes or roughness value to serve as input for a two dimensional hydrodynamic model.
Lastly, the durability of the measure can also affect the selection. Embankment relocation has a long lasting effect on the flood levels. In contrast, roughness lowering can be reversed in years due to vegetation succession, and floodplain lowering can be undone in decades due to increased sediment deposition (Baptist, 2005;Makaske and Maas, 2007). The seven extensions would add to Fig. 11. Mean flood stage reduction at design discharge as a function of material displacement volume and newly created water volume at design discharge for seven landscaping measures and six intensities of application. When the material volume equals the water, the lines are superimposed in the figure and the dashes are invisible. Black dots indicate the modelled measure intensities.
the completeness of the decision support, but can not replace stakeholder interaction with experts to make the final decision.
In the Netherlands, flood risk management has been based on a design flood with an average return period 1250 years (10,165 m 3 s -1 for the River Waal). For the future, we should consider the new risk-based approach, which takes the economic value of the protected hinterland and the number of lives at risk into account in designing the type, size and location of flood protection measures (Broekx et al., 2011). This should include a cost-benefit assessment of the measures that addresses the flood risk of the protected area as well (Brouwer and van Ek, 2004). In addition, it might be required to increase the conveyance capacity of the River Waal from 10,165 to 11,436 m 3 /s (Te Linde et al., 2010), which can be combined with ecological restoration. Which measures should be carried out and in what order? A large scale setback of the main embankment clearly had the strongest effect on flood levels. Should we continue with small-scale measures in the floodplain, or would it be more cost-effective and ecologically superior to set back the embankments, and rebuild houses on raised mounds? Making these trade-offs visible could revitalize the public debate on flood-proofing river and delta systems. Much of the time during the planning of flood alleviation measures is spend on negotiations between stakeholders, decision making processes and exploration of alternatives, which is typically done in a sequential manner rather than in parallel. There is no guarantee that the outcome of the policy arena will be the same if the process is repeated. Evolving decisions depend on timing, the individuals involved, interdependencies between people and organizations involved, and the larger context we are operating in. A systematic inventory of intervention options and their costs and benefits (hydrodynamic, financial costs, biodiversity, ecosystem services) would provide high-dimensional feature space that makes the choices transparent and numerically underpinned. Pareto-optimal solutions in this feature space represent the rational, numerically optimized cost-benefit solution, which can be compared against solutions driven by desires of specific stakeholders, or political optimization. Given a fully operational and interactive tool, we believe the decision process can be shortened by years. This approach turns the usual planning process around by starting at the effects, and working backwards towards the spatial implementation, thereby transferring more significant information on flood, cost, and ecology while fundamental information on geometry is still available. Nonetheless, the modeling results should always be interpreted by an expert panel.

Conclusions
In this study, we aimed to (1) develop a tool to automatically position and parameterize seven flood hazard reduction measures and (2) evaluate these measures on hydrodynamic effects plus the required volume of displaced material. The RiverScape toolbox automatically positions and parameterizes typical landscaping measures and updates a hydrodynamic model accordingly. The ranking of flood hazard reduction in terms of transported material from high to low effectiveness is: vegetation roughness smoothing, main embankment raising, groyne lowering, minor embankment lowering, side channel construction, floodplain lowering and relocating the main embankment. This provides an integrated assessment at river-reach scale rather than many disconnected measures for individual floodplains as is the current practice for the study area. Water level reductions of more than 0.5 m could only be achieved with floodplain lowering, or embankment relocation for the Waal River. The modelled reduction in flow velocities in the main channel served as a proxy for morphological tendencies, which suggested that the trend of ongoing bed degradation could be slowed down due to lower flow velocities in the main channel. We applied all measures in isolation to determine the endmembers of river management options. However, the routines are flexible in their application. Spatial subsets could be used for local planning, or a combination of measures could be tested to optimize specific solutions with respect to biodiversity, or long term flood safety. Given a fully operational and interactive tool and an expert panel to interpret the results, we believe the decision process can be shortened by years.