Urban Flood Modeling Using 2D Shallow-Water Equations in Ouagadougou, Burkina Faso

: Appropriate methods and tools accessibility for bi-dimensional ﬂow simulation leads to their weak use for ﬂoods assessment and forecasting in West African countries, particularly in urban areas where huge losses of life and property are recorded. To mitigate ﬂood risks or to elaborate ﬂood adaptation strategies, there is a need for scientiﬁc information on ﬂood events. This paper focuses on a numerical tool developed for urban inundation extent simulation due to extreme tropical rainfall in Ouagadougou city. Two-dimensional (2D) shallow-water equations are solved using a ﬁnite volume method with a Harten, Lax, Van Leer (HLL) numerical ﬂuxes approach. The Digital Elevation Model provided by NASA’s Shuttle Radar Topography Mission (SRTM) was used as the main input of the model. The results have shown the capability of the numerical tool developed to simulate ﬂow depths in natural watercourses. The sensitivity of the model to rainfall intensity and soil roughness coefﬁcient was highlighted through ﬂood spatial extent and water depth at the outlet of the watershed. The performance of the model was assessed through the simulation of two ﬂood events, with satisfactory values of the Nash–Sutcliffe criterion of 0.61 and 0.69. The study is expected to be useful for ﬂood managers and decision makers in assessing ﬂood hazard and vulnerability.


Introduction
In human history, natural disasters have always existed as a common and recurrent phenomenon. Among them, floods are the most widespread worldwide disaster and appear at the top ranking due to important fatalities and damages caused [1][2][3]. According to the United Nations International Strategy for Disaster Reduction [4], flood disasters account for about one-third of all-natural disasters, ranking by economic losses, between 1994 and 2013. The occurrence of flood events in many parts of the world has raised public, political, and scientific interest in flood risks and prevention strategies [5,6]. Despite considerable technologies developed for forecasting, management, and defense, floods are still a significant threat around the world [7]. In a recent study, Petrucci et al. [7], showed that 812 floods recorded killed 2466 people, with important losses, based on 39 years (between 1980 and 2018) of recorded data over nine study areas located in Europe.
In recent decades, Sahelian countries have been strongly affected by high-impact weather events. For example, heavy rains in 2018 affected more than 200,000 victims and caused 52 deaths, where Niamey and its vicinity was the hardest hit area [8]. In addition, this extreme event killed thousands of livestock and destroyed more than 170,000 houses and almost 8000 hectares of crops. In 2009, flooding affected 600,000 people in 12 West African countries, resulting in 193 deaths. The hardest-hit countries were Burkina Faso, Senegal, Ghana, and Niger [9]. Unfortunately, future climate projections The performance of the numerical tool developed was appreciated through the simulation of some real flood events observed in a peri-urban watershed in the city of Ouagadougou. This paper is structured as follows: Section 2 gives a description of the study area, and Section 3 describes the mathematical formalism of flow dynamics in the shallow water, the numerical resolution scheme used, and the input data for application. Section 4 presents the results of the sensitivity analysis and outputs from simulations of real cases of flood events. The conclusions of this study are summarized in Section 5.

Study Area
Ouagadougou, the capital city in Burkina Faso (Figure 1), records important damages linked to floods almost each year. Ouagadougou town has three (3) reservoirs (dams) built since 1963 over Masili river, which drain into the Nakanbe river. A small peri-urban watershed in Ouagadougou was selected for this study because of important damages linked to floods each year and mainly because for two years (2016-2017), flow data monitored at a sub-hourly time scale (15 min) is available. This watershed at the outlet (Latitude 12 • 18 , Longitude −1 • 26 ) covers a surface area of 75 km 2 , with almost 47% of the total area covered by urban settlements (Figure 1b). It is drained by intermittent tributaries of the Massili River. The topographic elevation ranges from 287 to 334 m, with a global slope of 1.45 m/km. The study area is located in a Soudano-Sahelian climate, and the mean annual precipitation is about 788 mm . For that same period, the average daily minimum temperature was 16 • C (Harmattan in December) and the daily maximum temperature was 40 • C in April. In this watershed, there are presidential, industrial, and commercial zones. The population in this zone is estimated at 387,000 inhabitants (estimation in 2020).

Two-Dimensional (2D) Flood Inundation Modeling
Overland surface-flow routing in an urban area is usually simulated by a two-dimensional model [5,32,33] in a wide range of various flow conditions [30] including floods events due to dam failures, heavy rainfall, tidal hydrodynamics, and ocean tsunami hazards. For flood management purposes, this study focused on the estimation of flood patterns (flood extent and flood depth) in a peri-urban area. To achieve this objective, the 2D model seems to be appropriate for this application. According to Alcrudo et al. [36], for such an application, four (4) main steps should be fulfilled for flood propagation model development. These are: (i) understanding of the flow processes relative to the problem; (ii) formulation of appropriate mathematical laws; (iii) development of numerical techniques to solve them and; (iv) validation of model output against experimental and real-life data. The developed model processes are described as follow based on assumptions made by Alcrudo et al. [36].

Mathematical Model Description
The Saint-Venant equations are derived from the laws of conservation of mass and momentum [31,36]. These equations in vector formulation are stated as ∂U ∂t • The functions F and G represent the advective and pressure fluxes, respectively, in the x and y directions: (2) • The terms S z and S f designate, respectively, the bottom slope and the friction terms: g [m/s 2 ] represents the gravitational acceleration and z(x, y) [m] define the soil elevation. Friction terms are computed using Manning's formula [38]: where n [m 1/3 /s] is the Manning-Strickler coefficient, which characterizes the soil roughness.
• S p represents the source term stated as: where P [m/s] and I [m/s] are, respectively, the rainfall intensity and infiltration rate. The infiltration rate was estimated using Horton's model [38]. It is the process of water entrance through the soil profile, and its rate decreases during the rainfall event.
where i 0 [m/s], i f [m/s], and r [s −1 ] are, respectively, the initial infiltration rate, the final infiltration rate, and a constant depending on soil properties.
In this study, the real evapotranspiration is neglected due to the fact that it is an event-based model. It is considered that during the extreme rain event, the real evapotranspiration is very low compared to the volume of precipitation. The model developed uses a structured grid to describe distributed watershed characteristics including topography, urbanization, soil type, and so on. In this model, physical processes of runoff generation are integrated, and particularly, rainfall runoff and infiltration processes are considered over each single grid cell.

Model Numerical Resolution
The model is solved numerically with the finite volume method in a structured grid, hence ensuring proper mass and momentum conservation [45,46]. The geometry of domain Ω is replaced by a structured mesh, where the nodes (grid cell centers) are indexed by integers i and j, respectively, in the x and y axis. The spacing between successive nodes is equal and denoted by ∆x in the x-axis and ∆y in the y-axis.

Grid and Numerical Scheme
The computational grids of domain consist of interior cells that are labeled by i = 1; 2; · · · ; m x and j = 1; 2; · · · ; m y ( Figure 2). The grid is extended by setting a set of two ghost cells at all sides, for i = −1; 0 and i = m x + 1; m x + 2, and forj = −1; 0 and j= m y + 1; m y + 2. At the beginning of each time step, the ghost-cell values are filled based on data in the interior cells and the given boundary conditions. The computational domain is assumed to be large enough, so that there will be only outgoing waves on the boundary of this domain. Hence, along the boundary edge, we set: The resolution in time is split into a sequence of one-dimensional problems.

(a)
In the x-sweeps, we start with a solution U n i,j along each row of cells C i,j , with j fixed. The numerical solution at time t n+ 1 4 is obtained from the resolution of the following scheme: where F n i−1/2,j is an appropriate numerical flux between cells C i−1,j and C i,j .
In the y-sweeps, we use the values of U n+1/4 i,j as data along each column of cells, with i fixed. That results in the calculation of the solution at time t n+2/4 , from the scheme where G n+1/4 i,j−1/2 is an appropriate numerical flux between cells C i,j−1 and C ij . (c) Friction and source terms (rainfall, infiltration, ...) are taken into account using the two following schemes:

Hydrostatic Reconstruction
The HLL approximation [46] presented above is used to build the numerical fluxes at the cell interfaces. Moreover, the water height is reconstructed in order to preserve equilibrium states [44].
(i) The fluxes in the x-axis at the right and left of cell C i,j are defined by where e y = (0, 1, 0) T . The reconstructed solution at the interfaces is defined by with (ii) The fluxes in the y-axis above and under cell C i,j are defined by where e z = (0, 0, 1) T . The reconstructed solution at the interfaces is defined by

Flux at Cell Interfaces
The fluxes F and G are computed by HLL approximation developed by Harten, Lax, and Van Leer [44,[47][48][49]. Following this HLL scheme, the upwind direction of each term of the fluxes is simply dictated by the sign of the flow velocity reconstructed at the interface of the cells. The HLL fluxes at cell interfaces are expressed as follows: where U G , U D , U B , and U U are, respectively, the left, right, down, and up approximations of the solution at the current cell. The coefficients c x i , c y i , i = 1, 2 are defined by The terms λ x i and λ y i , i = 1, 2, 3 are, respectively, the eigenvalues of the Jacobian of F and G. They are defined by The numerical solution at the grid cell indexed by (i, j), at time t n , is the approximation

Application to a Peri-Urban Watershed in Ouagadougou
The model described above is applied to a small peri-urban watershed located in Ouagadougou, Burkina Faso ( Figure 1). Some input data are needed to build and to evaluate a numerical model. In general, four types of information are required: i-topographical data for the description of the geometry of the study area; ii-hydrometric observations for the definition of initial and boundary conditions; iii-effective friction and infiltration parameters for each grid cell; and iv-hydrographs for model validation.

Digital Elevation Model or Topographic Data
The calculation of overland flow and the estimation of the flooded areas require elevation information of the watershed. To simulate the hydrologic processes on the watershed, topographic data or the Digital Elevation Model (DEM) was downloaded from USGS (U.S Geological Survey) website (https://earthexplorer.usgs.gov/), with a spatial resolution of 30 m. From the DEM, regular square grids were built including the area of interest. The coordinates system has been reset based on the lowest values of longitude and latitude coordinates. Thus, 409 × 409 nodes were obtained over the study domain. This topographic data with grid nodes is presented in Figure 3.

Roughness and Infiltration Parameters
Roughness and infiltration rates are influenced by many factors such as soil type, soil moisture, land cover, and land use, among others. Land cover and land use classification data were provided by the Burkina Faso National Institute of Geography (IGB, http://www.igb.bf/). After re-projection, re-sampling, and extraction of land-use classes, these data were overlayed with the watershed grid system. Manning's roughness values were selected from the published literature [50][51][52][53]. Roughness parameters and those of infiltration were collected from the scientific literature because the land-use types considered for our study area are similar to those found in the literature. The values used could be similar to those from field measurements. The roughness values for the overland surface grid cells and interception storage capacity were specified based on land-use classes in Table 1. Likewise, the infiltration map was elaborated by taking into account the maximum infiltration rate of different land cover given in Table 2. Model scenarios were built upon information from Tables 1 and 2.

Rainfall and Water Depth Data
To build the model, synthetic rainfall data were prepared to analyze the model sensitivity to this input variable with the hypothesis of rainfall uniformity over the watershed. To evaluate the model, two real cases of rain events in 2016 were used. Rainfall and flow data used in this study were collected in the framework of AMMA-2050 (https://www.amma2050.org/) project activities.

Model Sensitivity Analysis
Model sensitivity analysis is widely applied in environment modeling problems such as flooding. It consists of testing how model output can be linked to uncertainties from each input parameter [33][34][35]54]. According to Paquier et al. [33] and Kumar et al. [54], sensitivity analysis to model parameters and input data are required to cope with associated uncertainty ranges and also to prove the accuracy and stability of the model concerning the study area. In this study, the model sensitivity analysis was done through three (3) scenarios described as follows: * Scenario 1: Sensitivity to rainfall intensity; * Scenario 2: Sensitivity to Manning's roughness coefficient; * Scenario 3: Sensitivity to imperviousness due to urbanization.
In scenario 1, the topographic data without any settlement (bare soil) was considered, while scenario 2 takes into consideration different land cover types and their corresponding infiltration and Manning's roughness coefficients. The current urbanized surface effect was compared to bare soil conditions in scenario 3, for which the value of roughness coefficient considered were n = 0.015 (urban) and n = 0.020 (bare). The duration for all simulation scenarios was 12 h, while 120 mm was taking for rain intensity during the first hour for both scenario 2 and scenario 3. The sensitivity analysis to rainfall intensity (scenario 1) was carried out by considering four cumulative rainfall values per hour: 50 mm, 80 mm, 100 mm, and 120 mm. For this purpose, it was assumed that the soil type was uniform over the watershed (bare soil for scenario 1), four types of land cover (scenario 2), and bare and urbanized land cover (scenario 3). Thus, Manning's and the infiltration parameters are presented in Table 3 for the first scenario, in Table 4 for the second scenario, and in Table 5 for the third scenario.

Model Validation
After testing the model's robustness through the sensitivity analysis, validation was carried out with two major rains events monitored in 2016 (3rd and 9th of August). The rainfall event on 3 August 2016 occurred within 5 h with a cumulative rainfall of 94.5 mm, and the maximum intensity of a 15 min time-step was 146.5 mm. The 9 August event was two times more important than the 3 August event and generated severe floods in Ouagadougou [28]. In fact, for this event, 210 mm of total rain was recorded within six hours, with a maximum intensity of 122.1 mm at 15 min time steps. Two types of land use were considered: urbanized areas (47% of the total surface area) and natural areas (53% of the total surface area), for which the values of Manning's coefficient were set at 0.015 (urban) and 0.020 (Bare soil). The infiltration parameters are presented in Table 2. The model performance was assessed by using the Nash-Sutcliffe Efficiency [42,55,56] criterion. It is a numerical criteria to measure the difference or gap between observed and simulated data. For this criteria, simulation outputs are accepted as satisfactory if the NSE varies between 0.5 and 1.
where Hsim i and Hobs i represent, respectively, the simulated and observed depth at a given time i, while Hobs is the average observed depth.

Sensitivity Analysis to Rainfall Intensity
The effect of increasing rainfall intensity on water depth in the watershed is shown in Figures 4 and 5. It can be observed in Figure 4 that the extent of flooded areas increases with an increase in rainfall intensity. Figure 5a indicates that the areas flooded under the four synthetic rainfall intensities varies between 4.3 km 2 and 12.7 km 2 when considering depths exceeding 0.25 m. The surface flooded is not proportional to the rainfall intensity. In fact, the flooded area for the 50 mm synthetic rainfall is 4.3 km 2 , but by multiplying the rainfall intensity by 2 (100 mm), the surface flooded is more than twice 4.3 km 2 (10.5 km 2 ). A similar result is observed at the outlet of the watershed, where the highest water depth (1.1 m) obtained corresponds to the largest rainfall intensity (120 mm) ( Figure  5b). The peaks of simulated hydrographs appear between 2 and 3 h, with values varying between 0.32 m for the lowest rainfall (50 mm) and 1.05 m for the highest (120 mm).

Sensitivity Analysis to Manning's Coefficient
The effect of a roughness coefficient change or different types of land cover is shown in Figures 6  and 7. It can be observed (Figures 6 and 7a) that there was a slight difference between the extent of flooding for urban (n = 0.015) and bare (n = 0.02) land-cover types. An important gap is observed between cultivated (n = 0.050) and vegetated (n = 0.065) land cover for the same rainfall intensity (120 mm) when considering flow depths exceeding 0.25 m. For this scenario, the shape of hydrographs (Figure 7b) is almost similar in the first hour for all simulations. Then a slight shift appears between the peaks for urban and bare soil conditions, while an important gap was realized between cultivated and vegetated land cover.

Sensitivity Analysis to Imperviousness Due to Urbanization
Two simulations were carried out to compare the effects of urbanization on flooding in the watershed. Flood extents for both land-use types (urban and bare soil) are presented in Figure 8. Figure 9a presents the cumulative surface flooded over the watershed and flow depth at the outlet for different land uses. It can be observed that water flows between house blocks (Figure 8a) and follows a natural topography slope in Figure 8b. In urbanized surface conditions, the hydrograph peak ( Figure  9b) is observe after 2 h, with a peak depth of 0.57 m, and 3 h for the case of bare soil. The peak depth is more important (1 m) in bare soil conditions compare to the urban surface case.

Observed Events Simulation
Burkina Faso is frequently hit by tropical heavy rainfall events that cause damage in urban cities such in Ouagadougou. In this study, two monitored rainfall events in 2016 (August 3 and 9) were used for the simulation of real events. Flood maps for both events are shown in Figure 10. Observed and simulated hydrographs at the outlet of the watershed are presented in Figure 11. Simulated water depth for both events showed good agreement with observed data at the respective NSE values of 0.61 and 0.69, showing that the developed model reproduces observed data with acceptable quality.

Discussion
This paper was devoted to the development of a 2D shallow water model for urban inundation extent using numerical simulation in Ouagadougou. The developed model's robustness was tested and some real cases simulated. In this study, the developed model parameters (Manning's roughness coefficient, Horton infiltration model parameters) were collected from the literature due to the difficulties in determining them with accuracy. Topographic data was also constructed from coarse resolution (30 m) satellite information. It was therefore then necessary to assess the model for construction robustness or sensitivity to changes in the selected parameters [38,54]. Sensitivity to rainfall intensity (main input) was first tested through flood spatial extent and water level at the outlet of the catchment due to the considerable spatial variation of rainfall intensity in tropical regions [13,14]. The simulations showed that the more rainfall intensity increases, the greater the extent of flooding, and consequently resulting in an increase in water depth at the outlet. This finding corroborates the results obtained by Laouacheria et al. [57], where the impact of the designed rainfall was investigated using hydraulic tools based on the complete Saint-Venant equations within the Storm Water Management Model (SWMM).
The test of model sensitivity to Manning's coefficient was carried out by using different values of Manning's friction coefficient (different land cover) and a unique rainfall intensity for all simulations (120 mm). According to Descroix et al. [25], the rapid land-cover change over the few decades may explain the high frequency [14,26] of floods in most Sahelian regions of West Africa. Increasing values of n were been considered at 0.150, 0.020, 0.050 and 0.065, which correspond, respectively, to urban, bare, cultivated, and vegetated land cover. Simulation results demonstrate a weak impact of increasing n values on flood extent ( Figure 6) but an important shift in the time of the water depth peaks. It is then well established that the sensitivity of the developed model to Manning's friction coefficient variation is weak in flood extent but high in flow depth at the outlet. These results are similar to those obtained by Kumar et al. [54], where the increase (+10% and +20%) or decrease (−10% and −20%) of Manning's roughness coefficient collected from published literature were found to be affecting water depth.
Another important finding this study demonstrates is the fact that water depth at the outlet is almost similar for urban and bare land covers, with high depth as compared to cultivated and vegetated land-cover type (Figure 7b). For the purpose of flood management, the results demonstrate the important impact of vegetating urban cities or to adopting "green cities" for urban planning to help increase infiltration and thus reduce the stagnation of water over a long period of time.
Model sensitivity to actual urban settlement was tested by comparison with the case where the catchment was considered as uniform, bare land cover. Urbanization and settlement was observed to impact flow by reducing the flood extent and increasing the flood peak at the outlet.
Available data concerning two rainfall events and their corresponding water levels observed at the outlet of the basin in 2016 were used for model calibration and validation. For those rainfall events, the lack of information about flood extent constrained us to compare only the simulated and observed water level as a validation process for the developed model. A good relationship between the observed and simulated water level was shown, with NSEs of, respectively, 0.61 and 0.69 for events 1 and 2.

Conclusions
This study focused on the development of a 2D numerical model for urban flood modeling in West Africa. The case study considered is a peri-urban area in Ouagadougou (Burkina Faso). The model is based on shallow-water equations. The numerical solutions are obtained from the finite volume by using structured grids. A sensitivity analysis of this model to input data (rainfall intensity) and its parameters (roughness coefficient and infiltration) revealed the capability of the model developed to simulate flow. Increasing rainfall intensities in the bare soil condition result in increasing water depth within the watershed and in the outlet, and an increase of n affected the hydrograph peak more than the spatial extent of flooding. Imperviousness due to urbanization reduced the flood extent and increased water depth at the outlet. Even though the model developed based on the coarse-resolution DEM (30 m) showed some gaps between observed data and simulated values, the capability of the model to simulate flow over an urban drainage network was demonstrated. The precision of model simulation may depend on DEM accuracy and field soil characteristics but also relies on hydrologic and hydraulic model structure. Even though this model in construction necessitates more improvement (for instance, taking into account drainage network characteristics and testing in various other conditions), this tool may be of great interest for decision makers (urban planners, flood-impact managers, etc.) in Burkina Faso in general, and particularly in Ouagadougou, as it allows the production of maps of maximum water level and flood extent, which is indispensable for effective flood management. This study lays out two important knowledge gaps recommended for further studies: • Flood spatial extent has not been validated because of the lack of field data for that process, and thus, high-resolution aerial photographic data of flood extent for flood map validation deserve to be investigated; • Sensitivity analysis to digital elevation model resolution for topographic input generation needs to be carried out with a high resolution of 1 m or less to improve accuracy.