Flood risk assessment due to cyclone-induced dike breaching in coastal areas of Bangladesh

Bangladesh, one of the most disaster-prone countries in the world, has a dynamic delta with 123 polders protected by earthen dikes. Cyclone-induced storm surges cause severe damage to these polders by overtopping and breaching the dikes. A total of 19 major tropical storms have hit the coast in the last 50 years, and the storm frequency is likely to increase due to climate change. The present paper presents an investigation of the inundation pattern in a protected area behind dikes due to floods caused by storm surges and identifies possible critical locations of dike breaches. Polder 48 in the coastal region, also known as Kuakata, was selected as the study area. A HEC-RAS 1-D–2-D hydrodynamic model was developed to simulate inundation of the polder under different scenarios. Scenarios were developed by considering tidal variations, the angle of the cyclone at landfall, possible dike breach locations and sea level rise due to climate change according to the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC). A storm surge for a cyclone event with a 1-in-25-year return period was considered for all the scenarios. The primary objective of this research was to present a methodology for identifying the critical location of dike breaching, generating a flood risk map (FRM) and a probabilistic flood map (PFM) for the breaching of dikes during a cyclone. The critical location of the dike breach among the chosen possible locations was identified by comparing the inundation extent and damage due to flooding corresponding to the developed scenarios. A FRM corresponding to the breaching in the critical location was developed, which indicated that settlements adjacent to the canals in the polders were exposed to higher risk. A PFM was developed using the simulation results corresponding to the developed scenarios, which was used to recommend the need of appropriate land use zoning to minimize the vulnerability to flooding. The developed hydrodynamic model can be used to forecast inundation, to identify critical locations of the dike requiring maintenance and to study the effect of climate change on flood inundation in the study area.

The frequency and intensity of the cyclones around the world are likely to increase due to climate change, which will require resource-intensive improvement of existing or new protection structures for the deltas. The identification and prioritization of the maintenance of critical locations of dike breaching can potentially prevent a disaster. The use of non-structural tools such as land use zoning with the help of flood risk maps and probabilistic flood maps has the potential to reduce risk and damage. The method presented in this research can potentially be utilized for deltas around the world to reduce vulnerability and flood risk due to dike breaching caused by cyclone-induced storm surge.

1 Introduction

Bangladesh is situated in a low-lying delta of three major rivers: Ganges, Brahmaputra and Meghna. A total of 80 % of the country's land is located below 10 m a.m.s.l. (above mean sea level) (Heitzman and Worden, 1989), and it is formed of sediments carried by the above-mentioned rivers. The population of Bangladesh was about 131.5 million by the year 2000 (World Bank, 2018), of which about 49 % were living in coastal zones (Neumann et al., 2015). The coastal areas of Bangladesh are flooded frequently due to cyclone-induced storm surges and occasionally due to high water levels in the rivers caused by heavy rainfall in the upstream catchments of Ganges, Brahmaputra and Meghna. The coast was hit by five severe cyclones between 1995 and 2010, causing flooding, huge damage and loss of life (Dasgupta et al., 2014).

Bangladesh has 123 polders in the coastal area, each surrounded by earthen dikes, which are designed to protect the inland from flooding due to high tides. The existing crest level of these dikes is only adequate enough to protect the coastal area from cyclones with 5- to 12-year return periods (Islam et al., 2013). These dikes usually get damaged and sometimes breached by tropical cyclones of high intensity, which causes flooding inside the protected areas, damages to properties and loss of life. For example, Cyclone Sidr hit the coast of Bangladesh in 2007, affecting 8.9 million people and causing USD 1.7 billion of damages (GOB, 2008; Dasgupta et al., 2014). In 2009, Cyclone Aila affected 3.9 million people, with estimated damages of USD 270 million (EMDAT, 2009).

Crest levels of the coastal dikes were recently designed for an event with a 25-year return period under the Coastal Embankment Improvement Project (CEIP) (BWDB, 2013). A storm surge event with a 25-year return period was considered in this study for the generation of different scenarios. Under the CEIP, the crest levels of the dikes were designed considering wave actions, astronomical tides and the required freeboard. Raising the crest level was considered the only mitigating measure. Various studies on the coastal areas of Bangladesh (e.g. Karim and Mimura, 2008; IWM, 2005; Azam et al., 2004; Madsen and Jakobsen, 2004; CSPS, 1998; Flather, 1994) considered flooding only due to overtopping of the dikes during storm surges. The effects of breaching of the dikes due to piping and scouring on the landside during cyclones have not been studied. The coast of Bangladesh is frequently hit by severe cyclones (five cyclones between 1995 and 2010, Dasgupta et al., 2014). The Bangladesh Water Development Board (BWDB) is responsible for the operation and maintenance of these dikes and lacks fund to conduct proper repair of damaged dikes subsequent to any severe cyclone. As a result the dykes remain vulnerable to breaching. Identifying the critical location(s) of dike breaching and prioritizing the repair of the critical location are likely to reduce the breaching possibility.

Moreover, non-structural measures for flood risk management such as land use zoning using a flood risk map (FRM) and a probabilistic flood map (PFM) to locate the vulnerable areas are currently unavailable for the coastal areas of Bangladesh. Flood zoning can be a useful risk mitigation measure as land use governs the exposure and may aggravate the hazard (Barredo and Engelen, 2010).

Furthermore, the intensity and frequency of these tropical cyclones are likely to increase in the future due to climate change. It is projected that by the year 2100, the frequency of the most intense cyclones will increase substantially and the intensity of tropical cyclones will increase by 2 % to 11 % due to global warming (Knutson et al., 2010). Flooding by tropical cyclones will also increase in the future as a result of sea level rise (SLR) (Woodruff et al., 2013). SLR and sea surface temperature (SST) will affect the cyclone-induced storm surge height in the Bay of Bengal (Karim and Mimura, 2008). With increasing SST, the storm surge height may increase from 21 % to 49 %, and with SLR, the flood depth due to storm surges may increase by 30 %–40 % (Karim and Mimura, 2008). The land subsidence in the delta will exacerbate the effect of SLR. By the year 2100 the annual estimated damage due to tropical cyclones may increase by USD 53 billion (Mendelsohn et al., 2012).

At present, a flood forecasting system is not available for the coastal region of Bangladesh. The BWDB, which is mandated to protect the area, does not have a clear picture about the inundation patterns corresponding to various climatic conditions. Moreover, identifying zones in the embankment critical to flooding in the polder will help BWDB in prioritizing their maintenance. This paper presents a methodology to identify the critical location of dike breach due to cyclones, generating a FRM and a PFM for the breaching of dikes by cyclone-induced storm surges. Different scenarios of storm surges were formulated by considering storms of different frequencies with varying tidal conditions, the angle of the cyclone at landfall and SLR. A cyclone event with a 25-year return period was considered in this research. A coastal polder (Polder 48) of southern Bangladesh was selected as the study area.

Resource-intensive adjustment of the protective structures for the deltas around the world will be required as the frequency and intensity of cyclones will increase with climate change. Along with structural measures, non-structural tools, such as land use zoning with the help of flood risk maps and probabilistic flood maps, have the potential to reduce risk and damage. The identification of the critical locations of breaching and intensification of the maintained effort for these locations can potentially prevent a disaster. The method presented in this research can be utilized for vulnerable deltas around the world, even though the coastal region of Bangladesh was selected as the case study area.

2 Study area

Polder 48, which was considered as the study area for this research, is surrounded by dikes and has a sea-facing dike of about 20 km length on the southern side of the polder. The polder is located on the south-western coast of the Bangladesh Delta (Fig. 1), stretching from 21^∘50^′28^′′ N 90^∘05^′17^′′ E to 21^∘50^′06^′′ N 90^∘14^′14^′′ E. The outline of the study area in Fig. 1 (also in Fig. 3) depicts the dike alignment of the study area as well. The area is also known as Kuakata, and it is in the administrative zone of the Kalapara Sub-district (Upazilla) of Patuakhali District. It has an area of 50.75 km² with 24 240 inhabitants according to the 2011 census (BBS, 2012). Most of the inhabitants are farmers and fishermen (Nasreen et al., 2013). Shrimp culture and tourism are also part of the economic activities. The land use is classified by the Ministry of Land of Bangladesh into the following four classes: rice fields, settlements, shrimp ponds and water bodies (rivers and canals). Climate and agricultural practices of Kuakata are similar to the climate and agricultural practices of the country (Bangladesh). The average yearly rainfall in Kuakata is 2590 mm (Climate-Data, 2016), and the annual average temperature is 25.9 ^∘C (Climate-Data, 2016). Rabi (November–February), Kharif-I (March–May) and Kharif-II (June–October) are the three seasons for growing crops (DAE, 2009). The elevation of 80 % of the area is 1.55 m below PWD, the vertical datum established by the Public Works Department of Bangladesh, which is 0.46 m below mean sea level. The land level surveys at different times have indicated that this polder is facing land subsidence issues. Brown and Nicholls (2015) reported the estimated mean subsidence rate of the Ganges–Brahmaputra–Meghna (GBM) delta to be 5.6 mm yr⁻¹, with an overall median of 2.9 mm yr⁻¹.

Figure 1 Location map of the study area, Polder 48 (Kuakata).

The area was severely affected by the recent storms Sidr, Aila and Mohasen in 2007, 2009 and 2013 respectively. For example, during Cyclone Sidr, 94 people died and 45 % of the crops were lost in the Kalapara Sub-district (Ahamed, 2012).

Figure 2 Methodological approach followed in this study.

Andharmanik, Galachipa and Khaprabhanga rivers are in the east, west and north of the study area respectively, whereas the Bay of Bengal is on the southern side of the study area (Fig. 1). Galachipa River is the widest among the rivers surrounding the area. On the southern side, the study area has a seashore of 20 km width, which is partly protected by the mangrove forest at several locations. There is a narrow sea beach on the south-western side of the area. The western part of the sea-facing dike was overtopped during Cyclone Sidr, causing flooding inside the polder (Hasegawa, 2008). The loss of livestock and food grains was such that it created partial deficiency of food in Kuakata (TANGO International, 2010). The average crest level of Polder 48 on the northern side is 4.5 m PWD, and on the southern side (sea-facing side), it is 6 m PWD (Islam et al., 2013). The existing embankments of 17 polders of the region, including Polder 48, were redesigned and rehabilitated during the first phase of the CEIP (Islam et al., 2013). The CEIP proposed a crest level of 7.36 m PWD for the dike of Polder 48 (Islam et al., 2013).

3 Methodology

The methodology followed is presented in Fig. 2 and described in the following sections.

3.1 Setting up of 1-D–2-D coupled model

In order to build a 1-D–2-D inundation model, field measurements (land level surveys, observed water levels, canal alignments and cross sections of the river and canals) and information from remote sensing (satellite imagery) were gathered (Fig. 3). The Institute of Water Modelling (IWM) of Bangladesh collected hydraulic, hydrologic and land-level data of the study area (along with other polders) in the framework of the feasibility study of the Coastal Embankment Improvement Project (CEIP). The IWM has kindly provided the measurement data for the study area.

The digital elevation model (DEM) was generated by combining the land-level surveys conducted by the IWM and FINNMAP. The land-level survey by the IWM (conducted in 2012) did not cover the whole study area. FINNMAP conducted a topographic survey of the study area in 1988 (MIWF, 1993). The differences in elevation between land surveys by IWM and FINNMAP indicated the land subsidence. An average subsidence was computed, which was used to update the elevations of the FINNMAP survey for the areas within Polder 48 for which survey data from the IWM were not available. The combined DEM has a resolution of 50 m. The same DEM was used for the simulations of the year 2100 without any corrections for further subsidence. Subsidence of the coast in the past has been reported by Brown and Nicholls (2015) and was verified with the survey data from the IWM and FINNMAP. Subsidence may continue in the future, but in the absence of scientific studies it was not considered for the future scenarios in this research. It is noteworthy that if subsidence continues, then the effect of the SLR may be increased, and the results reported in this research should be treated to some extent as underestimated values.

The bathymetry of the sea near the coast was collected from the Global Bathymetric Chart Of the Oceans (GEBCO) (Smith and Sandwell, 1997). The land use data were collected from the Ministry of Land of Bangladesh. MODIS reflectance data were used for the analysis of previous flood events. The methodology and equations suggested by Hoque et al. (2015) were used to analyse the MODIS reflectance data to determine flood extents during previous flood events. The intention was to utilize the flood extent generated from MODIS reflectance for calibration of the hydraulic model. However, no flood images from MODIS were available during the simulation period.

The river analysis tool HEC-RAS (version 5.0) from the US Army Corps of Engineers was used to develop the 1-D–2-D coupled inundation model. The flow in the river was modelled in 1-D, whereas the flow over the floodplain was modelled in 2-D. HEC-RAS 5.0 is a free tool which can simulate 1-D, 2-D and 1-D–2-D coupled models for steady and unsteady flow. The 2-D module of HEC-RAS provides the option to simulate flow of water either with the diffusion wave equation or with the full shallow water equation (St. Venant equation). The availability of irregular flexible mesh in HEC-RAS and the option for faster simulations led to the selection of HEC-RAS 5.0 as the modelling tool. Data utilized for developing the model and their sources are presented in Table 1.

Table 1 The data used in developing the mathematical model and their sources. IWM, GEBCO and JSCE stand for the Institute for Water Modelling, the General Bathymetric Chart of the Oceans and the Japan Society of Civil Engineers respectively.

The 1-D section of the model was developed and calibrated using the information shared by the IWM. The 1-D part of the developed model was calibrated for non-flood conditions as measured discharge and water level data during a cyclone event were unavailable. The model was simulated using discharge as the west boundary and water level as the east boundary conditions (Fig. 3). The calibrated 1-D model was then coupled with the 2-D model of flow over the floodplain using the DEM of the study area.

For the 1-D–2-D inundation model, a computational mesh with a flexible shape, was developed in HEC-RAS (Fig. 3). HEC-RAS generates meshes with irregular shapes. The rectangular cells of the developed 2-D mesh had a resolution of 25 m, and the non-rectangular cells had areas ranging from 625 to 1282 m². The roughness coefficient (Manning's n varying from 0.025 to 0.05) was provided according to the land use of each cell. A sensitivity analysis as suggested by Hall et al. (2005) was carried out by varying Manning's roughness coefficient n before the calibration of the 2-D inundation model.

Figure 3 Schematic diagram of the study area with location of control structures and gauges and the considered breach locations.

Building and calibrating the 1-D model was the preliminary step for developing the 1-D–2-D coupled model. The water bodies surrounding the study area were included in the 1-D model. The study area has Khaprabhanga River on the northern side and the sea on the southern side (Fig. 3). The connection of Khaprabhanga River with other rivers was not considered in the model. This was due to the fact that storm surges are observed during the pre- or post-monsoon periods, whereas fluvial floods are observed during the monsoon. Flow through rivers did not play a major role during the previous cyclones. The western and eastern side of the embankment have mangrove forests between the rivers and the embankment (Fig. 3).

For the river, the surveyed cross sections were used in the 1-D model (Fig. 4). The storm surge on the sea was conceptualized as a water surface profile in a 1-D channel on the southern side of the study area (Fig. 3). The GEBCO bathymetry (Fig. 5) was used for the channel. An alternative was to develop a 2-D model for the coastal hydrodynamics. However, as the coast of Bangladesh is flat and shallow a large area of the sea would have been included in the model. As the focus was on studying the inundation of Polder 48 and not the coast, we followed a simpler representation of the storm surges using a 1-D model. The synthetic water level data for boundaries of the model were generated by following the tidal water level pattern and the storm surge height considered for all the scenarios (Table 2). The water surface profile corresponding to each scenario (Table 2, discussed in Sect. 3.2) was considered as the profile in the 1-D model of the seaside (Fig. 6).

The dense canal network of 122 km, inside the study area, is connected with the Khaprabhanga River, which regulates the in- and outflow into the river network through a system of 13 control structures. The regulators remain closed during cyclones, making the canal network isolated. Therefore, the canal network inside the polder was not included in the 1-D model. However, the simulation of the overland flow consequent on breaching of the dike will be affected by the canal geometry, and therefore, the wider and larger canals were included in the DEM.

Figure 4 A typical cross section of the Khaprabhanga River.

The geometry and propagation of the breach of the dike depend primarily on the storm surge height, the angle of landfall, soil properties and wave action. The coastal embankments of Bangladesh are usually earthen. The geometrical properties of the breaching of the dike and the time required for breaching were calculated following the instructions of the US Bureau of Reclamation. An S curve was used for breach propagation with time (Oumeraci, 2006). As the geometry of the breach is not independent, it was not considered as a parameter for scenario development.

Figure 5 Cross section for the 1-D network on the seaside.

In order to ensure model stability, a maximum spacing between the computational points was imposed and computed using Samuels' formula (1989), presented in Eq. (1):

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }x\le \mathrm{0.15}\times D/S_{\mathrm{0}},\end{array}$$

where Δx is the spacing between the computational points, D is the average bank full depth of the channel and S₀ is the average slope of the channel. The maximum spacing between cross sections was calculated to be 300 m. The river had a steeper bed slope than the long shore slope of the sea bathymetry, requiring smaller Δx to ensure stability, and the same Δx will reduce instability of the foreshore as well.

Figure 6 Variation of water level at three locations (chainage 0, 12 940 and 20 300) along the 1-D channel on the seaside according to a specific scenario (out of 72 scenarios).

As suggested by Fromm (1961), the Courant number was kept less than or equal to 1.0 to maintain the stability of the numerical model by controlling the time step. The Courant number was calculated using the following Eq. (2):

$$\begin{array}{}\text{(2)}& Cr=V\times \mathrm{\Delta }t/\mathrm{\Delta }x,\end{array}$$

where Cr is the Courant number, V is velocity, Δt is the time step and Δx is the spacing between the cross sections.

3.2 Cyclonic scenarios considered

Different scenarios were developed considering the probability of the occurrence of cyclones, the angle of landfall, SLR due to climate change, diurnal, semi-diurnal and seasonal variation of tides, locations of breaching of the dike and geometrical properties of the breach.

• Frequency of the cyclone. A cyclone with a 1-in-25-year return period was considered for all the scenarios as this is used as the design criteria for the dikes (BWDB, 2013). A total of 19 previous cyclones for different tidal conditions were simulated by the IWM using a 2-D model for the Bay of Bengal. A statistical analysis was conducted using these model results to generate the storm surge height corresponding to a cyclone with a 25-year return period (Islam et al., 2013). Due to a lack of data, change in the probability of the occurrence of cyclones in the future was not considered.

• Angle of landfall. The angle of landfall affects the height of storm surges. The storm surge height increases with angle of the storm to the coastline (Azam et al., 2004). The angle of attack governs the wind speed, which is one of the parameters for the height of cyclone-induced storm surges (Azam et al., 2013).

• Tides. The difference between the storm surge at high tide and low tide is 1.2 m for the study area (Azam et al., 2004). The average seasonal variation of the tidal range is 1.3 m.

• Sea level rise. The coast of Bangladesh may be severely affected by SLR, and one-quarter of the land may be lost due to SLR by 2100, which will directly affect 3 million people (Ericson et al., 2005). IPCC published their Fifth Assessment Report (AR5) in 2013. Among the scenarios considered in AR5, RCP2.6 (Representative Concentration Pathway 2.6) is the most optimistic one and RCP8.5 is the worst considering the carbon emission, rise in temperature and SLR. The mean SLR at the end of 21st century is estimated to be 0.4, 0.47, 0.48 and 0.63 m for RCP2.6, RCP 4.5, RCP6.0 and RCP 8.5 respectively (Stocker et al., 2013). For this study, RCP8.5 with SLR of 0.63 m was considered for developing the scenarios.

• Location of breach. The sections of the sea-facing dike of the study area protected by mangrove forest, sand dunes and a wide beach are least likely to be breached due to storm surges. The study considered breach locations with the least protection. The locations considered for dike breaching as well as the mangrove forest around the study area are shown in Fig. 3.

A scenario matrix consisting of 72 scenarios was generated by combining different phases of tides, angle of landfall, SLR and breach locations (Table 2). A single breach was considered for all the scenarios. The highest storm surge height among all the developed scenarios was 7.2 m PWD, considering the angle of landfall to be 230^∘, a high tidal phase during spring tides and SLR and dike breaching at any of the chosen locations. The highest storm surge height as the boundary condition with breaching in the western, central and eastern parts of the dike was considered to be the worst-case scenario and was denoted as Scenario S1, S2 and S3 respectively (Fig. 2). Flooding due to overtopping of the dikes was not considered as the crest level (7.36 m PWD) was higher than the highest storm surge height (7.2 m PWD).

Table 2 Storm surge heights corresponding to different scenarios considered. The bold values are the storm surge height for the worst case scenarios.

To identify the critical locations of breaching, results of the scenarios simulated with HEC-RAS corresponding to the three worst-case scenarios S1, S2 and S3 were compared based on the total area flooded and estimated damage due to flooding. Using the calculated damage and probability of occurrence of the event, a risk map was generated for the critical locations of the sea-facing dike. A probabilistic flood map (PFM) was generated from the flood maps of the 72 scenarios (Table 2; Output 2 in Fig. 2). As the storm surge height suggested by Islam et al. (2013) corresponds to an event with a 25-year return period, the PFM generated in this study corresponds to a 1-in-25-year return period.

3.3 Estimation of damage due to floods

A comprehensive damage calculation should involve both direct and indirect damage due to floods (Büchele et al., 2006). Direct damage is caused by physical contact of properties and human beings with floodwater. Indirect damage is caused by interruption of services, production and transportation and degradation of health due to floods. Due to a lack of data, only the direct damages to properties were calculated for the study area. The damage was considered a function of flood depth. The land use of the study area was classified by the Ministry of Land of Bangladesh as settlements, rice fields, shrimp ponds and water bodies (rivers/canals). Only the tangible damage was considered, and no environmental damage was calculated. Damage to the canal network was not considered. The damage in a flood event was calculated using Eq. (3).

$$\begin{array}{}\text{(3)}& D=\left({\sum }_{i=\mathrm{0}}^n{x}_i\times f\left(x_i\right)\right)\times A_i,\end{array}$$

where D is the total direct tangible damage in a flood event, n is the total number of computational cells within the flooded area, x_i is the flood depth of cell i, f (x_i) is the damage function for the land use of the flooded cell i and A_i is the area of cell i.

Depth–damage curves for different land classes for the study area were developed by adapting depth–damage curves found in the literature (Fig. 7). Reese et al. (2010) calculated flood damage as a percentage of the property value of buildings categorized based on the construction material. The buildings of the study area are primarily built of timber due to its low cost and easy availability. The depth–damage curve suggested by Reese et al. (2010) for buildings made of timber was used as a basis for generating the depth–damage curve for the settlements (residential area). Simple Action for the Environment (SAFE) carried out research on the average value of properties in rural areas of Bangladesh (SAFE, 2011). These property values were used to update the damage values used by Reese et al. (2010). Muktadir and Hasan (1985) reported that rural houses of Bangladesh are built with a large courtyard, and, as a result, houses have a lot of open and unoccupied space around buildings. The damage curve considered for the residential area was used for the damage to the buildings and not for the courtyard. Moreover, the satellite image of the area also indicated that about half of the settlement was without buildings. A satellite image from Google Earth was used for analysis. The satellite image of the area was downloaded and georeferenced. Then, the areas for buildings and open areas for households were manually calculated using ArcGIS. Therefore, 50 % of the settlement area was considered to have no damage.

The cultivation of rice involves flooding the rice field with water up to a few centimetres. However, if the height of water increases and the rice plant goes under water, then the productivity decreases. The damage to rice plants also depends upon the flood duration. If the rice plant is continuously under water for more than 2–3 days, then the damage can be up to 80 % (Chau et al., 2014). The simplified (with regards to flow velocity and flood duration) depth–damage curve for rice fields suggested by Chau et al. (2014) was used in this study (Fig. 7).

Shrimp ponds are surrounded by embankments so that there is no damage to shrimp ponds till the flood level crosses the embankment level. However, when the flood level is higher than the embankment level, shrimps escape causing a loss of the total investment. To take this into account, the investment made by farmers was assessed using a study conducted by Fatema et al. (2011). According to the study the investment for shrimp pond in the study area was about EUR 0.09 m⁻². Based on the practices in the study area, the banks of the shrimp ponds were considered to be 2 m above the adjacent land, and the depth–damage curve (Fig. 7) was modified accordingly.

Figure 7 Depth–damage curves for different land use classes.

The damage calculations were carried out using ArcGIS. The simulated flood depth and land use for each grid cell were used as input, and the damage in each grid cell was computed using the depth–damage curve corresponding to that land use. The damage for each scenario was estimated using this procedure.

3.4 Calculation of flood risk and generation of risk map

Flood risk assessment is an essential part of risk management. Spatial distribution of risk and areas requiring mitigation measures can be identified from flood risk maps. To examine the spatial variation of risk, flood risk analysis was carried out and a risk map was generated considering dike breaching at the critical locations. Van Manen and Brinkhuis (2005) and Klijn (2009), as part of the FLOODsite project, carried out research to quantify the flood risk for the polders in the Netherlands for dike failure defining the risk as a product of the probability of occurrence of the event and the consequences which was defined by Helm (1996). Equation (4) was used to calculate the risk due to flooding:

$$\begin{array}{}\text{(4)}& R=P_{\mathrm{F}}\times S,\end{array}$$

where R is risk, P_F the probability of occurrence of the flood hazard and S the consequences.

The exceedance probability (return period) of the cyclone-induced storm surge was used as the probability of occurrence of the hazard. The probability of flooding within a protected area is not the same as the probability of the hazard and depends also upon the probability of failure of the dike. It is a difficult probability to compute as the probability of dike failure also depends upon the dike maintenance, about which information was not available. Here we have assumed that the probability of occurrence of the hazard and the probability of failure of dike are the same.

3.5 Probabilistic flood map

Purvis et al. (2008) stated that the risk assessment for the most probable scenario cannot take into account the impact of the scenario of low probability, stressing the necessity of a probabilistic risk analysis. The equation suggested by Purvis et al. (2008) for probabilistic risk analysis was adjusted and used for this research to calculate the probability of flooding of each cell and is presented below in Eq. (5):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{P}_{i(i=\mathrm{1}\mathrm{to}N)}=\\ \text{(5)}& {\displaystyle }& {\displaystyle \frac{{\sum }_{j=\mathrm{1}}^M{F}_{ij}\times P_{fj}}{{\sum }_{j=\mathrm{1}}^M{P}_{fj}}}\text{and}{F}_{ij}=\left\{\begin{array}{ll}\mathrm{1},& \text{if flooded}\\ \mathrm{0},& \text{if dry}\end{array}\right.,\end{array}$$

where P_i is the probability of flooding at cell i; P_fj is the probability of reaching a certain storm surge level in simulation number j; F_ij is the binary value indicating if the cell i is flooded or not in simulation j; and $j=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },M$, where M is the number of scenarios considered (=72) and $i=\mathrm{1},\mathrm{2},\mathrm{3},\mathrm{\dots },N$ are the computational grid cells on the polder area and N is the number of cells.

Equation (5) was used in this study to calculate the probability of flooding in each cell. The probabilistic flood map (PFM) was calculated using the results of all the scenarios.

4 Results and discussion

The developed 1-D–2-D model for the present study was calibrated for the 1-D part by comparing the observed and simulated values for discharge and water level. The corresponding performance indicators used for evaluation were the coefficient of determination (R²), the root mean square error (RMSE) and the mean absolute error (MAE), for which values of 0.98, 2.15 and 1.68 m³ s⁻¹, respectively, were obtained for discharge, and 0.98, 0.09, and 0.08 m, were obtained for water level, respectively. The average values of the discharge and water level for the considered simulation period were 5.68 m³ s⁻¹ and 0.82 m respectively. The period of simulation for calibration coincided with the surges corresponding to Cyclone Sidr (from 14 to 17 November 2007). The simulation results indicate that the dike facing the seaside was overtopped and the area inside the polder was inundated. This conclusion is in line with the survey conducted by the Japan Society of Civil Engineers (JSCE) (Hasegawa, 2008).

The coupled 1-D–2-D model has not been calibrated because there were no flood maps showing flood extents available for recent cyclones. However, the 2-D part of the model was pseudo-calibrated considering MODIS reflectance data. Such data were used in order to analyse the inundation extent, though this also posed considerable challenges due to the cloud coverage during the cyclones. The survey conducted by JSCE after Cyclone Sidr aimed to investigate the flood extent and depth, but only provided flood depth for one location inside the study area. This location was used for the calibration of 2-D model. The difference between the reported and the simulated flood depth was 4.5 %. Prior to the calibration of the 2-D model, sensitivity analysis was carried out regarding the roughness coefficient (Manning's n). The analysis indicated that the inundation model is not highly sensitive to the roughness coefficient, and the areas of low flows (locations furthest from the dike breach) are most sensitive. The sensitivity analysis was done for the breaching in the western part of the dike only. It was considered that the breaching in other locations will have similar effects as the area inside the polder is flat and low-lying with mostly farmlands near the dike.

This 1-D–2-D model, which had limited calibration points, was further used in simulating the developed scenarios. The simulated results were used to analyse flood depth, extent and damage due to flooding. The FRM and the PFM were generated based on flood results of the model.

4.1 Inundation corresponding to three worst-case scenarios

Among the simulated scenarios, the results of three worst-case scenarios (Scenario S1, S2 and S3) were compared to identify the critical location of breaching. The corresponding flood maps for the worst-case scenarios are presented in Fig. 8.

Figure 8 Flood extent corresponding to three worst-case scenarios of dike breaching in the central, eastern and western section of the dike.

Flood extents corresponding to all different scenarios presented in Table 1 were compared to understand the effect of SLR, diurnal and seasonal tidal variation and the angle of cyclone at landfall. The flood extents of different scenarios considering the breaching in the central part of the sea-facing dike are presented in Fig. 9.

Figure 9 Comparison of flooded areas corresponding to different scenarios considering the breaching in the central part of the sea-facing dike.

Moreover different land classes were considered while computing the flood extent for the three worst-case scenarios (Table 3). The analysis of the flood extent for different flood depths, based on the considered land uses, is presented in Fig. 10. The highest storm surge height among all the developed scenarios was 7.2 m PWD (Table 2). This storm surge height with breaching at the western, central and eastern parts of the dike was considered as the worst-case scenario and was denoted as Scenario S1, S2 and S3 respectively.

Table 3 Flooded areas of different land classes corresponding to the three worst-case scenarios.

Figure 10 Flooded areas for different ranges of flood depths corresponding to different scenarios.

4.2 Comparison of calculated damages

The damage due to flooding was calculated using the depth–damage curves for different land classes. The calculated damage for different land classes and damage for different flood depths corresponding to the three worst-case scenarios are presented in Table 4 and Fig. 11.

Table 4 Calculated flood damages for different land classes corresponding to different scenarios.

Figure 11 Variation of estimated flood damages with varying ranges of flood depths corresponding to Scenario S1, S2 and S3.

Figures 10 and 11 correspond to the flooded area and damage due to different ranges of inundation respectively. The flood area and damage were highest for inundation depths of 0.5 to 1.0 m.

4.3 Risk map for the worst-case scenario

The flood risk map for the scenario with the critical locations of dike breaching is presented in Fig. 12. The risk map presents the assessed risk of flooding due to breaching at critical locations of the dike. Comparison of the flooded area and damage due to flooding for the three worst-case scenarios led to the identification of Scenario S1 as the critical location of breaching. The identification of the critical location of breaching is described in Sect. 4.5.

Figure 12 Flood risk map corresponding to the breaching at the critical location of the dike. The following three classes of risk are shown: high, medium and low. The four land uses considered are shown as well.

4.4 Probabilistic flood map

Although the inundation maps are widely used for spatial planning and flood mitigation measures, the uncertainty of mathematical modelling affects the output of inundation maps (Alfonso et al., 2016). In order to account for uncertainty, probabilistic flood maps are suggested to be used (Domeneghetti et al., 2013). The probabilistic flood map was calculated from the inundation maps corresponding to the 72 scenarios considered in the study. Probabilistic flood maps were calculated for a threshold of flood depth greater than 0.5 m. The developed damage curves suggest that the damage for flood depth below 0.5 m is minimal. Moreover, considering the widely accepted “living with floods” philosophy in Bangladesh, a threshold of 0.5 m was adopted. This threshold was used in developing the PFMs. This threshold was not considered while the estimation of damage due to flooding was conducted. The calculated probabilistic flood map is presented in Fig. 13. The probabilistic flood map indicates the likelihood of being flooded. This will assist the planning for future land use zoning, which can be used to restrict further developments in the floodplains.

Figure 13 Probabilistic flood map of the study area. Varying colours indicate probabilities of obtaining flood depths more than 0.5 m.

4.5 Discussion on the results

The flood extent for the simulated result of 72 scenarios was compared. Flood extent varies for different scenarios with different conditions such as the daily (high and low tide) and biweekly (spring and neap tide) tidal variation, sea level rise and the angle of landfall (Fig. 9).

Three worst-case scenarios (Scenario S1, S2 and S3) were compared by generating flood maps and calculating total flooded areas and total damages. The flood maps for Scenario S1, S2 and S3 (Fig. 8) demonstrated that a large area was flooded for all the breach locations. At least 25 % of the total area of Polder 48 was inundated for the three scenarios (Table 3). In the case of all three scenarios considered, the inundation area with flood depths from 0.5 to 1.0 m was larger than the inundation areas with other flood depths (Fig. 10). The inundation area with flood depths more than 1 m was largest for Scenario S3, due to the depressions close to the dikes (Fig. 10). The rice fields were flooded most, while the shrimp ponds were flooded least in all the scenarios (Table 3).

Flood risk was quantified with damage due to floods (negative consequences) and the probability of occurrence. The total estimated damages due to flooding for Scenario S1, S2 and S3 were EUR 10.7, 10.6 and 8.6 million, respectively (Table 4). For all the scenarios, a 1-in-25-year cyclone event was considered. The damages in the settlements were greater than other land classes for all the scenarios (Table 4). Rice fields were flooded most but they did not experience the highest damage compared to other land use classes (Tables 3 and 4). This can be explained by the high damages in settlements compared to rice fields (Table 4). The damage to crops depends on the flood depth, duration and overland flow velocity. For simplification, only the damage related to flood depth was used. As the probability of cyclones was considered the same for all the scenarios, the calculated damage governed the estimated flood risk; i.e., higher damage to the settlements translated as a higher risk of flooding. The primary economic activity of the inhabitants of the study area is farming (BBS, 2011), and most of the inhabitants are poor (with a poverty rate of 0.628) (Alamgir et al., 2018). Even though the estimated damage and risk of flooding to crops were much less compared to areas with other land uses, it will affect the people living in the study area most as they depend on the farming of rice for their livelihood (Nasreen et al., 2013). Hasan et al. (2004) found out that the dependence on fishing (in the sea) by the inhabitants of Polder 48 is increasing due to loss of crops by flood, loss of productivity, lack of jobs and poverty. Fishing in the coastal region of Bangladesh yields lower economic returns, leading to enhanced poverty (Hasan et al., 2004).

The damage was maximum for flood depths of 0.1 to 0.5 m for all the scenarios (Fig. 11). The damage due to inundation less than 0.1 m was small and insignificant. Damage is a function of flood depth, but the unit is per unit area (per m²). Therefore, if the flood extent for higher depths is lower, the damage due to flooding might be lower even though the flood damage increases significantly for inundation more than 0.5 m according to the depth–damage curves developed (Fig. 7).

Generated the PFM indicated that the areas adjacent to the dike facing the seaside have a higher probability of flooding and the rice fields are more prone to flooding (Fig. 13). Moreover, the areas protected by mangrove forest might also be flooded if the unprotected location of the dike is breached (Fig. 13), stressing the importance of proper maintenance of the dike everywhere.

The damage due to flooding was maximum for Scenario S1, which results in a higher risk of flooding for Scenario S1. The total flooded area for settlements of Scenario S1 was lower than Scenario S2 (Table 3), but the estimated damage for settlements of Scenario S1 was more than Scenario S2 (Table 4). This indicates that the settlements in Scenario S1 were exposed to greater flood depth and a higher risk of flooding than Scenario S2. Furthermore, Scenario S1 had similar total damage due to flooding, with a lower flood extent than Scenario S2 (Tables 3 and 4). Considering these facts, Scenario S1 was selected as the worst-case scenario, and breaching in the western part of the sea-facing dike was identified as the critical location for breaching during cyclones.

The scenarios with the effect of climate change (sea level rise) had more damage compared to the scenarios without climate change. Scenarios S1, S2 and S3 were associated with the highest storm surge height. With the same set of conditions without climate change (sea level rise), the storm surge height was 6.52 m PWD (Table 2). The damage corresponding to breaching of eastern, central and western locations of the dike due to the storm surge without climate change impact was 23.3 %, 20.5 % and 21.7 % lower than the damage with the climate change impact respectively. The corresponding values for the flood extent were 30.1 %, 21.67 % and 27.21 % lower than the flood extent areas with the climate change impact respectively.

The probability of occurrence of the storm surge and damage caused by inundation were taken into consideration for the risk calculation. In the case of breaching of the dike, the probability of flooding was considered the same as the probability of occurrence of storm surges. The depicted risk map (Fig. 12) shows the areas adjacent to the dike breach are at higher risk, and the risk reduces as the flood propagates towards the east.

Canals are used as a mode of transportation by the inhabitants of the area. Most of the economic activities and residential areas are near the canals. The risk analysis show that the areas at highest risk are the settlements by the canals (Fig. 12). Therefore, although canals play a crucial role in the economy and social life of the area, they also increase the risk of flooding and probability of higher damage to the adjacent areas.

Land use planning plays an important role in the reduction of vulnerability to disasters (Burby, 1998). Probabilistic flood maps (PFMs) can be used for land use planning (Alfonso et al., 2016). For better understanding of the area at risk of flooding due to the breaching of dikes, PFMs were generated for the study area (Figs. 12, 13). The results of 72 scenarios from the scenario matrix were used for the calculation of PFM. The areas adjacent to the sea dikes had a higher probability of flooding due to the breaching of dikes for both PFMs. The areas inland had a lower probability of flooding. Existing land use indicates that the areas with a lower probability of flooding are mostly rice fields (Figs. 12, 13). Land use zoning and management using the PFM can reduce the vulnerability.

5 Conclusions

A 1-D–2-D coupled model was developed to investigate the inundation pattern inside a polder due to breaching of a dike by cyclone-induced storm surges. Different scenarios were formulated and simulated using a 1-D–2-D coupled model. The results of these simulations were used to calculate the total flooded area and damage due to flooding. Simulated results of three worst-case scenarios, S1, S2 and S3, were compared based on the total flooded area and estimated damage. The comparison led to the identification of the critical location of dike breaching during a cyclone. The flood risk map and probabilistic flood map were generated for the dike breaching during a storm surge using the results of the developed scenarios to identify the areas at higher risk and higher probability of flooding.

Flood inundation for the three worst-case scenarios, S1, S2 and S3, indicated that the maximum flooded area was obtained corresponding to the breaching of the central part of the sea-facing dike. The highest depth was obtained corresponding to Scenario S2 (breaching in the central part). The damage for scenarios S1 (breaching in the western part) and S2 (breaching in the central part) was equal. From these findings it can be concluded that the flood extent, flood depth and damage depended on the breach location. Moreover, the comparison of the flood damage and flood extent led to the identification of Scenario S1 as the worst-case scenario and the western part of the sea-facing dike as the critical location for breaching.

The scenarios considering the effect of climate change (sea level rise) indicated that the flood extent and damage due to flooding will increase with sea level rise.

Flood risk was calculated as the product of the probability of occurrence of a flood event and negative consequences (damage). The generated flood risk maps indicated that for all the scenarios, areas adjacent to the dike and canals inside the polder had a higher risk of flooding. For better access to the canals, for transportation and for livelihood, development of infrastructure and households nearby the canals increases the vulnerability. Similarly, developing land for infrastructure and household on the landside of the dikes increases vulnerability. Combining the effects of increased vulnerability and higher flood depth results in an elevated risk of flooding due to the breaching of dikes during a cyclone.

Inundation maps of all 72 scenarios were compared to generate the probabilistic flood map, which indicated that the areas with rice fields are the least probable areas to be flooded, and the settlements are the most probable areas to be flooded. Although the inhabitants are mostly dependent on agriculture, the flooding of settlements will cause most damage and force relocation.

Measured storm surge levels for previous cyclones were unavailable. Therefore, for this research synthetic water level time series were generated considering the storm surge height presented by Islam et al. (2013), for a cyclone with a 25-year return period. The probability of flooding in a protected area is complicated, and it was assumed that the probabilities of storm surge occurrence and breaching are the same. A limited number of field observations were available to compare the results of the 2-D model. The limited calibration possibility of hydraulic models stresses the importance of field observations before, after and during flood events. As the future land use data were not available, current land use has been used for the future scenarios as well.

The primary objective of the research was to present a methodology for generating FRMs and PFMs for the breaching of dikes during a cyclone. Due to the lack of data on the existing conditions and previous history of breaching of the dike, the probability of the dike breaching could not be determined. Comprehensive surveys should be conducted to determine the physical conditions of the existing embankments and their breach history. Using these data, a joint probability of flooding due to storm surges and breaching may be considered in future studies. As the sea beach outside the dike on the seaside was not included in the 2-D model, the effect of mangrove forest could not be determined. A single breach location was considered for all the scenarios developed. The probability of multiple dike breaching for a polder should be studied as well. Moreover, due to a lack of data, the storm surge height for the present scenarios was used for future scenarios as well. As the sea surface temperature will change in the future due to climate change, the height and intensity of the storm surges will be affected as well. Research on the change of storm surge height and intensity due to climate change should be conducted in the future. Bathymetric data with a coarse grid resolution from GEBCO were used as the measured bathymetric data for the sea were not available. Furthermore, the study relied on the previous literatures for developing depth–damage curves. Conducting field surveys to generate these curves will provide a more reliable estimate of damages due to flooding. The developed and simulated model depended on the field measurements and logical assumptions which might be the source of errors. For damages, only direct damages were included. Inclusion of indirect damages will provide more realistic estimates.

Bangladesh is a hazard-prone country, and cyclone-induced storm surges are one of many natural disasters that affect the coast of Bangladesh. The storm surges cause severe damage to the earthen embankments/dikes protecting the coastal polders. The methodology presented in this paper to develop the 1-D–2-D inundation model, the PFMs and the risk maps and to identify the critical locations for breaching can assist in better preparedness against flooding and help in damage reduction through land use zoning and management. At present, the PFM and FRM due to storm surges and breaching of the dikes are not available for the coastal polders.

Climate change will likely cause increase in the frequency and intensity of cyclones around the world. This will call for large investments for the improvement of existing or new protection structures for the deltas. The identification and prioritization of maintenance of critical locations of dike breaching can potentially prevent a disaster. Non-structural tools such as land use zoning, with the help of flood risk maps and probabilistic flood maps, have the potential to reduce the risk and the damage due to dike breaching. The method presented in this research can potentially be utilized for the deltas around the world to reduce vulnerability and flood risk due to the breaching of dikes caused by cyclone-induced storm surges.