Inundation mapping based on reach-scale effective geometry
- 1Irstea, UR HYCAR, 1 Rue Pierre-Gilles de Gennes, 92160 Antony, France
- 2Sorbonne Universités, UPMC Univ Paris 06, CNRS, EPHE, UMR 7619 Metis, 4 Place Jussieu, 75005 Paris, France
Correspondence: Cédric Rebolho (firstname.lastname@example.org)
The production of spatially accurate representations of potential inundation is often limited by the lack of available data as well as model complexity. We present in this paper a new approach for rapid inundation mapping, MHYST, which is well adapted for data-scarce areas; it combines hydraulic geometry concepts for channels and DEM data for floodplains. Its originality lies in the fact that it does not work at the cross section scale but computes effective geometrical properties to describe the reach scale. Combining reach-scale geometrical properties with 1-D steady-state flow equations, MHYST computes a topographically coherent relation between the “height above nearest drainage” and streamflow. This relation can then be used on a past or future event to produce inundation maps. The MHYST approach is tested here on an extreme flood event that occurred in France in May–June 2016. The results indicate that it has a tendency to slightly underestimate inundation extents, although efficiency criteria values are clearly encouraging. The spatial distribution of model performance is discussed and it shows that the model can perform very well on most reaches, but has difficulties modelling the more complex, urbanised reaches. MHYST should not be seen as a rival to detailed inundation studies, but as a first approximation able to rapidly provide inundation maps in data-scarce areas.
Floods are a recurring phenomenon in France: in September 2014, intense rainfall affected the south of the country, leading to several deaths and about EUR 0.6 billion worth of damage. The following year, in October, about 20 people died in the south-east due to massive flooding, which caused a loss of EUR 0.5 billion. Then, in June 2016, large-scale flooding occurred over the Seine and Loire catchments, mainly affecting their tributaries and resulting in four deaths at a cost of EUR 1.4 billion. These are only examples which underline the value of flood inundation mapping to anticipate the impact of such events. Public authorities and insurance companies are showing a growing interest in the field of rapid inundation modelling, and for the development of simple methods, that would work for any river with easily available data.
Flood hazard assessment usually combines rainfall observations or simulations, a hydrological model, streamflow simulations or observations, and an inundation model in order to generate inundation extents, height maps and sometimes other information (e.g. velocities). Traditionally, flood inundation models are derived from the shallow water equations (SWEs) in one or two dimensions (the so-called hydraulic models), with various simplifications that have proved to give satisfying results. For instance, the Regional Flood Model (RFM), probably one of the most comprehensive approaches published so far, is made of four parts (Falter et al., 2014): a daily distributed rainfall–runoff model, a 1-D hydraulic model for channel routing, a 2-D hydraulic model for floodplain mapping and a flood loss estimation model. Its application on the Mulde catchment in Germany (Falter et al., 2015) showed mixed results concerning inundation extents, correctly predicting only 50 % of the flooded area for the August 2002 event. This underestimation was explained by dike breaches that were not accounted for within the model. A lack of observed data did not allow validation on other events.
Not all hydraulic models need to have this degree of complexity. It is indeed possible to neglect specific parts of the SWEs depending on the situation. Usually, 2-D models use the complete Saint-Venant equations while 1-D models often disregard one or several terms, leading, for instance, to the diffusive wave or kinematic wave approximations (Moussa and Cheviron, 2015). Some methods choose to couple 1-D and 2-D models, the former for streamflow routing and the latter for overbank flow (Morales-Hernández et al., 2016). Despite the accuracy of such models, studies often try to further simplify them because of the large computing time to simulate small areas and the lack of precise data required to run these models.
LISFLOOD-FP (Bates and De Roo, 2000), a hydraulic model developed to simulate floodplain inundation, was used in several studies (Biancamaria et al., 2009; Horritt and Bates, 2001; Hunter et al., 2005). The model offers different possibilities: using 2-D equations or 1-D equations decoupled on a 2-D grid with kinematic, diffusive or inertial approximations (Bates et al., 2010). Horritt and Bates (2002) published a comparison between different models with gradually increasing complexities (1-D, 1-D on 2-D grid and 2-D) and, surprisingly, showed that the 1-D model had a better ability to reproduce the two events that were used in validation. The subsequent analysis concluded that the reach studied was relatively narrow and could easily be modelled using simple methods, and the authors argued that the other models would be more appropriate for more complex reaches.
However, these examples concern relatively small and well-instrumented reaches and assessing flood hazard at a larger scale may require different approaches. Alfieri et al. (2014) applied LISFLOOD-ACC, an inertial version of LISFLOOD-FP with decoupled 1-D equations on a 100 m resolution grid over Europe in order to map flood hazards for a 100-year return period, assuming a constant return period along the reaches. Broadly speaking, the model splits rivers into small reaches, to apply the hydraulic models independently and to merge simulated maps together, but only for rivers with a catchment larger than 500 km2. The model was then validated against regional and national hazard maps for six catchments in Germany and the United Kingdom and showed a general over-prediction. Another variation of LISFLOOD-FP for large-scale flood inundation modelling was introduced by Neal et al. (2012), including a new sub-grid representation of channel networks for improved model accuracy (Neal et al., 2012; Schumann et al., 2013).
Le Bihan et al. (2017) developed an approach aimed at the forecasting context, in order to cope with excessive computing times. The solution chosen was to run a simple 1-D hydraulic model during a “pre-analysis phase” and create a catalogue of inundation extents corresponding to various return periods. These maps are then used, in a forecasting context, to give an estimate of the level of flooding, depending on the forecast discharge.
The lack of precise data (especially for channel cross sections) and the computing time required by numerical methods for solving the SWE motivated the development of potentially alternative methods, mostly based on DEM analysis. For instance, the rapid flood spreading method (Gouldby et al., 2008) was chosen to divide floodplains into impact zones of different elevations in order to explore the effects of dike breaches using a spilling algorithm based on water depth. Other methods derive inundation maps from topographic information only: one can cite EXZECO (Pons et al., 2010), which introduces elevation noise in the DEM in order to create a single map of “maximum flow accumulation” that can be seen as a potential inundation area, and HAND (“height above nearest drainage”), a descriptor originally used for terrain classification (Nobre et al., 2011; Rennó et al., 2008), which has recently been adapted to static flood inundation mapping (Nobre et al., 2016) and is increasingly used to produce flood maps (Afshari et al., 2018; McGrath et al., 2018; Speckhann et al., 2018). HAND calculates the difference between river cells' elevation and that of the connected floodplain cells, thus giving relative height information which can be compared to observed flood depths and the corresponding inundation extent.
MHYST, the method presented in this paper, is a simplified approach developed with the aim of rapidly producing inundation maps in data-scarce areas. It combines (i) concepts of hydraulic geometry to characterise channel geometry and (ii) DEM-derived relative elevations to characterise the floodplain; it does not work at the cross section scale but computes effective geometrical properties representative of the reach scale. Combining reach-scale geometrical properties with simplified steady-state hydraulic laws allows one to rapidly generate flood inundation maps while ensuring reach-scale coherence. After describing the method and the calibration dataset, MHYST is compared against the inundation extent observed for the major event that occurred in May–June 2016 in France. The last section discusses the spatial distribution of performance and the impact of uncertainties on the results obtained.
The MHYST model stands for Modélisation HYdraulique simplifiée en écoulement STationnaire, i.e. Simplified Steady-state Hydraulic Modelling. It is a flood inundation model which aims to map inundation extents at the reach scale. Where classic hydraulic models use cross sections, this method is based on an effective geometry representative of each river reach. Since no detailed geometric data were available to describe the shape and roughness of the channel river bed for this study, a sub-grid representation of the channel was derived from hydraulic geometry relationships linking the drainage area with bankfull width and height (Leopold and Maddock, 1953). When discharge exceeds bankfull capacity, the model computes a reach-scale relation between streamflow and the HAND defined by Nobre et al. (2016). This relation can finally be used to assess which height corresponds to the given streamflow, and thus to derive the corresponding inundation map.
2.1 Processing of DEM: from elevations to height above nearest drainage
The initial step consists of processing the digital elevation model (DEM) in order (i) to obtain a flowing drainage direction map, (ii) to identify the subcatchments (corresponding to the river reaches), and (iii) to compute the height above nearest drainage (HAND) in each subcatchment. This initial processing is the basis of the floodplain analysis in MHYST. To compute the drainage direction map, we used the D8 method from the Flow Direction function provided by ArcGIS 10.3. It computes the drainage direction by calculating the steepest slope from the eight possible directions for a given cell.
Figure 1 shows the procedure used to compute HAND values: for a given floodplain cell, it is the difference between its elevation and that of the closest river cell in terms of drainage direction. For instance, the cell of elevation 25 at the top of the figure is linked to (i.e. flows towards) the most upstream red river cell which has an elevation of 18: thus, its HAND value is . This relative height has been used as a proxy for inundation height by various studies (Afshari et al., 2018; Nobre et al., 2016). To derive an inundation map from HAND values, we must define a threshold height HT: the flooded area corresponds to all the cells whose HAND value is strictly lower than HT (Jafarzadegan and Merwade, 2017).
2.2 Model description
MHYST is mostly based on a DEM and its derivatives (drainage map and drainage areas) and on the hydraulic equations describing a steady uniform flow at the reach scale. This means that for a given time step (day in this case), at a given reach, we make the approximation that the flow is constant over time and space (this is obviously a strong simplification that we will discuss later). Table 1 sums up the variables used in the following equations as well as their respective units and interpretations. Table 2 describes the two free parameters of the model. The following equations show the path to build a reach-scale relation between HT and the streamflow Q by calculating, with hydraulic formulas, the discharge value corresponding to a given HT. Once this relation is known, the model can easily simulate a hydrological event by inverting the relation, and by searching for the HT which corresponds to the given Q (Fig. 2).
Other variables can be directly calculated from the DEM (Fig. 4): for a given threshold height HT at a reach of length L (L is a fixed parameter of the model), V(HT) is the sum of volumes above all flooded pixels and S(HT) is the area occupied by the flooded cells. A(HT), in Eq. (1), is the average cross section area of the flooded reach and it depends on V(HT) and on the bankfull cross section area of the channel (, Fig. 3). This variable can also be defined as the sum of the channel cross section area and the floodplain cross section area . B(HT), in Eq. (2), is the average surface width of the flooded reach, defined similarly from S(HT) and L.
The only unknown variables in these equations are sub-grid parameters hb (bankfull water level) and Wb (bankfull width), i.e. the bankfull geometry, which cannot be obtained from usual DEMs and are only available from detailed surveys for a small number of rivers. Indeed, in situ bathymetric data are quite scarce and red lasers cannot penetrate the water surface more than a few centimetres, which means that the real elevation of the river bed is mostly not correctly represented in the DEMs. This is why we chose to use downstream hydraulic geometry equations to estimate these geometric parameters, assuming a rectangular channel, the size of which depends on the upstream drainage area (Eqs. 3 and 4). To assess the coefficients α and β, we used satellite images from the French platform Géoportail in order to link observed bankfull widths and drainage areas. The values found for the Loing catchment are α=0.053 and β=0.822. The other coefficients, δ and ω, were taken from a study by Blackburn-Lynch et al. (2017), which attempted to regionalise these parameters in the US. We used the general values found for the whole set of catchments: δ=0.27 and ω=0.21. Although these values probably add uncertainties in the model, they are an accessible way to assess bankfull channel geometry and could still be improved by local bankfull studies when available.
The fundamental equations of the MHYST model come from an experimental study by Nicollet and Uan (1979) which defines the DEBORD formulation as in Eqs. (5) to (7). Building on the Manning–Strickler formula, these authors proposed an empirical parameterisation of turbulent momentum exchange between the channel and the floodplain. This formulation expresses the conveyance capacity depending on channel-related and floodplain-related variables. The coefficient C takes into account the interaction of flows between the fast-flowing channel and the slow-flowing floodplain, as well as the corresponding head losses.
The streamflow Q is finally defined from the conveyance capacity and the channel slope, since we hypothesise a uniform flow. Rch and Rfp can easily be calculated from the assumed reach geometry (Eqs. 8 and 9), which only leaves the Strickler coefficients as unknown variables.
The two Strickler coefficients add 2 degrees of freedom, and Kch is additionally used to calculate the bankfull flow from the Manning–Strickler formula (Eq. 10).
Here, we sum up the procedure, which operates at the reach scale:
For a given threshold height HT, we use the DEBORD formulation to calculate the corresponding discharge Q.
By repeating the operation for all possible HT, we obtain a reach-specific table matching values of HT and Q.
When working on an event where only Q is known, when it is greater than Qb (which means that the river overflowed), the model looks for the corresponding HT value in the table, by calculating a linear interpolation between two values if necessary, and then assigns to each cell in the subcatchment a flooded height .
Although this method and that of Zheng et al. (2018) were developed independently, they share a lot of similarities, both using HAND to derive a reach-scale geometry which is used as input for a simplified hydraulic model. However, in addition to HAND, MHYST uses downstream hydraulic geometry relationships to evaluate a sub-grid representation of the channel geometry. The hydraulic model is also different: Zheng et al. (2018) use the Manning–Strickler formula, while MHYST computes streamflow values from the DEBORD formulation.
2.3 Boundary conditions
MHYST can work with either simulated or observed flows. In this paper, observed data from 12 measurement stations of the French HYDRO database (Leleu et al., 2014) were used to create an observed distributed streamflow map by interpolating flows based on drainage area (Eq. 11) for river pixels between outlets:
where Q and AD are the streamflow and drainage area of any river cell between two outlets, Qup, Qdown, AD,up and AD,down are the direct upstream and downstream outlet discharges and drainage areas. This way, streamflow is coherently interpolated over the network, and then averaged at the reach scale.
3.1 Generic data
In this study, we used a 5 m resolution DEM with a vertical resolution of 0.01 m covering the Loing catchment (Fig. 5) from IGN (the French national institute for geographic information), which was filled and corrected to avoid depressions and to allow a strict coherence of flow directions, meaning that every pixel flows to the sea. Drainage directions and areas were derived from this DEM and used as model inputs along with elevations. The adaptations and modifications of the DEM were conducted using ESRI ArcGIS 10.3.
Daily observed discharges were obtained from the French HYDRO database (Leleu et al., 2014) and the stations were used to delineate the hydrological network over the catchment. Calibration data for the Loing catchment were obtained from the activation EMSN028 of the Copernicus Emergency Management Service (© 2016 European Union). The original Copernicus study covered a small part of the River Seine and half of the Loing catchment (Fig. 6). However, since the study area and the defined river network were smaller, we cropped the inundation extent to match the study area (Fig. 6). These calibration data are post-processed observed data, meaning that the original maps came from satellite observations but they were then modified to build a more homogeneous inundation extent, i.e. nearby areas whose elevations were below the observed flood level were added to the inundation extent and merged with all the others. The maximum flood extent was then validated by the European Service against reported flood damage and hydrological measurements (SERTIT, 2016).
3.2 Event of May–June 2016
Following an extremely wet month of May (namely the wettest on record for many stations), a heavy rainfall event started on 30 May 2016 over the centre of France, affecting the Upper and Middle Seine basin and the Middle Loire basin. This episode lasted until 6 June and, combined with highly saturated soils due to a series of preceding minor events, led to major flood inundations. Over this period, overall precipitation reached 180 mm in Paris and Orléans, while in some tributaries, such as the River Loing, peak flows largely exceeded those of the record 1910 flood event (Fig. 7). The flood resulted in 4 deaths, 24 people injured and EUR 1.4 billion worth of damage. A total of 1148 cities were declared to be in a state of natural disaster and insurance companies received about 182 000 claims (CCR, 2016).
Since calibration data were available for June 2016 event, we chose to use our model to simulate this episode and compare the results with observations. We conducted this study over the River Loing, tributary to the River Seine, with a catchment covering 3900 km2, a mean elevation of 148 m and a mean slope of 0.03 m m−1. This catchment was heavily impacted by the flood event and contains a significant proportion of the inundated area, making it a suitable area to carry out the study. Streamflow data were interpolated from measurements, so no hydrological model is involved in this paper.
4.1 Calibration procedure
To assess the model's performance, we used several criteria based on the contingency table in Fig. 8. These scores are presented in detail by Jolliffe and Stephenson (2003) and are defined as a ratio between members of the table where n1 is the number of hits, i.e. the number of flooded cells correctly forecast; n4 is the number of pixels correctly forecast as dry; n2 is the number of false alarms; and n3 the number of observed flooded cells missed by the model. Table 3 summarises the formulas and the interpretations of each score used in this study.
The POD (probability of detection), which is also called Correct (Alfieri et al., 2014) or M1 (Teng et al., 2015), calculates the percentage of observed inundated pixels intersected by the simulation map. Its main drawback is that it does not take into account the false alarms and thus it can give good results for a clearly overestimating inundation extent. On the contrary, the FAR (false alarm ratio) or M2 (Teng et al., 2015) computes the proportion of cells wrongly flooded by the model. But similarly, if the model does not flood anything, the FAR can reach its optimal value. The critical success index (CSI), also known as Fit, F index or FAI (Alfieri et al., 2014; Bates and De Roo, 2000; Falter et al., 2015), is a criterion which tries to give an overall performance of the simulation by calculating the percentage of correctly flooded cells above the total number of flooded cells (observed and simulated). In this way, the score is penalised by the over- and underestimation. However, this criterion does not specify if the model is over- or underestimating the observed extent. This is why we also looked at the BIAS, which computes the ratio between the number of simulated and observed flooded cells. If it is above 1, the model overestimates, and if it is below 1, it underestimates. However, a value of 1 does not equal a perfect simulation since there may be a balance between the misses and the false alarms.
These ratios are particularly reliable if they are used to compare simulations and exhaustive observations. This is almost the case with Copernicus calibration data, which represent a “maximum flood extent”. However, MHYST outputs are dated, which is not the case for the observed map. This is why all daily simulated inundation extents were merged into one maximum simulated extent, meaning that we did not try to validate the temporal dynamic of the flood, but only aimed to assess its largest area. Thus, the preceding scores will only evaluate MHYST's ability to reproduce the maximum flood extent.
MHYST has two free parameters (Table 2): Kch (the Strickler roughness coefficient for the channel) and Kfp (the Strickler roughness coefficient for the floodplains). Preliminary studies showed that, for the Loing catchment, a length of 1000 m was a good trade-off between accuracy and computation time; consequently L was fixed at 1000 m in the rest of this study. Kch and Kfp values were tested in the range [0.1;30] in order to explore a wide range of possibilities (121 combinations were tested).
To help make a decision on the optimal parametrisation of the model, we used the following graphs, on which each (Kch, Kfp) couple is characterised by one overall value:
two contour plots (Fig. 9) showing the impact of Kch and Kfp for the two main scores (BIAS and CSI);
a Pareto plot (Fig. 10a) showing the role played by the Kfp parameter in balancing the POD and the FAR;
a Pareto plot (Fig. 10b) showing that the CSI identifies the best compromises between the POD and the FAR.
Last, to be able to analyse the variability of results between reaches (we have a total of 90 reaches affected by the inundation), we also computed the CSI and BIAS reach by reach, and produced two cumulative distribution plots showing these results (Fig. 11). We found the following:
The fit criteria are very sensitive to the Kfp value and much less to the Kch value (Fig. 9): this should not be a surprise given that we deal with the maximum flood extents for calibration, where Kch only plays a minor role. Remember also that (i) we are modelling a very extreme event (T>500 years) with substantial overflowing, and (ii) that we are working with a channel geometry derived from hydraulic geometry relationships. All this contributes to making the estimation of the channel roughness coefficient more difficult.
The CSI clearly shows an optimal zone around Kfp=5 and . The best CSI values (greater than 0.66) correspond to combinations where with Kfp=5 or with Kfp=4. Given the equifinality, a good way to choose a combination in this range could be to use the most physical one, which, in this case, would be Kch=10 and . Indeed, over the catchment, floodplains mainly consist of 44 % non-irrigated arable land, 17 % broad-leaved forest and 10 % pastures with corresponding roughness coefficients reported in the literature of 8, 2 and , respectively (Grimaldi et al., 2010).
Another way to confirm the validity of this choice (Kch=10 and ) is to look at how this parametrisation behaves at the reach scale. Given a total of 90 reaches, we can compute the CSI and BIAS criteria for each of them and draw a distribution (Fig. 11): we observe that the “optimal” distribution is unbiased and that it represents a solution among the best available for each percentile, we can thus trust this parametrisation as a relatively “all-terrain” one for the Loing catchment.
Last, Fig. 10 provides a good illustration of how parameter sets interact with the FAR, POD and CSI criteria: choosing from the parameter sets with the best CSI makes it possible to find a compromise between a high POD and a low FAR.
4.3 Model behaviour
Figures 12 to 14 provide a further illustration with a colour-coded classification of each reach depending on its CSI and BIAS value. A total of 11 regions are highlighted and numbered because of their poor performance. The reasons of why MHYST was not able to reproduce the inundation extent in these regions are explained below.
For the downstream-most part of the Loing (Fig. 12), the reaches are red or orange because this area is only partially covered by the observation, which stops just after the confluence with the small tributary.
The small tributary (Fig. 12) is mainly red or orange for various reasons: downstream, at the confluence, the DEM is full of small high-elevation zones (not corrected in the DEM) which the model cannot reach, thus degrading the simulation. Along the tributary, the reason can be either the observed discharge values which seem small compared to the rest of the catchment or simply the effective geometry defined by the model, which does not correspond to the actual one. Finally, the upstream part of the tributary is not covered by the observation, which stops in the middle of what MHYST simulated. However, the study zone defined by the Copernicus Emergency Management Service goes further, so we cannot know whether it was not flooded or whether the service did not map this part because it was too insignificant.
The orange part in the middle of the BIAS map (Fig. 12) is due to the railway tracks which act like a wall in the DEM, preventing the model from reaching the other side (from east to west), where a small tributary, which looks like a partly subterranean urban stream, overflowed in its open air part.
Finally, the red and orange zones in the south of the presented map (Fig. 12) correspond to a part of the river where the Loing man-made waterway plays a major role, running parallel with the main river. This configuration is difficult for MHYST because we only consider the main river, defined by the DEM, with an effective reach-scale geometry and we cannot take into account such specificities, which would require a 2-D hydraulic model.
The area identified (Fig. 13) shows a slight underestimation leading to a moderate CSI. This issue can be explained by a motorway which is represented in the DEM by a more elevated area. This motorway separates the reach into two parts linked by artificial openings made by the producers of the DEM. This and the Loing waterway and another road act as dikes that prevent the model from reaching a further part of the reach. The parameterisation of the model is not suitable to address this difficulty.
Similarly to the previous area (Fig. 13), a railway crosses the DEM from north to south with only one opening for the water. Given the parameterisation of the model, it is not possible to go over the railway to flood the missed area.
In that case (Fig. 13), the model clearly overestimates the flood. The water fills a depression which looks like a tributary but is only a thalweg. Once more, the parameterisation of the model does not provide an adequate representation of this reach.
In this area (Fig. 14), MHYST underestimates the inundation extent due to a road that works like a dike. However, with another parameterisation, the model would be able to provide enough water to go over the road.
In the western part of the upstream area (Fig. 14), MHYST overestimates the flood because it is a relatively flat zone. The exceeding water, still due to the parameterisation, is thus spread over the area.
This area (Fig. 14) is special because the overestimation of MHYST is due to a non-continuous observation map, creating large parts of reaches that are observed to be dry. However, since MHYST works at the reach scale, it necessarily floods the whole river reach. Moreover, one tributary, the Solin, is not defined in the hydrographic network used by the model, because no observed discharges were available, whereas it appears in the observed map, leading to an underestimation of the flooded area.
The most upstream part of the simulated area (Fig. 14) suffers from an excess of water and a non-continuous observation, leading to similar effects. Moreover, several elevated roads appear in the DEM and force the model to flood the area using artificial openings across the roads.
In order to complete our interpretation of MHYST behaviour, we conducted two sensitivity analyses, one with the Morris method (Morris, 1991) and the other with the Sobol method (Sobol, 2001). We chose to assess the effect of six potential parameters, Kch, Kfp, α, β, δ and ω, that may play a major role in the computation of HT–Q relationships. In both analyses, we found that ω, which parameterises the regionalisation of bankfull heights, has the most substantial effect on the performance and that, surprisingly, Kfp has no influence at all. As a matter of fact, when we conducted the Sobol analysis with fixed hydraulic geometry parameters, we showed that Kfp is considerably more influential than Kch. We concluded that the previous results were due to the fact that these sensitivity analyses explore the parameter space in detail, and even with reasonable boundaries, they can reach values that may not be consistent with the characteristics of the catchment studied.
4.4 Influence of the DEM resolution
It is possible to assess the sensitivity to the DEM in two ways: first by aggregating our DEM from 5 m to various resolutions (10, 25, 50 and 100 m) and then by changing the source of the DEM. Figure 15 provides the CSI scores obtained by the model while changing the resolution. It shows that the resolution has relatively little effect on the optimal value, which varies between 0.65 and 0.69. However, the position of this optimal, i.e. the combination of parameters (Kch and Kfp) leading to it, changes. We can also see that for some resolutions, such as 25 or 50 m, the equifinality zone is much smaller than the one for the 100 m resolution, for example. If we also look at the “physical” set of parameters we previously identified (Kch=10 and ), we can see that the CSI reached by the model for this combination varies between the resolutions. Nevertheless, the result still seems satisfying, so it could be used as a “default” parameterisation, for instance for ungauged catchments. But this should be tested on other catchments with observed data to lead to a more comprehensive conclusion.
Before using the RGE 5 m DEM from IGN, we tried to use the 25 m EU-DEM from the European Environment Agency, and it showed poorer results, because it was not precise enough. Figure 16 shows the evolution of CSI for the same combinations of parameters as before. We see that the best combinations of parameters only lead to a 0.53 maximal CSI, which is more than 10 points below what we can obtain with the RGE DEM. There is also strictly no connection between the best values of BIAS and those of CSI, the latter being obtained for a clear overestimation of the flood extent (BIAS∼1.5). These results are due to the lack of precision of the EU-DEM, which does not distinguish the channel from the floodplain, leading to a 2 km wide channel in some parts of the river.
The objective of this paper was to present and validate a simple hydraulic model for rapid inundation mapping in data-scarce areas. MHYST is based on DEM analyses and simple hydraulic equations, creating a reach-scale relation between the average discharge and the average “height above nearest drainage” which can then be used to simulate any event, past or future, as long as streamflow information (observed or simulated) is available. This model was calibrated against an observed exceptional flood which occurred in 2016 on the Loing River near Paris and showed results that are certainly not perfect, but from our point of view and for our objectives quite encouraging. Furthermore, we compared our methodology with the traditional HAND approach, using a single threshold height of 4 m (measured height at the outlet) for the whole catchment. The simple HAND model reached CSI=0.49, BIAS=1.55, POD=0.84 and FAR=0.46. It is clearly penalised by the overestimation (almost 50 % of false alarms), which is not surprising according to other studies (Nobre et al., 2016).
The simple structure of MHYST allows it to be used almost anywhere with few data and only two parameters. The model can, however, be used in first approximation, when a lack of time and data restrains the use of a more complex method.
For the sake of honesty, we would like to specify the theoretical limits of the MHYST approach:
The model equations were solved by using the hypothesis of a reach-scale steady uniform flow (probably one of the most simplifying assumptions one can make). This simplification is probably too extreme for highly complex situations, especially in the presence of dikes and bridges. Indeed, on the one hand, the DEM resolution is too coarse to precisely take into account hydraulic structures, and on the other hand, the DEBORD formulation is not sufficient to describe the interaction between the flow and these structures.
The DEM is a critical part of the model, because geometrical relationships and variables are directly related to the shape and distribution of elevations. Another DEM was actually tested as model input and showed much poorer results.
Moreover, since the channel geometry was unknown, hydraulic geometry equations were used to assess bankfull height and width, with fixed parameters from another study in the case of height, which may not be the optimum for this catchment, adding its share of uncertainty.
Finally, there is at this point no continuity equation between reaches, since the calculations were made for each reach separately. Uncertainties may therefore be higher in areas around connection points between reaches, especially if it is a confluence of rivers. One way to address this issue could be to add a continuity equation between the reaches, which might increase the overall coherence of the flood. However, at this point of the development of the model, we have not included this specificity.
Thus, the maps produced by MHYST should be seen as a maximum extent of the flood which can be used as a first and rapid estimation. To further test this approach, we consider that attention should first be given to the following: assessing the impact of the DEM choice, resolution and quality; testing the approach on a range of (less extreme) events and catchments, to better assess the range and stability of its parameters and performance; and improving the treatment of possible discontinuities between reaches.
The IGN DEM cannot be freely downloaded. Copernicus Emergency Management Service data and the corresponding report can be downloaded at http://emergency.copernicus.eu/mapping/list-of-components/EMSN028 (last access: 20 November 2018). French observed discharges can be downloaded at http://hydro.eaufrance.fr/indexd.php (last access: 20 November 2018).
In order to assess the sensitivity of the model to its main parameters (Kfp, Kch, α, β, ω, δ), we conducted two sensitivity analyses, using different but complementary well-known methods: Morris (Morris, 1991) and Sobol (Sobol, 2001).
A1 Morris method
The Morris method (Morris, 1991) provides a qualification of the effect a parameter can have on the outputs. It is a OAT (one-at-a-time) methodology, which means that the effect of a parameter is measured by changing its value by adding ±Δ without modifying the other parameters and by comparing the outputs. In order to provide a relevant analysis, we generated 160 sets of parameters, using the Latin hypercube sampling method, which acts as starting points from where the Morris method can assess the significance of parameters by changing their values one-at-a-time. Thus, more than 1000 simulations are needed to conduct the analysis. By using the 5 m resolution DEM we used in this paper, this study would take several days, if not weeks, to complete. But since we showed that the performance of MHYST did not really change with the resolution, we chose to use a coarser version of our DEM, which was aggregated at a 50 m resolution, by simply averaging the elevations, allowing us to complete this sensitivity analysis in only a few hours. For each permutation and for each parameter, Di, the difference in CSI divided by the computing step, is calculated. The results in terms of means and standard deviations are presented in Fig. A1. The analysis shows that the model is very sensitive to changes of ω, the exponent in the calculation of the regionalised bankfull width (Wb). The most surprising part of the analysis is the fact that Kfp has little or no effect on the model, while Kch has a moderate effect. This is contradicted by Fig. 9, which clearly shows that for a given value of Kfp, the CSI value varies only slightly for a Kch between 0.1 and 20. Kfp is, contrary to what the Morris analysis shows, a significant parameter of the model, particularly in a major overflowing event such as the one studied here, where the channel only represents a fraction of the water.
The problem might be that despite the use of a Latin hypercube sampling method, the “good” values of the parameters never meet, i.e. when ω has a sensible value, Kfp has not and vice versa. And of course, if the ω value does not coherently represent the channel, the model is not able to conduct a correct simulation (i.e. little or no flooding), leading to little or no influence of the Kfp parameter.
Moreover, the issue with sensitivity analyses such as the Morris method is that the results can be very different depending on the catchment or the event modelled. Indeed, if the water is concentrated in the channel part for a very steep catchment, a very flat one will on the contrary rely on the floodplains, and so the parameterisation of the model will add more value to Kch or Kfp. Thus, the conclusions one can make by interpreting one analysis of an example do not necessarily reflect the global behaviour of the model.
A2 Sobol method
The Sobol method (Sobol, 2001) is a variance-based sensitivity analysis which aims to compute the fraction of the variance that can be attributed to each parameter. For this study, 2×500 sets of parameters were randomly chosen with a Latin hypercube sampling method, thus creating two 500×6 matrices, XA and XB. Each column of XA has sequentially been substituted by a column of XB, corresponding to one of the six parameters, leading to six other matrices. In order to limit the computation time, the interaction of several parameters (i.e. substituting two or more columns of XA by those of XB) has not been assessed. Indeed, MHYST has been launched with the 4000 sets of parameters, with a resolution of 50 m, which takes longer than the Morris method that only needed about a thousand simulations. The first-order Sobol indices Si, which indicate the contribution of one parameter to the total variance, and the total-effect indices S−i, which calculate the total contribution of one parameter to the variance, including the possible interactions between parameters, have been computed. Then, with a bootstrap re-sampling method, the distributions of Si and S−i have been assessed, allowing several characteristics such as the bias to be computed, the standard deviation and the confidence intervals.
The results of this analysis are presented in Table A1 for Si and Table A2 for S−i. The first-order indices confirm parts of what was concluded from the Morris analysis, interpreting ω as the most influential parameter, Kch and α as moderately influential and Kfp as not influential, despite the observations we made in the article when we calibrated the parameters. The total-effect indices complete the analysis and confirm the conclusions we made with the Morris method, adding β to the list of influential parameters.
The distributions of Si and S−i show that the values calculated are not biased, but the 95 % confidence interval is rather large, which means that in some cases, the interpretation may differ. This might explain why when we set values for all downstream hydraulic geometry equations parameters (α, β, δ, ω) from regionalised studies or observations, Kfp has a greater influence which is not highlighted by the sensitivity analyses. These methodologies (Morris, Sobol) indeed explore the parameter space in detail, and even with reasonable boundaries, they can reach values that may not be consistent with the characteristics of the catchment studied. Another limitation is the fact that these analyses are only valid for this particular example (the Loing catchment and the event of May–June 2016). They should ideally be used with a larger set of catchments and events to be reliably trusted.
In order to understand why Morris and Sobol give, contrary to our initial expectation, so little importance to Kfp, we conducted a quick Sobol analysis with fixed hydraulic geometry parameters, i.e. we considered the α, β, δ and ω values used in the original study and only made Kch and Kfp vary. This time, the results confirm what we observed: and , which means that Kfp is a major parameter in our situation, and that Kch has a smaller role.
The hydraulic geometry parameters are clearly important, but if they are fixed to legitimate values estimated by observations or tables of regionalised values, their impact becomes minor in front of the Strickler coefficients.
The model presented in this paper was developed and analysed by CR during his PhD work. He also wrote the paper, which was corrected by VA and NLM.
The authors declare that they have no conflict of interest.
The first author was funded by a grant from the AXA Research Fund. Thanks are
extended to Rafal Zielinksi, who helped us access data from the Copernicus
Emergency Management Service. We would also like to thank the AXA Global P&C
research team for their advice and our discussions on the development of
simple conceptual inundation models. The MHYST model was developed using R (R
Core Team, 2015) and GFortran, Gnu compiler collection (gcc) Version
Edited by: Roger Moussa
Reviewed by: Renata Romanowicz and one anonymous referee
Afshari, S., Tavakoly, A. A., Rajib, M. A., Zheng, X., Follum, M. L., Omranian, E., and Fekete, B. M.: Comparison of new generation low-complexity flood inundation mapping tools with a hydrodynamic model, J. Hydrol., 556, 539–556, https://doi.org/10.1016/j.jhydrol.2017.11.036, 2018. a, b
Alfieri, L., Salamon, P., Bianchi, A., Neal, J., Bates, P., and Feyen, L.: Advances in pan-European flood hazard mapping, Hydrol. Process., 28, 4067–4077, https://doi.org/10.1002/hyp.9947, 2014. a, b, c
Bates, P. D., Horritt, M. S., and Fewtrell, T. J.: A simple inertial formulation of the shallow water equations for efficient two-dimensional flood inundation modelling, J. Hydrol., 387, 33–45, https://doi.org/10.1016/j.jhydrol.2010.03.027, 2010. a
Biancamaria, S., Bates, P. D., Boone, A., and Mognard, N. M.: Large-scale coupled hydrologic and hydraulic modelling of the Ob river in Siberia, J. Hydrol., 379, 136–150, https://doi.org/10.1016/j.jhydrol.2009.09.054, 2009. a
Blackburn-Lynch, W., Agouridis, C. T., and Barton, C. D.: Development of Regional Curves for Hydrologic Landscape Regions (HLR) in the Contiguous United States, J. Am. Water Resour. As., 53, 903–928, https://doi.org/10.1111/1752-1688.12540, 2017. a
CCR: Inondations de mai-juin 2016 en France – Modélisation de l'aléa et des dommages, Tech. rep., Service R&D modélisation – Direction des Réassurances & Fonds Publics, 2016. a
Falter, D., Dung, N., Vorogushyn, S., Schröter, K., Hundecha, Y., Kreibich, H., Apel, H., Theisselmann, F., and Merz, B.: Continuous, large-scale simulation model for flood risk assessments: proof-of-concept: Large-scale flood risk assessment model, J. Flood Risk Manag., 9, 3–21, https://doi.org/10.1111/jfr3.12105, 2014. a
Falter, D., Schröter, K., Dung, N. V., Vorogushyn, S., Kreibich, H., Hundecha, Y., Apel, H., and Merz, B.: Spatially coherent flood risk assessment based on long-term continuous simulation with a coupled model chain, J. Hydrol., 524, 182–193, https://doi.org/10.1016/j.jhydrol.2015.02.021, 2015. a, b
Gouldby, B., Sayers, P., Mulet-Marti, J., Hassan, M. A. A. M., and Benwell, D.: A methodology for regional-scale flood risk assessment, Water Management, 161, 169–182, https://doi.org/10.1680/wama.2008.161.3.169, 2008. a
Grimaldi, S., Petroselli, A., Alonso, G., and Nardi, F.: Flow time estimation with spatially variable hillslope velocity in ungauged basins, Adv. Water Res., 33, 1216–1223, https://doi.org/10.1016/j.advwatres.2010.06.003, 2010. a
Hunter, N. M., Horritt, M. S., Bates, P. D., Wilson, M. D., and Werner, M. G. F.: An adaptive time step solution for raster-based storage cell modelling of floodplain inundation, Adv. Water Res., 28, 975–991, https://doi.org/10.1016/j.advwatres.2005.03.007, 2005. a
Jolliffe, I. T. and Stephenson, D. B.: Forecast Verification: A Practitioner's Guide in Atmospheric Science, John Wiley & Sons, Chichester, 2003. a
Le Bihan, G., Payrastre, O., Gaume, E., Moncoulon, D., and Pons, F.: The challenge of forecasting impacts of flash floods: test of a simplified hydraulic approach and validation based on insurance claim data, Hydrol. Earth Syst. Sci., 21, 5911–5928, https://doi.org/10.5194/hess-21-5911-2017, 2017. a
Leleu, I., Tonnelier, I., Puechberty, R., Gouin, P., Viquendi, I., Cobos, L., Foray, A., Baillon, M., and Ndima, P.-O.: La refonte du système d'information national pour la gestion et la mise à disposition des données, La Houille Blanche, 1, 25–32, https://doi.org/10.1051/lhb/2014004, 2014. a, b
Leopold, L. B. and Maddock, T.: The Hydraulic Geometry of Stream Channels and Some Physiographic Implications, Tech. Rep., 252, Washington, 1953. a
McGrath, H., Bourgon, J.-F., Proulx-Bourque, J.-S., Nastev, M., and Abo El Ezz, A.: A comparison of simplified conceptual models for rapid web-based flood inundation mapping, Nat. Hazards, 93, 905–920, https://doi.org/10.1007/s11069-018-3331-y, 2018. a
Morales-Hernández, M., Petaccia, G., Brufau, P., and García-Navarro, P.: Conservative 1D–2D coupled numerical strategies applied to river flooding: The Tiber (Rome), Appl. Math. Model., 40, 2087–2105, https://doi.org/10.1016/j.apm.2015.08.016, 2016. a
Moussa, R. and Cheviron, B.: Modeling of Floods – State of the Art and Research Challenges, in: Rivers – Physical, Fluvial and Environmental Processes, edited by: Rowiński, P. and Radecki-Pawlik, A., 169–192, Springer International Publishing, Cham, https://doi.org/10.1007/978-3-319-17719-9_7, 2015. a
Neal, J., Schumann, G., and Bates, P.: A subgrid channel model for simulating river hydraulics and floodplain inundation over large and data sparse areas, Water Resour. Res., 48, W11506, https://doi.org/10.1029/2012WR012514, 2012. a, b
Nobre, A. D., Cuartas, L. A., Hodnett, M., Rennó, C. D., Rodrigues, G., Silveira, A., Waterloo, M., and Saleska, S.: Height Above the Nearest Drainage – a hydrologically relevant new terrain model, J. Hydrol., 404, 13–29, https://doi.org/10.1016/j.jhydrol.2011.03.051, 2011. a
Nobre, A. D., Cuartas, L. A., Momo, M. R., Severo, D. L., Pinheiro, A., and Nobre, C. A.: HAND contour: a new proxy predictor of inundation extent: Mapping Flood Hazard Potential Using Topography, Hydrol. Process., 30, 320–333, https://doi.org/10.1002/hyp.10581, 2016. a, b, c, d
Pons, F., Delgado, J.-L., Guero, P., and Berthier, E.: EXZECO: A GIS and DEM based method for pre-determination of flood risk related to direct runoff and flash floods, in: 9th International Conference on Hydroinformatics, 7–11 September, Tianjin, China, 2010. a
Rennó, C. D., Nobre, A. D., Cuartas, L. A., Soares, J. V., Hodnett, M. G., Tomasella, J., and Waterloo, M. J.: HAND, a new terrain descriptor using SRTM-DEM: Mapping terra-firme rainforest environments in Amazonia, Remote Sens. Environ., 112, 3469–3481, https://doi.org/10.1016/j.rse.2008.03.018, 2008. a
Schumann, G. J.-P., Neal, J. C., Voisin, N., Andreadis, K. M., Pappenberger, F., Phanthuwongpakdee, N., Hall, A. C., and Bates, P. D.: A first large-scale flood inundation forecasting model: Large-Scale Flood Inundation Forecasting, Water Resour. Res., 49, 6248–6257, https://doi.org/10.1002/wrcr.20521, 2013. a
SERTIT: EMSN-028 Flood delineation and damage assessment, France, Technical Report, 2016. a
Speckhann, G. A., Borges Chaffe, P. L., Fabris Goerl, R., Abreu, J. J. d., and Altamirano Flores, J. A.: Flood hazard mapping in Southern Brazil: a combination of flow frequency analysis and the HAND model, Hydrol. Sci. J., 63, 87–100, https://doi.org/10.1080/02626667.2017.1409896, 2018. a
Teng, J., Vaze, J., Dutta, D., and Marvanek, S.: Rapid Inundation Modelling in Large Floodplains Using LiDAR DEM, Water Resour. Manage., 29, 2619–2636, https://doi.org/10.1007/s11269-015-0960-8, 2015. a, b
Zheng, X., Tarboton, D. G., Maidment, D. R., Liu, Y. Y., and Passalacqua, P.: River channel geometry and rating curve estimation using height above the nearest drainage, J. Am. Water Resour. As., 54, 785–806, https://doi.org/10.1111/1752-1688.12661, 2018. a, b