2.1 Numerical model

2.2 Topographic data and spatial resolution

When levee-breach-induced floods in lowland areas are modelled, a high-resolution mesh is often necessary to describe the relevant terrain features typical of man-made landscapes (e.g. roads, railways, channels and embankments). High-resolution topographic information for large areas can be obtained from new remote sensing techniques, such as light detection and ranging (LiDAR) and the Shuttle Radar Topography Mission (SRTM), which provide raw data for the generation of digital terrain models (DTMs). LiDAR-based DTMs are nowadays available for public access for most flood-prone areas in Europe, and meshes derived from these data (even coarsened) often provide the most accurate results for flood modelling (Md Ali et al., 2015).

DTM grids can include billions of cells; however, the number of cells (and thus the runtimes) can be decreased by performing a moderate downsampling (e.g. reducing the grid size from 0.5–1 to 2–5 m) and by adopting non-uniform meshes, either unstructured (Liang et al., 2008; Sætra et al., 2015) or structured (Vacondio et al., 2017). In fact, in urban and suburban areas, the presence of road and railway embankments can influence the flood dynamics significantly, and the bathymetry near these elements should be at a high resolution (2–5 m). On the other hand, for uniform rural areas a lower resolution (e.g. 10–50 m) can be used without impairing the overall accuracy.

It is relevant to note that high-resolution DTMs can be exploited thanks to the availability of parallel 2-D codes; until a few years ago, traditional 2-D and 1-D–2-D models usually adopted a low resolution of the order of 50-100 m for flood-prone areas in order to reduce the computational times (Aureli and Mignosa, 2004; Aureli et al., 2005; Vorogushyn et al., 2010; Masoero et al., 2013; Mazzoleni et al., 2013; Huthoff et al., 2015).

2.3 Upstream and downstream boundary conditions

Discharge hydrographs with a specific return period are assigned as the upstream boundary condition. Sometimes these hydrographs are already available from previous hydrological studies and can be provided by local river basin authorities; otherwise, they can be derived from rainfall–runoff modelling or from statistical analyses of recorded discharge hydrographs (e.g. Tomirotti and Mignosa, 2017).

For the purpose of this study, multiple hydrological conditions should be considered in order to cover possible configurations characterized by different breach triggering mechanisms, flood volumes, etc. At least two different discharge hydrographs should be considered for each breach location even though the simulation database can be extended with more hydrological inputs if needed. The first case (“inflow A”) corresponds to the condition for which the water surface elevation reaches the levee crest somewhere along the river, thus generating overtopping. The second configuration (“inflow B”) concerns a flood event with a lower return period, for which the levee is never overtopped, but other mechanisms might induce the levee collapse. In fact, earthen levees can also experience breaching for piping and internal erosion processes, even when water levels remain below the levee crest. Besides this, the dens of burrowing animals (e.g. porcupine, badger and nutria) were recently identified as another possible cause for breach triggering (Viero et al., 2013; Orlandini et al., 2015). Incidentally, the collapse of an embankment during a flood event with a relatively low return period can be very threatening for human lives because the highest warning thresholds may not be reached, and the population can be unprepared to face flooding.

The choice of the return period for inflows A and B is river-dependent because it is influenced by the design return period of the levee system, by the presence of other flood control structures, etc. In general, preliminary simulations of the propagation of flood waves with different return periods (e.g. 10, 20, 50, 100, 200 and 500 years) for each river should be performed, and the event that corresponds to incipient overtopping can be identified as inflow A. Then, a higher frequency hydrograph can be selected as inflow B in order to consider levee collapses due to piping or other mechanisms during a flood event that does not induce overtopping (for example, when a specific freeboard is guaranteed).

The discharge hydrographs thus obtained are imposed as the upstream boundary condition for the levee-breach scenarios. The downstream boundary condition also deserves special attention. In fact, downstream of the breach location, the discharge may reduce or even fall to zero and reverse, but the water depth does not necessarily reduce accordingly. Hydraulic gradient and inertia can play a significant role, and the stage–discharge relationships at downstream sections may be characterized by a strong loop effect (D'Oria et al., 2015). This leads to the suggestion that, if a single-valued rating curve is imposed as outflow boundary condition, it should be assigned at the farthest possible section downstream from the breach location (even if some criteria concerning the optimal location of a downstream unique rating curve were studied; e.g. Singh et al., 1997).

2.4 Levee-breach location and modelling

Several breach locations must be identified along the river levees so that a real event can be associated with one of the simulated scenarios. The distance between two consecutive breach positions should be selected considering the river dimensions, the presence of urban settlements, the flood-prone area topography, and the possible presence of roads and embankments influencing the inundation dynamics.

The description of the gradual opening of the levee breach must be somehow included in the 2-D modelling, since the hypothesis of instantaneous break is not realistic for river embankments. Among the available approaches in the literature, which also include the coupling of the SWEs with a sediment transport model (Faeh, 2007), or with an erosion law (Dazzi et al., 2019), the simple “geometric” approach is selected in this work for two reasons. First, the uncertainties about the geotechnical parameters of the embankment and the complexity of the breaching process (three-dimensionality, interactions between erosion, infiltration, bank stability, etc.) can be neglected. Second, this method was successfully applied to a real test case (Vacondio et al., 2016), and a similar approach is often adopted for 1-D and 2-D coupled models, especially in the context of flood hazard assessment (e.g. Vorogushyn et al., 2010; Mazzoleni et al., 2013). Thus, in the RESILIENCE project the breach evolution is modelled by varying the local topography in time assuming a trapezoidal opening that widens and deepens from the crest of the embankment to the ground elevation; the breach geometric dimensions (i.e. width) and opening time are given as input parameters.

The two breach parameters should be set consistently with historical data, if available (e.g. Nagy, 2006; Vorogushyn et al., 2010; Govi and Maraga, 2005), or hypothesized according to the river and embankment characteristics. A breach final width of the order of tens to hundreds of metres is often assumed in the literature (Apel et al., 2006; Kamrath et al., 2006). As for the breach development time, very few field observations are available, and values in the range of 1–3 h are often reported (Apel et al., 2006; Alkema and Middelkoop, 2005). The impact of the parameters' uncertainty on the results should be evaluated for at least one hydrological scenario, by means of a probabilistic treatment (Apel et al., 2006; Vorogushyn et al., 2010; Mazzoleni et al., 2013) or a sensitivity analysis (Kamrath et al., 2006; Huthoff et al., 2015).

For each breach location and return period, the triggering time for breach opening should be selected after preliminary simulations and corresponds to the moment when either overtopping or the peak water surface elevation is observed at the breach location.

Finally, simulations must be extended in time long enough to capture the flooding evolution in the lowland, which is actually the goal of the presented methodology. The selection of this temporal interval should be guided by considerations of the flood duration in the river, of the inundation dynamics, and of the fact that drainage operations and/or breach closure would start at some point.

2.5 Outputs

The outcomes of the modelled scenarios could help the civil protection activities for emergency planning and/or at the occurrence of an event similar to one of those already modelled. Arrighi et al. (2019) recently presented a framework that integrates hydrologic–hydraulic modelling with human safety and transport accessibility evaluations, applied to a small municipality near Florence (Italy), and confirmed the usefulness of detailed spatial and temporal flood data provided by numerical modelling for civil protection purposes. Thus, the first mandatory output concerns spatial and temporal information about the flood dynamics in the lowland area, which can be visualized as an animation of the inundation pattern or as a sequence of frames at selected times.

Other useful indicators for evaluating the flood severity for each scenario include the maps of arrival times of the wetting front, maximum water depths (and/or water surface elevations), maximum velocities, and a hazard index which combines simultaneous water depth and velocity. In general, these maps can be opened in a GIS environment and overlapped with layers of data (e.g. about population, transportation, buildings, critical and vulnerable structures, etc.) to analyse the possible flood impacts on the territory, with the aim of increasing the resilience in the area. First, information about the maximum water depth allows for the definition of the affected areas as well as the identification of escape routes and safety zones where assembly points can be organized. On the other hand, in areas where only shallow water levels can be expected, people can simply be instructed to protect their homes with anti-flood barriers to prevent water from damaging their property. Besides this, analysis of the inundation dynamics can reveal possible service interruptions and disturbance to road accessibility, which can also prevent rescue operations; in this way, alternative routes can be identified.

Moreover, the maximum flow velocity map should not be neglected because high velocities can reduce the stability of vehicles and pedestrians, increasing the hazard for human lives (Milanesi et al., 2015). In general, the most important results for quantifying the hazard for human lives are the maximum simultaneous water depth and velocity, which can be represented in terms of the Froude number, total head, total force and/or total depth. This last indicator, which represents the water depth at rest D whose static force is equivalent to the total force of the flow, is evaluated as follows (Aureli et al., 2008; Ferrari et al., 2019):

where h(t) represents the water depth and Fr(t) the Froude number at time t.

Finally, the map of the arrival times of the wetting front can be useful for early warning and for establishing the timeline for the evacuation procedures, in particular for vulnerable people (older adults, hospitalized patients, etc.). The offline analysis of the simulation results, on the other hand, can help in the identification of possible strategies to reduce the damage, as for example the adoption of movable defence systems (e.g. flood barriers) to preserve lowland urban settlements from flooding or other emergency interventions to drain water (e.g. pumping and relief cuts). Selected strategies can also be tested numerically to verify their effectiveness.

3.1 Set-up

The computational domain was built based on the available 1 m resolution DTM of the riverbeds and of the floodable lowlands derived from LiDAR surveys carried out in the years 2008, 2015 and 2016. In order to avoid the excessive memory requirements and heavy computational costs (even for a fast GPU-parallel model), related to the adoption of a uniform 1 m mesh (which would require billions of cells), the DTM was downsampled to a resolution of 5 m. This operation did not affect the crest elevation of the artificial embankments. In fact, each 5 m×5 m cell crossed by an embankment was identified, and its elevation was set equal to the maximum value among the original 25 points belonging to that cell; otherwise, its elevation was simply computed as the average of the 25 terrain data comprised in that cell. For other urban features that were not captured correctly by the LiDAR survey, additional corrections were introduced manually.

Then, the study domain was discretized by means of a non-uniform BUQ grid: the maximum resolution (5 m) was forced along the riverbed, the levees, the main artificial embankments and channels, and at the breach location, whereas for rural areas it gradually decreased up to 40 m according to the algorithm presented by Vacondio et al. (2017). The resulting computational grid (Fig. 2), whose spatial resolution is considered suitable for the detailed modelling of the river and the lowland area, consists of roughly 13×10⁶ cells, and the number of cells is reduced by approximately 70 % compared to a uniform 5 m mesh.

Figure 2Multiresolution computational grid for the pilot area.

The study area also includes several urban settlements: Modena, with about 185 000 inhabitants, is the most populated one, and its old city centre comprises narrow streets, which cannot be accurately described with a 5 m resolution. Therefore, limited to the few scenarios concerning the flooding of Modena, simulations were performed using a finer mesh (up to 2.5 m) in the town, and the urban layout was modelled according to the “building hole” method (Schubert and Sanders, 2012), in order to capture the flow field inside the urban area correctly.

Buildings were explicitly resolved only for the city of Modena, whereas, according to previous findings (Vacondio et al., 2016), the presence of the other (smaller) settlements was taken into account by means of a higher resistance parameter (“building resistance” method; Schubert and Sanders, 2012). In particular, the roughness coefficient for the urban areas was calibrated based on the event that occurred in 2014 (the flood arrival times at selected locations were known) and was set equal to 0.14 m${}^{-\mathrm{1}/\mathrm{3}}$ s, while for rural areas a Manning coefficient of 0.05 m${}^{-\mathrm{1}/\mathrm{3}}$ s was chosen. In the absence of data for calibration concerning past flooding events, land use maps can be exploited to assign standard roughness values from the literature.

As for the river roughness, the calibration for the Secchia River was performed in previous works (Vacondio et al., 2016) based on the water levels recorded at the available gauging stations during recent flood events. A Manning coefficient equal to 0.05 m${}^{-\mathrm{1}/\mathrm{3}}$ s provided the best results. The roughness of the Panaro River was also subject to calibration with a similar procedure, and the value 0.04 m${}^{-\mathrm{1}/\mathrm{3}}$ s was selected.

With regard to the selection of the hydrological scenarios for the two rivers, the synthetic design hydrographs (Tomirotti and Mignosa, 2017) with assigned return periods were considered. After preliminary simulations, inflow A (overtopping-induced flooding) was identified as the 50-year return period hydrograph for the Secchia River and as the 100-year one for the Panaro River. Then, the inflow hydrograph with a 20-year return period was selected as inflow B for both rivers in order to consider an event with higher frequency. These discharge hydrographs were assumed as upstream boundary conditions during the simulations.

The downstream section was located at the confluence (of the Secchia and Panaro rivers) into the Po River, and a constant water level in this (much larger) river, which did not affect the breach outflow even for the most downstream breach location, was assumed as the boundary condition.

As mentioned before, the pilot area can be flooded by hypothetical failures occurring along the Secchia and Panaro rivers. Assuming a distance of about 2 km from one breach to the other, 30 breach locations were identified along the right levee of the Secchia River (for the flooding scenarios related to this river) and 26 ones along the left levee of the Panaro River. Figure 1 reports all the 56 simulated breach sites. Among these, eight breach scenarios on the Secchia River and four on the Panaro River involve the city of Modena. Based on past observations, the breach final width was assumed equal to 100 m for all scenarios, while the opening time was set equal to 3 or 6 h for inflows A and B, respectively.

All the 112 simulations (56 breach locations and two hydrological scenarios) were prolonged for 72 h (3 d) because at that point the outflow from the breach was almost null, the flooded area had reached its maximum extension and no significant flow dynamics can be observed.

3.2 Outcomes

Figure 3 shows an example of the results for one scenario on the Secchia River. The maps of maximum water depths shown in Fig. 3a reveals that the flooding involves the northern portion of the domain, partially affecting some urban settlements (S. Possidonio and Mirandola), while villages in the east are safe from this inundation and could temporarily host the evacuated population. The maximum velocity magnitude (Fig. 3b) remains below 1 m s⁻¹, except for the surroundings of the breach and close to some road embankments that are overtopped by water (see detail in Fig. 3b), where drivers can be in grave danger. The combination of water depth and velocity (Fig. 3c) highlights that lowland areas are mostly affected with low (green $\mathrm{0}\le D<\mathrm{0.5}$ m) and medium (yellow 0.5 m ≤ D<1 m) total depth values, even if higher values (orange 1 m ≤ D<1.5 m and red D≥1.5 m) are reached where high water depths are observed. The map of flood arrival times in Fig. 3d shows that, for this scenario, approximately 45 % of the flooded area is affected 24 h after the breach opening (purple contour line), guaranteeing a considerable amount of time for emergency activities. It is relevant to note that also during the levee breach occurred in 2014 on the Secchia River (Vacondio et al., 2016) one of the most affected villages was flooded the day after the opening, but no countermeasures were taken at that time, since a priori knowledge of the inundation dynamics was not available. Finally, a video showing the flooding evolution for this scenario is provided as additional material.

Figure 3Example of the resulting maps concerning the maximum (a) water depth, (b) velocity, (c) total depth and (d) flood arrival time for a given scenario on the Secchia River with inflow B. The main roads, railways and urban settlements are identified, and the breach location is indicated with a black cross. The topographic information in background is made available by Emilia-Romagna region at https://geoportale.regione.emilia-romagna.it/it (last access: 8 January 2020).

In addition to the maps representing specific hydraulic indicators for each scenario, further information can be obtained by analysing the results of the whole database. In particular, the most affected parts in the pilot area and those never hit by flooding can be investigated. Therefore, for both inflows A and B, the maps of the inundated areas were combined in order to quantify the number of scenarios affecting each computational cell in the domain. Figure 4 shows the resulting map for inflow B: in this case, about 50% of the pilot area is affected by at least one breach scenario. In particular, two areas can be identified as being most affected by the possible flooding induced by levee breaching, since up to 21 breach scenarios (from both the Secchia and the Panaro River) involve these regions. On the other hand, it can be observed that there is a large zone, in the middle of the pilot area, that is not affected by any of the considered breach scenarios thanks to the favourable terrain topography; hence, it is potentially recommended for evacuation purposes (e.g. organization of assembly points).

Figure 4Number of flooding scenarios affecting each cell of the pilot area with inflow B. Two portions of the domain (in red) are flooded for 21 scenarios (from either the Secchia or the Panaro River), while the uncoloured zones are never inundated. The topographic map in background is made available by Emilia-Romagna region at https://geoportale.regione.emilia-romagna.it/it (last access: 8 January 2020).