Modeling complex flow dynamics of fluvial floods exacerbated by sea level rise in the Ganges–Brahmaputra–Meghna Delta

Global warming is likely to exacerbate future fluvial floods in the world’s mega-delta regions due to both changing climate and rising sea levels. However, the effects of sea level rise (SLR) on fluvial floods in such regions have not been taken into account in current global assessments of future flood risk, due to the difficulties in modeling channel bifurcation and the backwater effect. We used a state-of-the-art global river routing model to demonstrate how these complexities contribute to future flood hazard associated with changing climate and SLR in the world’s largest mega-delta region, the Ganges-Brahmaputra-Meghna Delta. The model demonstrated that flood water in the main channels flows into tributaries through bifurcation channels, which resulted in an increase in inundation depth in deltaic regions. We found that there were large areas that experienced an increase in inundation depth and period not directly from the SLR itself but from the backwater effect of SLR, and the effect propagated upstream to locations far from the river mouth. Projections under future climate scenarios as well as SLR indicated that exposure to fluvial floods will increase in the last part of the 21st century, and both SLR and channel bifurcation make meaningful contributions.


Introduction
The world's mega-delta regions are vulnerable to floods because they lie in a low and flat topography, and experience both overflow from rivers and storm surges [1,2]. In addition, mega-delta regions are major centers of population (about half a billion people [3]) and agriculture that hold both ecological and economical value [4]. Floods cause adverse effects such as heavy casualties and loss of assets, restricting industrial development of countries in mega-delta regions. Bangladesh contains the world's largest mega-delta region, which is formed by the Ganges-Brahmaputra-Meghna river systems (the GBM Delta). Severe flooding events have frequently occurred in the country (e.g. , 1987, 1988, 1998, and 2007) due to heavy and intensive rainfall in monsoon season over the GBM river basins [5]. For example, annual flood losses in 2007 were approximately 1.1 billion United States Dollars (USD) [6].
Previous studies have predicted an increase in the frequency of fluvial floods associated with changing climate globally throughout the 21st century [7][8][9] including in mega-delta regions. Moreover, the global sea level is projected to rise during this century and this will continue beyond 2100 [10]. Deltaic lowlands will become more vulnerable to inundation, flood events, and erosion as the sea level rises (SLR) [3,11,12]. Given that fluvial floods and SLR are common issues among deltaic regions worldwide, methods applicable to global-scale analyses are preferable. Therefore, to simulate land-coast interactions including both changing climate and SLR has considerable implications for future projections of fluvial floods and their associated hazards in the world's mega-delta regions. However, to the best of our knowledge, no global-scale assessments of future flood risk have simulated the effects of SLR on fluvial floods in mega-delta regions. Some studies have globally projected changes in fluvial floods under warming climate by 2100 [8,13], but they did not take into account SLR (hereafter we focus on fluvial floods and omit the word 'fluvial' in indicating fluvial floods). Some regional-scale hydrodynamic models have been used to show the effect of complex river networks [14] or the impact of SLR on storm surge floods [15], but they have not yet been applied to global-scale analysis due to the considerable amount of boundary and topography data processing required for global simulations.
Complex water flow in deltaic lowlands is a major difficulty in flood modeling in mega-delta regions because water flow often bifurcates and changes its direction. To simulate future changes in fluvial floods as well as the effects of SLR, it is indispensable to utilize basin-scale hydrologic modeling together with the complex flow dynamics of the lowland area. In particular, a continental-or global-scale river routing model is required to project flood risk for a mega-delta such as the GBM Delta, whose catchment area is 1.6 million km 2 [16]. However, since previous global river routing models have assumed only one downstream flow direction for each grid cell [17,18], bifurcating water flow in mega-delta regions cannot be reproduced using those models.
The representation of the backwater effect adds further complexity to flood projections with SLR. The backwater effect is the propagation of the effect of SLR upstream and can be expressed as changes in water levels in the river flow calculation. While some regional-scale inundation models use a hydrodynamic flow equation containing the backwater effect [14], most global river routing models (e.g., PCRGLOB-WB Dyn-Rout [16], LISFLOOD [19]) use a kinematic wave approximation of the momentum equation, which prevents an explicit representation of the backwater effect from SLR in global-scale analyses.
In this study, we analyzed the effects of channel bifurcation and SLR on fluvial floods in the world's largest mega-delta region, the GBM Delta. Yamazaki et al recently developed a method to represent channel bifurcation in mega-delta regions within a framework of global hydrodynamic models [20]. Applying this model to our study, we overcame two complexities mentioned above: complex water flow through divergent channels and the backwater effect from SLR. In addition, the future impacts of a changing climate including SLR on fluvial floods in Bangladesh were estimated as flood exposure using atmosphere-ocean general circulation models (AOGCMs).

Global river routing model: CaMa-Flood
A state-of-the-art global river routing model (Catchment-based Macro-scale Floodplain Model: CaMa-Flood [20][21][22]) was used to integrate runoff input along a river network and to calculate river discharge, flood inundation extent, and inundation depth at a 0.1°spatial resolution (approximately 11 km at the equator). The modeled inundation depth was diagnostically downscaled onto a 9 arc-second spatial resolution (approximately 250 m at the equator) digital elevation model (DEM) using a similar method to Winsemius et al [16] (for detailed method, see supplementary information S1). The latest version of the model implements channel bifurcation to achieve a reasonable reproduction of the complex flow dynamics in deltaic regions, and demonstrated good performance in a historical simulation of the Mekong River basin [20]. Bifurcation channels represent downstream river flow in multiple directions (red lines in figure 1). In this study, bifurcation channels were automatically delineated following the same method as Yamazaki et al [20] (see more detail in this reference and supplementary S1 and S2).
A hydrodynamic flow equation (i.e., the local inertial equation [22]) was used as a momentum equation of CaMa-Flood instead of the traditional kinematic wave equation, so that the model facilitated flood simulation with backwater effects caused by SLR. Runoff input to CaMa-Flood was obtained from the offline simulation of a land surface model, the MATSIRO [23] with explicit groundwater representation (MAT-SIRO-GW [24]), driven by climate forcing [25] at a sub-daily time step over the period of 1979-2010 and at spatial resolution of 1°over the domain (73-98°E, 21-32°N). The modeled river discharge was validated against observations. River discharge was reasonably consistent with observations (correlation coefficient was 0.87 and 0.90 for the Ganges at HB and the Brahmaputra at BD, respectively. Locations are shown in figure 1). The variation in river water level was also well reproduced (correlation coefficient was 0.87 at BB, 0.84 at GL, and 0.69 at CD). Details of the validation are presented in supplementary S2.

Experiment design
Using CaMa-Flood, we conducted a series of simulations to determine how channel bifurcation and SLR affect flood inundation in the GBM Delta. First, to investigate how bifurcation channels influence the simulation of floodplain inundation, we performed simulations with and without channel bifurcation (Bif and NoBif, respectively) assuming no SLR. Second, to assess the influence of the backwater effect of SLR on flooding, we conducted simulations with and without the backwater effect caused by SLR (BWSL and NoBWSL, respectively. Schematics of these simulations are shown in figure 3(a)). In the BWSL simulation, SLR was represented as a change in the downstream water level boundary condition at the river mouth in river discharge calculation. In the NoBWSL simulation, we first conducted a flood simulation in which the sea level is 0 m, and afterwards raised the sea level to the height of the assumed sea level increase. NoBWSL simulation ignored the backwater effect caused by SLR and just indicates land submergence due to SLR. Note that both BWSL and NoBWSL simulations include the backwater effect in river discharge calculation. To highlight the impact of SLR, we conducted sensitivity experiments for past floods where different SLR experiments (1 m and 2 m) were applied. Because the target period of flood exposure estimation (section 4) extends to 2100, we chose a 1 m SLR scenario in reference to IPCC. This value is the approximate number of the worst scenario within the range of RCP 8.5 (0.98 m) [10]. In addition, to make the physical impact clear, we added 2 m SLR to the scenario, which is one extreme case of SLR projections of previous studies cited in IPCC [26].
Here we neglected the dynamic changes in coastal topography due to gradual changes in SLR. The 2007 flood was used for both channel bifurcation and SLR analyses, because it was one of the heaviest floods in recent decades in Bangladesh [6]. Figure 2(a) shows the modeled annual maximum inundation depth in 2007 [5]. The results show inundation areas in and around the main channels of the Ganges, Brahmaputra, and Meghna, which can be seen in the satellite-derived (MODIS) inundation area (figure 2(c). See also supplementary information S2 and figures S3-S5) [27]. When bifurcation channels were incorporated into the model, inundation depth decreased in the main channels and increased in the tributaries (figure 2(b)). A significant difference was seen in the midstream of the Ganges and the upstream Meghna, where the bifurcation channels were delineated ( figure 1(a)). The maximum increase and decrease in inundation depth was 9.2 m and 3.9 m, respectively. A wide area of the lowland deltaic region was inundated via bifurcation channels, which did not occur in the original river model (NoBif).

Effects of channel bifurcation
The large difference in inundation depth between the Bif and NoBif simulations was because floodwater in the main channels flowed into tributaries or the deltaic region through bifurcation channels. Even in large flooding events, if the bifurcation channels are not taken into account, the water level of the tributaries and the deltaic region, which are not connected to the main channels, remains low. Furthermore, the areas where the inundation depth increased due to the bifurcation channels were located in accordance with the flood-prone areas of satellite-derived inundation map ( figure 2(c)). Hence, the channel bifurcation scheme is a reasonable approach because it expresses the dynamic floodwater allocation between adjacent independent sub-basins.  Figures 3 and 4 show how the backwater effect caused by SLR would affect the flood hazard in the GBM Delta. In BWSL simulation, the impact of the SLR spreads upstream due to the backwater effect ( figure 3). For example, the inundation depth at the river mouth in the Ganges increased by 1.4 m. In contrast, the increase in inundation depth in the NoBWSL simulation was limited to the coastal region. These characteristics were also found in the 1 m SLR experiment (supplementary figure S6). Figure 4 shows observed and modeled water surface elevation (elevation of a water surface in a river channel) with the backwater effect (BWSL simulation). At BB and CD seasonal variation in water surface elevation without SLR is similar to that of observations, however, the rise in water surface elevation is delayed in the simulations. This may be due to an underestimation of snow melt runoff calculated by a land surface model. The large difference in water surface elevation at GL may be due to parameterization of river channel depth by an empirical equation (parameterization of river channel depth and validation of water surface elevation are available in supplementary S2). Despite the limitations of this model simulation (supplementary S5), we assume that the overall performance of the model is adequate to quantify the effects of complex water flow in deltaic lowlands under SLR conditions.

Effects of backwater caused by SLR
The changes in water surface elevation due to SLR via the backwater effect varied among three different locations ( figure 4). The site near the river mouth (Chandpur: CD, modeled elevation is 0 m) was directly affected by the SLR, with the increase in the water surface elevation being the largest of the three locations. At all three locations, the water surface  elevation increased by the same amount as the difference between each elevation and SLR in the dry season. However, in the wet season, the increase in the water surface elevation was greater downstream (CD) and became smaller upstream (BB and GL). Near the river mouth (CD), the difference in the annual maximum water surface elevation between the 2 m and 0 m SLR experiments was 1.3 m, whereas in the confluence of the Ganges and Brahmaputra (Goalondo: GL, modeled elevation is 1 m) it was 0.90 m. In the upstream point of the Meghna (Bhairab Bazar: BB, modeled elevation is 1 m) the difference was 0.49 m.
The SLR in the BWSL simulation was represented by a change in the downstream water level boundary condition at the river mouth. The change in the river discharge at the river mouth was then propagated upstream through the backwater effect, causing a large increase in water surface elevation as well as a longer duration of high water surface elevation in the BWSL simulation (shown as blue and red areas in figure 4). On the other hand, water surface elevation in the dry season was close to sea level even upstream due to the backwater effect, because in CaMa-Flood, the elevation of a grid cell is defined as that at the top of a river channel reservoir and channel depth is determined by upstream river discharge; hence, the elevation of a river bed is often below sea level in low and flat regions such as the GBM Delta. Therefore, water surface elevation approaches the sea level even in the dry season. This result implies that when the sea level rises, water surface elevation upstream could rise and may cause issues such as poor drainage. These findings could not be represented in the NoBWSL simulation due to the lack of a backwater effect. Hence, an adequate hydrodynamic flow equation (the local inertial equation [22]) is required for an analysis of the backwater effect, instead of the kinematic wave equation traditionally employed in previous river routing models.

Application: future flood exposure in Bangladesh
The results reported here can be transferred to the implication of future risk of flooding. To quantify future changes in flood risk under changing climate, we computed the number of people (flood exposure) in Bangladesh affected by flooding. In addition, the contributions of channel bifurcation and SLR to flood exposure calculations were analyzed.

Method
Past and future CaMa-Flood simulations, with channel bifurcation and the backwater effect, forced by daily runoff outputs calculated from AOGCMs were used to compute flood exposure. The runoff projections of 11 AOGCMs participating in the Coupled Model Intercomparison Project (CMIP5 [28]) were used to obtain the potential range of future runoff projections under two Representative Concentration Pathways (RCPs). These AOGCMs were selected due to the availability of daily runoff data and because their host institutions were independent, since different versions of AOGCMs from the same institution may not be considered independent (details on the AOGCMs selected in this study are available in supplementary table S4).
We conducted historical AOGCM simulations (hereafter called 'Historical' run, 1960-2005) and future simulations (2006-2099) forced by RCP 4.5 and 8.5. They were interpolated from the original spatial resolution of a land surface model implemented in each AOGCM, and then integrated to discharge on a river network map (spatial resolution is 0.1°). As AOGCMs are prone to biases in precipitation including weak extreme precipitation, the absolute magnitude of runoff data from the AOGCMs may not be able to replicate a reasonable inundation depth under flooding. Hence, the input runoff data from the AOGCMs were adjusted using the long-term (1986-1995) mean of the ensemble of runoff reanalysis by GSWP2 [29]. The daily runoff of AOGCMs were adjusted using a simple ratio method such that the long-term (1960-2005) monthly means of AOGCM runoff are the same as those of GSWP2. For future SLR, 0 m and 1 m scenarios in 2100 were selected, which included the worst scenario in the IPCC AR5 [10].
To compute flood exposure, we first downscaled the modeled inundation depth to 9 arc-second spatial resolution and produced time series of annual maximum inundation depth, showing the potential extent of maximum inundation area. The flood exposure was then estimated by overlaying the modeled inundation area onto a gridded annual population density (derived from the History Database of the Global Environment [30]). The spatial resolution of population data was 5 arc-minutes and it was uniformly disaggregated to 9 arc-second to overlay it onto the modeled inundation depth. In this process, population in the areas inundated over the entire period of 1979-2010 was not accounted to exclude river channel pixels, because CaMa-Flood assumes that each grid cell has a river channel. For all simulations the population was fixed at the 2010 level to highlight the changes due to warming climate and SLR. Note that the estimated flood exposure did not take into consideration any flood-protection measures such as dams and levees.  5(b)). Despite the large interannual variability, the flood exposure in RCP8.5 was significantly higher than that in RCP4.5 after the middle of the 21st century. SLR of 1 m contributed an additional increase of 3.48%-4.86%. The contributions of bifurcation channels (NoBif) or backwater effects caused by SLR (NoBWSL) on this increase were 5.39%-14.9% and 2.61%-3.60%, respectively. These results reveal that these effects play significant roles in flood risk projections in deltaic regions.

Conclusions
In this study, we determined the impact of SLR on fluvial flood hazards in the GBM Delta, through the explicit representation of the complex water flow due to bifurcating rivers and backwater effects caused by the SLR. It was apparent that flood water in the main channels flows into tributaries through bifurcation channels, which resulted in increases in the inundation depth (>9 m) in low-lying deltaic regions. There were also large areas that experienced an increase in inundation depth and duration not from the SLR itself, but from backwater effects of SLR, and the effect propagated upstream. These two schemes (channel bifurcation and an adequate hydrodynamic equation to represent backwater effects) are required for flood simulations in mega-delta regions.
Estimation of future flood risk indicated that SLR will increase the flood exposure in Bangladesh (about 5%) in addition to the runoff increase associated with a warming climate (1.5-2 times). Without bifurcation channels or backwater effects from SLR, the calculated flood exposure could be underestimated by about 15% and 3.6%, respectively. The results emphasize that the explicit representation of these effects is not negligible for flood risk projections in deltaic regions.
The proposed method with our inundation scheme has the potential to facilitate new scientific advances. First, while our study focused on the GBM Delta, the method that we employed is applicable to any deltaic region worldwide. In addition, while we focused only on static SLR, it can also be applied to dynamic coastal tides and storm surges. Given that recent studies have indicated the possibility of dependence between fluvial floods and coastal storm surges [31,32], global-scale analysis compounding these two events is a possible avenue for further research. Finally, geomorphological changes in deltaic regions due to sediment transport and coastal erosion should be taken into account to achieve a more comprehensive assessment of flood disasters. Critical future studies include simultaneous simulation of future flood inundation and topographical changes along with the rising sea levels with river routing models, sediment transport models, and GCMs.