Itzï (version 17.1): an open-source, distributed GIS model for dynamic flood simulation

Abstract

Worldwide, floods are acknowledged as one of the most destructive hazards. In human-dominated environments, their negative impacts are ascribed not only to the increase in frequency and intensity of floods but also to a strong feedback between the hydrological cycle and anthropogenic development. In order to advance a more comprehensive understanding of this complex interaction, this paper presents the development of a new open-source tool named “Itzï” that enables the 2-D numerical modelling of rainfall–runoff processes and surface flows integrated with the open-source geographic information system (GIS) software known as GRASS. Therefore, it takes advantage of the ability given by GIS environments to handle datasets with variations in both temporal and spatial resolutions. Furthermore, the presented numerical tool can handle datasets from different sources with varied spatial resolutions, facilitating the preparation and management of input and forcing data. This ability reduces the preprocessing time usually required by other models. Itzï uses a simplified form of the shallow water equations, the damped partial inertia equation, for the resolution of surface flows, and the Green–Ampt model for the infiltration. The source code is now publicly available online, along with complete documentation. The numerical model is verified against three different tests cases: firstly, a comparison with an analytic solution of the shallow water equations is introduced; secondly, a hypothetical flooding event in an urban area is implemented, where results are compared to those from an established model using a similar approach; and lastly, the reproduction of a real inundation event that occurred in the city of Kingston upon Hull, UK, in June 2007, is presented. The numerical approach proved its ability at reproducing the analytic and synthetic test cases. Moreover, simulation results of the real flood event showed its suitability at identifying areas affected by flooding, which were verified against those recorded after the event by local authorities.

Introduction

Worldwide, several records point towards an increase in the number of reported flood disaster events, which in due course have raised the magnitude of economic losses associated with their occurrence. For example, only in the last decade, 870 million people were directly affected by floods (59 000 deaths), with associated economic losses up to USD 340 billion . At the same time, the world is becoming increasingly urbanised, with more than 50 % of the global population already living in urban areas . This represents an aggravation of already existing stresses, which in combination with projected climate-induced changes will increase the expected impacts of flooding.

Indeed, population growth in urban areas is one trend that has been reported in connection with floods, which has led to more people living in potentially hazardous areas. Unless there are means to reduce the overcrowding of urban spaces located in flood-prone areas, people affected, observed damages and economic losses are set to rise further. The observed global trend of population growth in the 21st century, in particular in low-to-middle-income developing countries, opens the door to poorly planned urbanisation making their population even more vulnerable to floods . This situation endangers the future sustainability of urban environment, especially in flood-prone regions of the world. It requires adaptation strategies for urban development and drainage infrastructure to a standard that may be greater than the design level originally defined in the construction of current settlements.

Due to the observed and expected flood impacts in urban areas, it is highly necessary to develop the numerical tools for these environments that can represent the involved physical processes at an adequate level of complexity. This is the reason why urban flood modelling has recently received an increased level of attention e.g.. One of the key requirements for an adequate urban flood modelling is the ability to handle large datasets at a high spatial resolution . This involves the ability of a model to avoid both numerical instabilities resulting from the complexity of urban areas and extremely high computational times .

It is acknowledged that the natural candidate to simulate bi-dimensional surface flows in urban areas is the numerical solution of the set of non-linear shallow water equations (NLSWEs) . However, it has also been pointed out that the application of these equations in urban cases is hampered by the high computational cost resulting from a much-needed high spatial resolution and large, whole-city-scale domains . In the last 15 years, as an alternative to overcome this limitation, numerical approximations based on a diffusive wave scheme have become increasingly popular e.g.. These studies have reported that, for large domains and coarse resolution, the diffusive wave scheme is computationally more efficient than the complete solution of the NLSWEs; however, it has also been noted that at higher resolutions the numerical schemes for the diffusive systems become less efficient than those resolving the NLSWEs . Indeed, in practice, this characteristic prevents the application of diffusive wave models to urban cases.

In recognition of this limitation, presented a flood inundation model based on a partial inertia numerical scheme. In this solution, water flows are estimated by solving the inertial momentum equation using a single explicit finite-difference formulation. This modification enabled longer stable time steps and smaller computational times through the reduction of numerical operations in comparison to those required by the solution of the NLSWEs . However, this solution reported numerical instabilities when friction in the floodplain was defined with values of Manning's roughness coefficient lower than 0.03 s m${}^{-\mathrm{1}/\mathrm{3}}$, which is a value commonly found in impervious urban areas .

More recently, in order to improve the numerical stability of the local inertia solution, presented two modifications of the partial inertia momentum equation. Firstly, they introduced the use of information from both neighbouring cell interfaces to form a one-dimensional three-point stencil; secondly, they resolved the friction term in two dimensions, making use of a five-point stencil to overcome the limitations of the staggered-grid data model. These two modifications have been shown to have a better numerical stability than the original local inertia scheme .

Despite these advances, to the best of the author's knowledge, there is no open-source model based on this approximation that integrates seamlessly with geographic information system (GIS) software. It should be noted that CAESAR-LISFLOOD  is an open-source numerical tool, which applies a similar approach but using the partial inertia scheme described by . Meanwhile, the presented implementation uses the damped partial inertia form of the equation, described by and , and allows a straightforward integration within a GIS environment.

Hence, this investigation presents an open-source implementation of the latest advances of the damped partial inertia approximation, under the general public license (GPL). The numerical model Itzï is written in Python, and its aim is to simulate surface flows induced by intense rainfall in urban settings . The model is tightly integrated with the open-source geographical information system known as GRASS , which allows the easy use of space–time varying raster datasets both as inputs and outputs. Moreover, it enables the automatic integration of geographic datasets from sources that have different spatial resolutions (e.g. elevation, land use, soil types) to adequately describe floodplain topographies.

The model verification will be carried out through the comparison of numerical results against three different test cases. Firstly, results will be compared against two analytical solutions of the non-linear shallow water equations. Secondly, to ensure that the model is able to adequately predict inundation flows, a reproduction of a widely accepted standard benchmark test case (no. 8a), published in a report from the UK Environment Agency , will be implemented. In this case, numerical results will be compared against those obtained from a well-established model based on the same equations . Finally, in order to test the ability of Itzï to reproduce a real flood event, the 2007 flood registered in the city of Kingston upon Hull, UK, will be presented. Affected areas identified by the numerical tool will be compared against those surveyed by local authorities.

This paper is organised as follows: Sect.  introduces the numerical scheme used for the solution of the partial inertia shallow water equations; Sect.  describes the numerical implementation in an open-source GIS platform, while numerical results for the test cases are illustrated in Sect. .

Numerical scheme

The developed programme uses an explicit finite-difference scheme to solve the simplified partial inertia shallow water equations described by  and . Figure  illustrates the variables used by the scheme in the $x$ dimension and their variations in time. On the other hand, Fig.  introduces a complete 2-D view of the same staggered grid and variables utilised in the numerical scheme. As shown in Fig. , water surface elevation $y$ and water depth $h$ are evaluated in the centre of cells, while the water flow $q$ (or velocity $u$) variables are evaluated at the cell interfaces.

Grid and variables used in the numerical scheme. Time is shown in the vertical axis. Variables at the bottom are used at time $t$ to calculate variables at the top at time $t+\mathrm{\Delta }t$.

[Figure omitted. See PDF]

The mass flux (e.g. water flow) is obtained by solving the 1-D simplified momentum equation at interfaces between cells using the value of $q$ at these interfaces (rather than in the centre of cells). To provide a bi-dimensional representation of the flow, momentum itself is updated at the cell interfaces with an explicit discretisation of the momentum equation in each direction separately. The numerical method is simple and extremely efficient from a computational point of view. For simplicity, in this section, we will present only the flow equation for the $x$ dimension. The exact same principle applies for water flows in the other direction, which is represented by the $y$ dimension.

Adaptive time step

In a similar vein to previous developments, an adaptive time-stepping method is used to estimate the suitable model time step based on the standard Courant–Friedrichs–Lewy (CFL) condition. The time step $\mathrm{\Delta }t$ is calculated at each time step by means of Eq. (). $$\mathrm{\Delta }t=\mathit{\alpha }{\displaystyle \frac{min\mathit{\{}\mathrm{\Delta }x,\mathrm{\Delta }y\mathit{\}}}{\sqrt{g\times h_{max}}}},$$ where $h_{max}$ is the maximum water depth within the domain, $g$ the acceleration due to the gravity and $\mathit{\alpha }$ an adjustment factor because the CFL condition is necessary but not sufficient to ensure stability. propose a value of $\mathit{\alpha }=\mathrm{0.7}$ as a default, as this has been shown to allow the appropriate simulation of subcritical flooding conditions. When $h_{max}$ tends to $\mathrm{0}$, $\mathrm{\Delta }t$ is set to a user-defined time step $\mathrm{\Delta }{t}_{max}$, which represents the maximum value for $\mathrm{\Delta }t$. Here, the default for this value has been set to 5 s. It could be adjusted by the user to optimise computation time while preserving numerical stability.

Flow calculation

The flow at each cell interface is calculated with Eq. (). $$q_{i+\mathrm{1}/\mathrm{2}}^{t+\mathrm{\Delta }t}={\displaystyle \frac{\left(\mathit{\theta },q_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}},+,(,\mathrm{1},-,\mathit{\theta },),{\displaystyle \frac{q_{i-\mathrm{1}/\mathrm{2}}^{\mathrm{t}}+q_{i+\mathrm{3}/\mathrm{2}}^{\mathrm{t}}}{\mathrm{2}}}\right)+g{h}_{\mathrm{f}}\mathrm{\Delta }tS}{\mathrm{1}+g\mathrm{\Delta }t{n}^{\mathrm{2}}\left|\right|q_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}\left|\right|/{h_{\mathrm{f}}}^{\mathrm{7}/\mathrm{3}}}},$$ where subscripts $i$ and $t$ denote space and time indices (see Fig. ).

The flow depth $h_{\mathrm{f}}$ is the difference between the highest water surface elevation $y$ and the highest terrain elevation $z$. It is calculated at the cell face using Eq. (). This value is used as an approximation of the hydraulic radius.

$\mathit{\theta }$ is a coefficient defining the importance taken by the average of upstream and downstream flows over the flow at the considered cell face ($q_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}$). proposes to set this weighting factor to 0.9. If $\mathit{\theta }$ is set to 1, neighbouring flows are not taken into account, being equivalent to the formula proposed by . In some rare cases, especially when $\mathit{\theta }$ is low, the flow term could end up with a different sign to the slope term. When this happens, the weighting scheme is dropped and the numerator of the equation becomes equal to the formulation presented by .

The slope $S$ is calculated using Eq. (). With the flow being calculated at cell interfaces, Manning's $n$ is obtained by averaging the neighbouring values, as shown in Eq. (). $$\begin{array}{}& {\displaystyle }& {\displaystyle }{h_{\mathrm{f}}}_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}=max\left\{y_i^{\mathrm{t}},,,y_{i+\mathrm{1}}^{\mathrm{t}}\right\}-max\left\{z_i^{\mathrm{t}},,,z_{i+\mathrm{1}}^{\mathrm{t}}\right\}& {\displaystyle }& {\displaystyle }S={\displaystyle \frac{y_i^{\mathrm{t}}-y_{i+\mathrm{1}}^{\mathrm{t}}}{\mathrm{\Delta }x}}& {\displaystyle }& {\displaystyle }{n}_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}={\displaystyle \frac{\left(n_i^{\mathrm{t}},+,n_{i+\mathrm{1}}^{\mathrm{t}}\right)}{\mathrm{2}}}\end{array}$$

The vector norm $\left|\right|q_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}\left|\right|$ is calculated using Eq. () given by . $$\left|\right|q_{i+\mathrm{1}/\mathrm{2}}^{\mathrm{t}}\left|\right|=\sqrt{{\left(q_{y,i+\mathrm{1}/\mathrm{2},j}^{\mathrm{t}}\right)}^{\mathrm{2}}+{\left(q_{x,i+\mathrm{1}/\mathrm{2},j}^{\mathrm{t}}\right)}^{\mathrm{2}}}$$

A 2-D view of the staggered grid used by the numerical scheme. The variable in bold will be calculated at $t+\mathrm{\Delta }t$ using all the displayed variable values at time $t$.

[Figure omitted. See PDF]

Inconveniently, due to the use of a staggered grid, $q_{y,i+\mathrm{1}/\mathrm{2},j}^{\mathrm{t}}$ is not being calculated by the model. To overcome this, the values of the neighbouring cells are used instead as shown in Eq. (). The positions of the given points are shown in Fig. . $$q_{y,i+\mathrm{1}/\mathrm{2},j}^{\mathrm{t}}={\displaystyle \frac{q_{i,j-\mathrm{1}/\mathrm{2}}+q_{i,j+\mathrm{1}/\mathrm{2}}+q_{i+\mathrm{1},j-\mathrm{1}/\mathrm{2}}+q_{i+\mathrm{1},j+\mathrm{1}/\mathrm{2}}}{\mathrm{4}}}$$

Water depth calculation

The new water depth at each cell is calculated using Eq. (). It consists of the sum of the current depth $h^{\mathrm{t}}$, the external source terms (rainfall, user-defined flow, etc.) $h_{\mathrm{ext}}^{\mathrm{t}}$ and the flows passing through the four faces of each cell. If the new calculated water depth is negative, it is set to 0 and the additional volume is accounted for. $$h^{t+\mathrm{\Delta }t}=h^{\mathrm{t}}+h_{\mathrm{ext}}^{\mathrm{t}}+{\displaystyle \frac{{\sum }^{\mathrm{4}}{Q}_{i,j}^{\mathrm{t}}}{\mathrm{\Delta }x\mathrm{\Delta }y}}\times \mathrm{\Delta }t$$

Rain routing

In order to maintain stability during events with direct rainfall, a rain-routing mechanism is implemented using a simple method described by . It consists of applying a constant velocity to the flow when water depth is below a user-given threshold. Before the simulation begins, the software calculates the draining direction of each cell of the domain. The drainage direction is determined by the highest slope out of the four neighbouring cells. During the simulation, the routing scheme is applied at each cell interface when each of the following conditions are true:

• $h_{\mathrm{f}}<h_{\mathrm{fmin}}$,

• the considered direction is allowed for routing according to the routing map and

• the slope $S$ is in the same direction as the above routing direction.

The routing flow is then calculated using a constant user-given velocity. According to , a depth threshold of 5 $\mathrm{mm}$ and a routing velocity of 0.1 $\mathrm{m}$ ${\mathrm{s}}^{-\mathrm{1}}$ gives good results.

Infiltration

In Itzï, the infiltration could be represented by a map or a time series of maps containing a fixed value for the infiltration rate in $\mathrm{mm}{\mathrm{h}}^{-\mathrm{1}}$. Alternatively, the Green–Ampt method can be used, as shown in Eq. (), where $f$ is the infiltration rate ($\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}$), $K$ the hydraulic conductivity ($\mathrm{L}{\mathrm{T}}^{-\mathrm{1}}$), ${\mathit{\theta }}_{\mathrm{e}}$ the effective porosity ($\mathrm{L}{\mathrm{L}}^{-\mathrm{1}}$), $\mathit{\theta }$ the initial water soil content ($\mathrm{L}{\mathrm{L}}^{-\mathrm{1}}$), ${\mathit{\psi }}_{\mathrm{f}}$ the wetting front capillary pressure head ($\mathrm{L}$) and $F$ the infiltration amount ($\mathrm{L}$). $$f=K\left(\mathrm{1},+,{\displaystyle \frac{({\mathit{\theta }}_{\mathrm{e}}-\mathit{\theta }){\mathit{\psi }}_{\mathrm{f}}}{F}}\right)$$

Implementation in Python

The software is written in the Python programming language and integrates tightly with the open-source GIS GRASS . It employs the libraries PyGRASS  to access the geographical functions and TGRASS  for the temporal management of both the input and the output data. Additionally, further optimisation of the numerical code was carried out by means of a Python profiler that recorded the call stack of the executing code, thus accounting for the time spent in the solution of each function within the code. This enabled the parallelisation through Cython  of those functions with the highest computational cost, reducing the overall computing time by taking advantage of the multicore CPU. The integration of this numerical model within GRASS provides Itzï with the following relevant characteristics:

• The spatiotemporal data management is straightforward as the integration within a GIS platform reduces the time spent on preparation of entry data and the analysis of results. Modifying the spatial extent and resolution of the simulation is done by simply changing the GRASS computational region, without the need for changing the entry data.

• Forcings could be of heterogeneous resolutions. For example, elevation at 5 m, rainfall at 1 km and friction coefficient at 30 m. Itzï will automatically read the data at the resolution defined by the computational region and uses the data seamlessly, without user intervention.

• Input data can vary in space and time (i.e. raster time series), permitting the use of, for example, spatially distributed rainfall or time-varying friction coefficients.

• The ability to use absolute time references in the form of date and time for start and end of the simulation facilitates the usage of historical rainfall. It is therefore possible to have several years of rainfall data stored in the GIS and to simulate just one specific event, without further data preprocessing.

Itzï is operated by a command-line interface taking a parameter file as input. If several input files are given, they are run in batch mode. The user can ask the software to output the following raster time series:

• water depth ($h$) and surface elevation ($h+z$);

• flow velocity magnitude and direction;

• volumetric flows in $x$ and $y$ directions;

• average volume added or subtracted to the domain by the action of infiltration, rainfall, user-defined inflow, drainage capacity or the application of boundary conditions; and

• volume created due to numerical instability (see Sect. ).

Additionally, the software can produce a CSV file that summarises the statistics mentioned above.

Verification and evaluation

Analytic test cases

For the analytic test cases, we utilise numerical experiments aimed at testing subcritical flow simulations recently published in a compilation of shallow water analytic solutions for hydraulic and environmental studies  . Both cases described here are constituted by a 1 $\mathrm{km}$ long channel of the MacDonald type , discretised at 5 $\mathrm{m}$ resolution. These test cases were generated with the free software SWASHES available at https://sourcesup.renater.fr/projects/swashes/.

The first case corresponds to a constant upstream flow of 2 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, while the second one combines an upstream flow of 1 ${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ and a uniform rainfall with an intensity of 0.001 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. In the model, the input flow is given as a mass addition. This creates an artificially high water level at the most upstream cell, where the input flow is added. Given that the goal of the analytic tests is to verify the validity of the numerical scheme, we determine the root mean squared error (RMSE) omitting the very first cell of the domain. Figures  and illustrate the performance of Itzï at reproducing results from the analytic solution, reporting RMSEs of 0.002 and 0.03 m, for each case, respectively. Those RMSE values are 1 to 2 orders of magnitude lower than the vertical accuracy of airborne lidar , demonstrating the suitability of the implemented simplified scheme to simulate subcritical flow conditions.

Comparing Itzï with an analytic solution (the case of the one-dimensional MacDonald long channel).

[Figure omitted. See PDF]

Comparing Itzï with an analytic solution (the case of the one-dimensional MacDonald long channel with rain).

[Figure omitted. See PDF]

Direct rainfall and sewer overflow in an urban setting

In order to further test the numerical model, previously published benchmark test cases for 2-D flood inundation modelling tools by the UK Environment Agency  were implemented. These cases correspond to a benchmarking exercise assessing the latest generation of 2-D hydraulic modelling tools for a variety of purposes in Flood and Coastal Risk Management (FCRM) to support Environment Agency decision making. This dataset is available online at http://evidence.environment-agency.gov.uk/FCERM/Libraries/FCERM_Project_Documents/Benchmarking_Model_Data.sflb.ashx.

In particular, one hypothetical test case was utilised to verify the proper implementation of the numerical scheme to simulate physical processes controlling flood movement across a floodplain. The test case (test case number 8a in the Environment Agency study) corresponds to a synthetic event which does not relate to any real event . The modelled area is in the city of Glasgow, UK (Cockenzie Street and surrounding streets), and is approximately 400 by 960 m. Ground elevations span a range of 21 to 37 m. The flood is assumed to arise from two sources: a uniformly distributed rainfall (applied only to the modelled area) and a point inflow representing a sewer outflow from a surcharging culvert. For completeness, Fig.  shows both forcings described by the hyetograph and the hydrograph specified at the point inflow. The digital elevation model (DEM) has a spatial resolution of 0.5 m, which is resampled to 2 m resolution for modelling purposes. This represents the terrain model with no vegetation or buildings and was created from a lidar dataset provided by the UK Environment Agency. The roughness coefficient was determined following the classification of the area with two land cover roughnesses: roads and pavements ($n=\mathrm{0.02}$) and everywhere else ($n=\mathrm{0.05}$). The model was run until time $t=\mathrm{83}$ min, as this was considered enough to allow the flood to pond in the lower parts of the modelled domain.

Rectangular hyetograph and point inflow hydrograph for the EA test 8a.

[Figure omitted. See PDF]

Numerical results obtained with Itzï have been compared to those obtained with the implementation of the acceleration solver from LISFLOOD-FP . This is done as the latter is considered the reference implementation of the numerical scheme here employed. For this comparison, eight different locations within the numerical domain were selected to compare water depths estimated by both numerical tools. Figure  illustrates the utilised digital elevation model, along with the position of the inflow point and selected control points for the comparison of model results.

DEM of EA test 8a showing the numbered control points (crosses) and inflow point (triangle).

[Figure omitted. See PDF]

Simulation parameters for the EA test case 8a.

Parameter	Value
$\mathrm{\Delta }{t}_{max}$	5 $\mathrm{s}$
$\mathit{\alpha }$	$\mathrm{0.7}$
$\mathit{\theta }$	$\mathrm{0.7}$
$v_{\mathrm{r}}$	0.1 $\mathrm{m}$ ${\mathrm{s}}^{-\mathrm{1}}$
$h_{\text{fmin}}$	0.005 $\mathrm{m}$

$v_{\mathrm{r}}$: routing velocity

The simulation is run for 83 min with both LISFLOOD-FP and Itzï using the same parameters, shown in Table . Figure  illustrates the time series of water level produced by both numerical models at the eight selected locations, as well as the time series of differences between the results of the two models. It is shown that in all eight selected locations, numerical results from both models are similar with small differences identified at the arrival time of the flood wave in each location. These differences are ascribed to the way that LISFLOOD-FP handles entry data in comparison to Itzï. In the first case, a temporal interpolation is performed during the simulation at each time step, while in the second case, this process should be carried out during the preparation of input space–time raster datasets. The RMSEs at the eight locations range from 0.2 to 10.6 mm (see Fig. ). This indicates that the numerical solution of the partial inertia approximation implemented in Itzï generates results with the same level of skill as the reference model.

Comparison of water depths at control points. The blue line represents the differences of water depth between Itzï and LISFLOOD-FP.

[Figure omitted. See PDF]

Real flood event: Kingston upon Hull, UK

In order to evaluate the ability of the presented tool to reproduce a real-life flood, we use an event that occurred on 25 June 2007 in the city of Kingston upon Hull, UK. Widespread flooding was registered due to an extreme precipitation event (hp $=$ 110 mm) that lasted for more than 12 h. The estimated 24 h rainfall was estimated between a $\mathrm{1}/\mathrm{150}$ and $\mathrm{1}/\mathrm{200}$ year return period . As reported by , the major forcing of this event was surface water runoff both locally in the urban area and through the rural lands surrounding the city.

Figure  introduces the location of the study area within Great Britain. Figure  presents the utilised DEM that was obtained for the area by means of a lidar system and has a spatial resolution of 5 m. The study area measures 87.7 ${\mathrm{km}}^{\mathrm{2}}$ and the numerical domain is comprised of 3.5 million cells. It can be seen that most of the city is located in a low-lying region, with some parts that are at or even below sea level.

Location of the Hull study area (satellite imagery Copernicus Sentinel 2016).

[Figure omitted. See PDF]

The numerical forcing condition was defined with measured rainfall recorded by a pluviometer located at the University of Hull, with the assumption that rainfall was uniform over the whole numerical domain. Figure  introduces the recorded hyetograph, which has a temporal resolution of 1 h. The whole event lasts for 24 h.

Digital elevation model of Hull.

[Figure omitted. See PDF]

Hyetograph at the University of Hull.

[Figure omitted. See PDF]

Multiple observation datasets for this flood event were obtained from the UK Environment Agency and Hull City Council. They have been collected by aerial photography and by carrying out a poll among the residents, respectively . Water depths in some local places were reported to be up to 3 m, but for most affected areas the depth was observed as being less than 1 m . Flood extent maps were produced by each of these authorities and are shown in Fig. . It could be noted that the two zones classified by the two administrations show significant differences highlighted in Table . Notably, less than half of each observation could be validated by the other one. Furthermore, due to the limitations of the data collection methods, it is highly possible that some areas that were actually affected might not be identified .

Reliable results of flood-routing and inundation simulation rely on an accurate estimation of the resistance coefficient. The hydraulic resistance of open-channel and overland flow results are typically represented through Manning's friction coefficient. It is acknowledged that the level of detail at which this process can be represented is dependent on the scale of the simulation. Usually, when modelling water flows in river floodplains, this drag may be conceptually divided into several zones, namely, main channel, soil grain roughness and vegetative roughness e.g.. Therefore, it is important to represent the spatial variability of this parameter. For this, data from the global land cover product (with resolution of 30 m – GLC30) provided by National Geomatics Center of China are employed. Figure  shows the repartition of those land cover classes in the study area. The land cover classes are then related to Manning's $n$ using values proposed by , as shown in Table .

Observed flood extent in Hull. Green indicates the UK Environment Agency. Orange indicates the Hull City Council. Purple indicates both the EA and City Council.

[Figure omitted. See PDF]

Comparison of identified flood extents.

Collecting entity	Area (km${}^{\mathrm{2}}$)
Environment Agency	5.16
City Council	6.18
Both City Council and EA	2.33

Relationship between land cover classes and Manning's $n$ in Hull.

Land cover class	Manning's $n$ (s m${}^{-\mathrm{1}/\mathrm{3}}$)
Urban area	0.019
Cultivated land	0.04
Forest	0.05
Grassland	0.03

According to , the drainage capacity in the urban and rural areas could be estimated to 70 and 15 mm per day, respectively. Therefore, a drainage capacity map has been generated using the land cover map (see Fig. ), where the artificial surfaces have been assigned a constant value of 2.917 mm h${}^{-\mathrm{1}}$ and the remaining areas a value of 0.625 mm h${}^{-\mathrm{1}}$.

In this presented case, the ability of Itzï to employ entry data with heterogeneous resolutions is advantageous. Although the simulation is run at 5 m, all the maps are kept at their native resolutions, which are 5 m for the DEM and 30 m for rasters deducted from land cover. The integration of the numerical tool within the GRASS GIS environment automatically provides Itzï with all the maps at the selected modelling resolution. This process does not require any preprocessing of the data by the user, which is seen as an advantage of Itzï over other numerical tools. This ability saves time during the preparation of entry data for modelling purposes.

Figure  shows five snapshots of the development of inundation in the floodplain (panels a to e), along with the maximum depth attained during the whole simulation time (panel f) during the simulation. Indeed, the result shown in panel (f) of this figure can be used for comparison against the identified flooded areas reported in Fig. .

Map of the land cover classes in Hull from GLC30 .

[Figure omitted. See PDF]

Water depth in Hull. Panels (a) to (e) show the evolution of water depth (m). Panel (f) shows the maximum water depth (m).

[Figure omitted. See PDF]

Since in this case, the map of the measured flood extent is dichotomous, i.e. the area is either inundated or not, numerical results will be evaluated in this way too. To achieve this, a wet/dry condition was employed and defined through a threshold water depth over which the numerical cell is considered flooded. Here, the maximum computed water depth has been compared with the union of the flooded areas presented in Fig.  using different water depth thresholds. Common dichotomous skill scores have been used for this task. A succinct description of those scores is provided in Table . Among them, the critical success index is commonly used by the community of hydraulic modelling as an indication of the fit .

Short description of the skill scores used in this paper .

Name	Description	Perfect score
Equitable threat score	Adequation between the observed flooded area and the simulated flooded areas, taking into account the number of hits due to chance	1.0
Heidke skill score	Informs about the accuracy of the model compared to random chance	1.0
Probability of detection	Fraction of flooded areas that were correctly modelled	1.0
Bias score	Informs how the frequency of computed flooded areas compares to the frequency of observed flooded areas	1.0
Success ratio	Fraction of the modelled flooded areas that were actually observed	1.0
Odds ratio skill score	Indicates how much the model improves the prediction over random prediction	1.0
Probability of false detection	Fraction of the observed non-flooded area that was modelled as inundated	0.0
Critical success index	Indicates how well the computed flooded areas correspond to the observed inundated areas	1.0
False alarm ratio	Fraction of the computed flooded area that actually did not inundate	0.0
Accuracy	Fraction of the computed flooded and non-flooded areas that were correctly predicted	1.0

Table  shows the different scores used to compare the computed area and the measured one. It could be noted that when the depth threshold is increasing (and therefore the computed flooded area decreases), the probability of false detection (PFD), probability of detection (POD), bias score and false alarm ratio (FAR) are decreasing; meanwhile the accuracy and success ratios are increasing. This could be explained by the predominance of the correct negative area that increases together with the threshold. The scores that take this bias into account, like equitable threat score, Heidke skill score and critical success index, are higher when the threshold depth is 20 cm.

Verification of the flood prediction in Hull.

	Threshold depth (m)
Score	0.05	0.10	0.20	0.30
Equitable threat score	0.209	0.242	0.277	0.258
Heidke skill score	0.346	0.390	0.434	0.410
Probability of detection	0.722	0.655	0.539	0.414
Bias score	1.958	1.557	1.038	0.690
Success ratio	0.369	0.420	0.520	0.599
Odds ratio skill score	0.780	0.795	0.830	0.849
Probability of false detection	0.243	0.178	0.098	0.054
Critical success index	0.323	0.344	0.360	0.324
False alarm ratio	0.631	0.580	0.480	0.401
Accuracy	0.751	0.795	0.842	0.858

The reproduction of this real flood event, with the considered assumptions, demonstrates the suitability of Itzï to simulate urban inundations. However, it is acknowledged that for this case there are several sources of uncertainty that may prevent a better skill in the numerical results. These are

• uncertainty in the observed flood extent maps, as shown by differences in the areas identified by the UK Environment Agency and the Hull City Council (see Fig. );

• a simplistic representation of the drainage system;

• the utilisation of uniform rainfall in the numerical domain;

• the lack of consideration for the infiltration processes; and

• an inadequate consideration of some terrain features, like walls and buildings, that might not be well represented in a 5 m DEM.

Lastly, we can evaluate the software performance of Itzï for this test case in terms of computational cost and numerical stability. The results of this analysis are presented in Table . It shows that the run time is about 3 h using an Intel^® Core^™ i7-4790 processor. Apart from using a faster computer with more cores, possible paths to improve running times include further optimisation of the code, making use of CPU vector extensions (like AVX and similar) and reprogramming the computationally intensive parts of the software to run on graphical processing units (GPUs). Additionally, it could be seen in Table  that the volume created due to computational instability (see Sect. ) represents only 0.03 % of the water volume in the domain at the end of the simulation, denoting a rather stable numerical scheme.

Evaluation of software performance of Itzï in the test case of Hull. $V_{\mathrm{error}}$ is the total volume created due to numerical instability during the simulation. $V_{\mathrm{final}}$ is the total volume in the domain at the end of the simulation.

Variables	Values
Computational time	03:06
$V_{\mathrm{error}}$ (m${}^{\mathrm{3}})$	1030
$V_{\mathrm{final}}$ (m${}^{\mathrm{3}})$	4 027 749
$V_{\mathrm{error}}/V_{\mathrm{final}}$	0.03 %

Conclusions

This paper presented Itzï, a new dynamic hydrologic and hydraulic model made for simulation of surface flows in two dimensions. The numerical tool is written in Python and uses a partial inertia numerical scheme coupled with a simple rainfall-routing system. Notably, the tool is tightly coupled with the open-source GIS GRASS, which allows easy management of input and output data of heterogeneous resolution and the use of data varying in space and time in the form of raster time series, like precipitation or friction coefficients.

The validity of the numerical scheme has been demonstrated by two analytic benchmarks, resulting in RMSEs 1 to 2 orders of magnitude lower than the airborne lidar vertical accuracy . The implementation has also been compared to a reference implementation of the same numerical scheme (LISFLOOD-FP) for the test case number 8a of the UK Environment Agency two-dimensional hydraulic model benchmark . In that test, simulating a urban flood from direct rainfall and sewer overflow, Itzï gives results nearly equal to LISFLOOD-FP. The presented model has then been evaluated with the numerical reproduction of a flood event that occurred in June 2007 in the city of Hull, UK, and showed Itzï's ability to identify the main flooded areas.

Future work will focus on (1) applying the advantages of the GIS-integrated model by using global datasets for flood modelling and (2) coupling this overland flow model with a drainage network model. The latter will enhance Itzï's capability of performing more complete hydrologic and hydraulic simulations in the urban environment.

The source code of Itzï is freely available under the GNU GPL . The link to the online repository and a user manual are accessible on the project website (http://www.itzi.org).

The authors declare that they have no conflict of interest.

Acknowledgements

We warmly acknowledge Dapeng Yu from Loughborough University, UK, for kindly providing the DEM, map of flooded areas and rainfall data for the numerical reproduction of the flood event in Kingston upon Hull. Laurent Courty is grateful for the support given by a full PhD scholarship from UNAM's Coordination of Postgraduate Studies. He extends his warm thanks to Jeffrey C. Neal from the University of Bristol for casting light on the numerical scheme implementation in LISFLOOD-FP and Mark Trigg from the University of Leeds for his insights on testing the early versions of Itzï.Edited by: B. Jackson Reviewed by: two anonymous referees