Interactive comment on “ Itzï ( version 16 . 8 ) : An open-source , distributed GIS model for dynamic flood simulation

Explications of the volume error calculation has been added to the numerical scheme description. To make clearer the values exported by the model, we added a listing of the map time-series that can be output by Itzï. Those include the volume error. However, the version 16.8 of the software described in the first submitted manuscript was suffering from a bug in the volume balance reporting. This bug has been corrected in the version 17.1, which has been published in January this year. Therefore, we decided to change the described version to the last published one, the 17.1. All the results presented in the paper have been updated to reflect the new software version. Furthermore, instead of presenting two simulations that give almost equal results for


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 to 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 340 billion US Dollars (IFRC, 2015).At the same time, the world is becoming increasingly urbanised, with more than 50 percent of the global population already living in urban areas (Zevenbergen et al., 2010).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 (United Nations, 2014).This situation endangers the future sustainability of urban environments, 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 this 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.Hsu et al., 2000;Mark et al., 2004;Schmitt et al., 2004;Sampson et al., 2012;Yu and Coulthard, 2015).One of the key requirements for an adequate urban flood modelling is the ability to handle large datasets at a high spatial resolution (Yu and Lane, 2006;Fewtrell et al., 2011).This involves the ability of a model to avoid both numerical instabilities resulting from the complexity of urban areas and extremely high computational times (Chang et al., 2015).
It is acknowledged that the natural candidate to simulate bidimensional surface flows in urban areas is through the numerical solution of the set of non-linear shallow water equations (NLSWE) (Hunter et al., 2008).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 (Neal et al., 2012b).In the last 15 years, and as an alternative to overcome this limitation, numerical approximations based on a diffusive wave scheme have become increasingly popular (e.g., Bates and De Roo, 2000;Bradbrook et al., 2004;Chen et al., 2005;Yu and Lane, 2006;Leandro et al., 2014;Guidolin et al., 2016).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 NLSWE; however, it has also been noted that at higher resolutions the numerical schemes for the diffusive systems become less efficient than those resolving the NLSWE (Hunter et al., 2008).Indeed, in practice this characteristic prevents the application of diffusive wave models to urban cases.
In recognition of this limitation, Bates et al. (2010) presented an : 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 NLSWE (Neal et al., 2012a).However, this solution reported numerical instabilities when friction in the floodplain was defined with low values of the Manning's roughness coefficient :::: lower :::: than :::: 0.03 :::::::: s.m −1/3 , which is an intrinsic characteristic of : a ::::: value ::::::::: commonly ::::: found :: in : impervious urban areas :::::::::::: (Chow, 1959).
More recently, in order to improve the numerical stability of the local inertia solution, De Almeida et al. (2012) 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; and 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 showed to have a better numerical stability than the original local inertia scheme (De Almeida and Bates, 2013).
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 a GIS software.It should be noted that CAESAR-LISFLOOD (Coulthard et al., 2013) is an opensource numerical tool, which applies a similar approach but using the partial inertia scheme described by Bates et al. (2010).
Meanwhile, the presented implementation uses the damped partial inertia form of the equation described by De Almeida et al.
(2012) and De Almeida and Bates (2013) and allows a straight forward 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 GPL license.The numerical model called Itzï is written in Python and its aim is to simulate surface flows induced by intense rainfall in urban settings (Courty and Pedrozo-Acuña, 2016a, b).The model is tightly integrated with the open-source Geographical Information System known as GRASS (Neteler et al., 2012), which allows the easy use of spacetime 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 adequatly predict inundation flows, a reproduction of a widely accepted standard benchmark test case (no.8a), published in a report from the UK Environment Agency (Néelz and Pender, 2013) will be implemented.In this case, numerical results will be compared against those obtained from a well established model based on the same equations (Bates et al., 2013).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. 2 introduces the numerical scheme used for the solution of the partial inertia shallow water equations; Sect.2.6 describes the numerical implementation in an open-source GIS platform, while numerical results for the test cases are illustrated in Sect.3.
staggered grid and variables utilised in the numerical scheme.As shown in Fig. 1, 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.
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 bidimensional representation of the flow, momentum itself is updated at the cell interfaces with an explicit discretization 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 ∆t is calculated at each time-step by means of equation Eq. ( 1).
Where h max is the maximum water depth within the domain, g the acceleration due to the gravity and α an adjustment factor because the CFL condition is necessary but not sufficient to ensure stability.Following the proposed value by De Almeida et al.
(2012), α = 0.7 as a default value, as this has been used in an appropriate manner to simulate subcritical flooding conditions.
When h max tends to 0, ∆t is set to a user-defined time-step ∆t max , which represents the maximum value for ∆t.Here the default for this value has been set to 5 seconds.It could be adjusted by the user to optimize computation time while preserving numerical stability.

Flow calculation
The flow at each cell interface is calculated with Eq. ( 2).
where subscripts i and t denotes space and time indices (Cf.Fig. 1).
The flow depth h 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. ( 3).This value is used as an approximation of the hydraulic radius.
θ is a coefficient defining the importance taken by the average of upstream and downstream flows over the flow at the considered cell face (q t i+1/2 ).De Almeida et al. (2012) proposes to set this weighting factor to 0.9.If θ is set to 1, neighbouring flows are not taken into account, being equivalent to the formula proposed by Bates et al. (2010).In some rare cases, especially when θ 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 Bates et al. (2010).
The slope S is calculated using Eq. ( 4).Flow being calculated at cell interfaces, n is obtained by averaging the neighbouring values, as shown in Eq. ( 5).
Inconveniently, due to the use of a staggered grid, q t y,i+1/2,j is not being calculated by the model.To overcome this, the value of the neighbouring cells are used instead as shown in Eq. ( 7).The positions of the given points are showed on Fig. 2.

Water depth calculation
The new water depth at each cell is calculated using Eq. ( 8).It consists of the sum of the current depth h t , the external source terms (rainfall, user-defined flow etc.) h t ext and the flows passing through the four faces of each cell.

Rain routing
In order to maintain stability during events with direct rainfall, a rain routing mechanism is implemented using a simple method described by Sampson et al. (2013).It consist :::::: 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: the considered direction is allowed for routing according to the routing map, 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 Sampson et al. (2013), a depth threshold of 5mm and a routing velocity of 0.1m.s−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 mm.h −1 .Alternatively, the Green-Ampt method can be used, as shown in Eq. ( 9), where f is the infiltration rate (L.T −1 ), K the hydraulic conductivity (L.T −1 ), θ e the effective porosity in (L.L −1 ), θ the initial water soil content (L.L −1 ), ψ f the wetting front capillary pressure head (L) and F the infiltration amount (L).

Integration in Python
The software is written in the Python programming language, is operated by a command line interface and integrates tightly with the open-source GIS GRASS (Neteler et al., 2012).It employs the libraries PyGRASS (Zambelli et al., 2013) to access the geographical functions and TGRASS (Gebbert and Pebesma, 2014) 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 records 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 (Behnel et al., 2011) of those functions with the highest computational cost, reducing the overall computing time by taking advantage of multi-cores CPU.The integration of this numerical model within GRASS provides Itzï with the following relevant characteristics: -The spatio-temporal 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 5m, rainfall at 1km, friction coefficient at 30m etc. 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 form of date and time for start and end of the simulation, facilitating 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 pre-processing.
3 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 (Delestre et al., 2013).Both cases described here are constituted by a 1km long channel of MacDonald's type (MacDonald et al., 1997), discretized at 5m 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 2m 2 .s−1 , while the second one combines an upstream flow of 1m 2 .s−1 and a uniform rainfall with an intensity of 0.001m.s−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 3 and 4 illustrate the performance of Itzï at reproducing results from the analytic solution, reporting RMSE of 0.003 and 0.03 metres, for each case respectively.Those RMSE values are one to two orders of magnitudes lower than the vertical accuracy of airborne LiDAR (Hodgson and Bresnahan, 2004), demonstrating the suitability of the implemented simplified scheme to simulate subcritical flow conditions.

Direct rainfall and sewer overflow in an urban setting
In order to further test the numerical model, previously published benchmark test cases for 2D flood inundation modelling tools by the UK Environment Agency (Néelz and Pender, 2013), were implemented.These cases correspond to a benchmarking exercise assessing the latest generation of 2D 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/stefan_zipfile.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 (Néelz and Pender, 2013).The modelled area is in the City of Glasgow, UK (Cockenzie Street and surrounding streets) and is approximately 400 by 960 metres.Ground elevations span a range of 21m to 37m.While 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 20 completeness, Fig. 5 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.5m, which is resampled to 2m 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.The model was run until time t = 83 min as this was considered enough to allow the flood to pond in the lower parts of the modelled domain.
Numerical results obtained with Itzï have been compared to those obtained with the implementation of the acceleration solver from LISFLOOD-FP (De Almeida et al., 2012;De Almeida and Bates, 2013).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 6 illustrates the utilised digital elevation model, along with the position of the inflow point (blue triangle) and selected control points for the comparison of model results(red crosses).
The simulation is run for 83 minutes with both LISFLOOD-FP and Itzï using the same parameters, shown in Table 1. Figure 7 illustrates the time series of water level produced by both numerical models and at the eight selected locations.Moreover, for clarity differences between each model result are also shown in all panels, where the blue solid line represents the numerical results of the LISFLOOD-FP, while the green dashed line illustrates results from Itzï, :: as :::: well ::: as ::: the :::: time ::::: series :: of ::::::::: differences.
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 LISLFLOOD-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 comparison of numerical results from both models, :::::: RMSE :: at ::: the :::: eight ::::::: location ::::: range ::::: from ::: 0.1 :: to ::: 5.9 :::: mm :::: (Cf.::::

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 the 25th of June 2007 in the city of Kingston upon Hull, UK.Widespread flooding was registered due to an extreme precipitation event (hp=110mm) that lasted for more than 12 hours.The estimated 24 hour rainfall was estimated between a 1/150 and 1/200 year return period (Coulthard and Frostick, 2010;Hanna et al., 2008).As reported by Yu and Coulthard (2015), the major forcing of this event, was surface water runoff both locally in the urban area and through the rural lands surrounding the city.and the numerical domain is comprised of 3.5millon ::::::: 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.
The numerical forcing condition was defined with measured rainfall recorded by a pluviometer located at Hull University :: the ::::::::: University ::: of :::: Hull, with the assumption that rainfall was uniform over the whole numerical domain.Figure 10 introduces the recorded hyetograph, which has a temporal resolution of one hour.The whole event lasts for 24 hours.Due to the lack 10 of information related to the drainage system within the city, the capacity of the urban drainage to store water was simulated following the work presented by Yu and Coulthard (2015), where this system is represented through a numerical loss which is considered constant in time and uniform in space.The specified value for this loss was uniformly set to 55mm/day.
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 the 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 flood plains, this drag may be conceptually divided into several Grassland : ::: 0.03 zones, namely, main channel, soil grain roughness and vegetative roughness (e.g., Pedrozo-Acuña et al., 2012).Therefore it is important to represent the spatial variability of this parameter.In this sense, the ability of Itzï to employ entry data with heterogeneous resolutions is handy ::::::::::: advantageous.To illustrate this, two numerical experiments are presented, firstly one that considers a uniform Manning's n value of 0.03 for the whole numerical domain, following the work presented by Yu and Coulthard (2015).Secondly, a spatially variable roughness map is defined in terms of the Manning's n value within the grasslands.
It should be noted that although the simulation is run at 5m, the friction map is kept at 30m.The automatic 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 the need of any pre-processing of the data by the user, which is seen as an advantage  Multiple observation data sets for this flood event were obtained from the UK Environment Agency and Hull City Council.
Water depths in some local places were reported to be up to 3 m ::: 3m, but for most affected areas the depth was ::::::: observed : as being less than 1 m :: 1m ::::::::::::::::::::::::: (Coulthard and Frostick, 2010).Flood extent maps were produced by each of these authorities and are shown in Fig. 12, with clear differences between the identified affected areas being noticeable.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 scores (Stanski et al., 1989) has been used for this task.:: A ::::::: succinct ::::::::: description :: of ::::: those :::::: scores : is :::::::: provided :: in ::::: Table :: 3. ::::::: Among Tables 4 and 5 shows the different scores used to compare the computed area and the measured one for each friction case.

5
As mentioned by Yu and Coulthard (2015), this event is not sensible ::::::: sensitive to friction, explaining the similar results in both the uniform and variable Manning's n.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 that ::::: where :::::: correctly :::::::: predicted. :: 1.0 : False Alarm Ratio (FAR) are decreasing; meanwhile the Accuracy and Success ratio are increasing.This could be explained by the predominance of correct negative area that increases together with the threshold.The scores that take into account this 10 bias like Equitable Threat Score, Heidke Skill Score and Critical Success Index are higher when the threshold depth is 20cm.
The reproduction of this real flood event, with the considered assumptions, demonstrate the 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 : .::::: These : are: -Uncertainty in the measured ::::::: observed : flood extent maps, as shown by differences in the areas identified by the UK Environment Agency and the Hull City Council (Cf.Fig. 12), -A simplistic representation of the drainage system, -The utilisation of uniform rainfall in the numerical domain, -The lack of consideration for : of : the infiltration processes, -An inadequate consideration of some terrain features, like walls and buildings, that might not be well represented in a 5m DEM.

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 proven ::::::::::: demonstrated by two analytic benchmarks, resulting in RMSEs one to two orders of magnitude lower than airborne LiDAR vertical accuracy (Hodgson and Bresnahan, 2004).The implementation has also been compared to a reference implementation of the same numerical scheme (LISFLOOD-FP) for test case number 8a of the UK Environment Agency two-dimensional hydraulic model benchmark tests (Néelz and Pender, 2013).On that test simulating a urban flood from direct rainfall and sewer overflow, Itzï gives results nearly equal to LISFLOOD-FP.The presented model has been then ::: 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 works will focus on 1) applying the advantages of the GIS integrated model by using global dataset :::::: 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 urban environment.

Code Availability
The source code of Itzï is freely available under the GNU GPL license (Courty, 2016).The link to the online repository and a user manual is :: are : accessible on the project website (https://www.itzi.org).
support given by a full PhD scholarship from UNAM´s Coordination of Posgraduate 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 10 Univeristy of Leeds for his insights on testing the early versions of Itzï.
Fig. ::: 7).:::: 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.

Figure 8 Figure 8 .Figure 9 .
Figure8introduces the location of the study area within the United Kingdom.Fig.9presents the utilised DEM that was obtained for the area by means of a LiDAR system and has a spatial resolution of 5 metres.The study area measures 87.7km 2

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

Table 1 .
Simulation parameters for the EA test case 8a

Table 4 .
Verification of the flood prediction in Hull using uniform friction

Table 5 .
Verification of the flood prediction in Hull using variable friction