Impact Assessment of a Major River Basin in Bangladesh on Storm Surge Simulation

Abstract

A two-dimensional bay and river coupled numerical model in Cartesian coordinates was developed to find the impact of the river on the simulated water levels associated with a storm along the coast of Bangladesh. The shallow water models developed for both the bay and river were discretized by the finite difference method with forward in time and central in space. The boundaries for the coast and islands were approximated through proper stair steps representation and solved by a conditionally stable semi-implicit manner on a staggered Arakawa C-grid. A one-way nested scheme technique was used in the bay model to include coastal complexities as well as to save computational costs. A stable tidal condition was made by forcing the sea levels with the most energetic tidal constituent, M₂, along with the southern open boundary of the bay model omitting wind stress. The developed model was then applied to foresee the sea-surface elevation associated with the catastrophic cyclone of 1991 and cyclone MORA. A comparative study of the water levels associated with a storm was made through model simulations with and without the inclusion of the river system. We found that the surge height in the bay-river junction area decreased by 20% and the surge height reduced by about 3–8% outside the junction area from this study. The obtained results were found to have a satisfactory similarity with some of the observed data.

1. Introduction

A low-lying riverine country, Bangladesh is affected by various environmental disasters like tropical cyclones. In 2013, the World Bank Group reported that Bangladesh would be the most affected country due to global warming. With climate change, a strong rise in the most horrific storms is evident. The concern is that Bangladesh will face various catastrophic cyclones every year due to global warming [1]. According to the Bangladesh Meteorological Department, in the last 44 years, 269 depressions and storms formed during January 1974 to May 2018. Due to this reason, huge economic losses and deaths occurred and this damage is not entirely due to the storm, the coastal structure of Bangladesh and its geological factors are also responsible. A huge number of people died as a result of past cyclone 1970 Bhola (167,000), cyclone 1991 (138,000) and cyclone Nargis 2008 (140,000) within this area. People believe that the storms of Bangladesh are a cruel cannibalistic monster. However, now, due to the extensive change in the study of storm surge research, there has been a decrease in the risk of death. Most of the research on storm surge has been conducted in the Atlantic, North Sea, Pacific, Gulf of Mexico, Gulf of Carpentaria and Bay of Bengal. For the Atlantic region and some other regions, some oceanographers [2,3,4,5,6] have conducted potential and pioneering works. The coastal zone of Apalache Bay and the northeastern Gulf of Mexico are at risk of acute severe storms. Some notable work has documented the storm surge in these coastal marshes [7,8]. Some 2D and 3D surge simulation models have been applied in this region to accurately investigate the surge height. The well-known surge simulation models SLOSH (Sea, Lake and Overland Surges from Hurricane) [9], CH3D [10,11], POM (Princeton Ocean Model) [12,13], FVCOM (Finite Volume Coastal Ocean Model) [14,15] and ADCIRC [16] have been used for this coast. For the SLOSH and POM model, a curvilinear orthogonal grid is used for the numerical computation. A non-orthogonal curvilinear grid is used in the CH3D model. After developing the grid system, an unstructured triangular grid is used for the ADCIRC and FVCOM models. In our study, a rectangular grid system was used in the bay-river coupled numerical model to compute the water level elevation due to a storm with a reduction in the time and cost. According to Alam [17], there are 6.69 cyclones and depressions crossing the coast of Bangladesh per year, which is 86.14% of the total formation of storms and depressions. As a result, floods occur around the coastal areas, which is also a major disaster for Bangladesh. Some numerical models for storm surge simulation in coastal areas have shown a little discrepancy. The reasons behind this may be the coastal complexities, wherein some worth mentioning can be pointed out as shallow bathymetry, long continental shelf, existence of offshore islands of different shapes, sediment accretion, curvilinearities of the coastal and island boundaries, high tidal range with the highest along Sandwip channel (Figure 1) and the high river discharge of the combined flow of Ganges–Brahmaputra–Meghna through lower Meghna and so forth. [18,19]. Besides these facts above-mentioned some other factors, namely, river runoff, torrential rain and so forth, also contribute to the rise in sea level caused by storm surges. The other factor, the north–west direction movement of cyclones striking the Bangladesh and Indian coast create serious disasters every year. The study of cyclone modeling in the Indian seas has been confined to the Bay of Bengal, especially the Coast of Bangladesh [20,21].

Thus, a complete understanding of the factors affecting surge development and its accurate prediction along the coast of Bangladesh is highly desirable. Although a considerable number of studies [18,20,21,22,23,24,25] have been conducted for the coast of Bangladesh, these studies have been insufficient at accurately incorporating the effect of the river system in storm surge modeling. The authors in [20] were the first to incorporate the effect of an idealized lower Meghna in a storm surge simulation, finding a contribution of the river on storm surge simulation. As the river was considered as an idealized river system in their study, their result may be somewhat different from reality. The work was further developed by Johns [26] to resolve the surge development along the coastal belt of Orissa in India with the use of a continuous contracting grid system. Dube et al. [27] conducted their study of the Bay of Bengal, adding the effect of the lower Meghna River for the simulation of storm surge along the Bangladesh coast. In their study, they also assumed an idealized river system and used a fixed water depth, which is also unrealistic. In subsequent work by Dube [28], the study coupled a one-dimensional river model with two-dimensional bay model with a curvilinear representation of the natural shoreline for the Orissa coast in India. In their study, a tidal solution was not used to provide the initial conditions for tide and surge interaction in the bay. Similarly retaining the limitations of the study of Agnihotri [29] introduced a two-dimensional bay-river coupled model for the simulation of storm surges along the east coast of India. The study was also limited by the lack of a representation of proper geometry of the river system.

The present study is a step towards the development of an operational surge forecasting bay and river-coupled model that incorporates the river system with realistic geometry. A non-linear interaction of tide and surge was also incorporated in this study. We also intended to investigate how the river influences the sea surface elevation when we incorporated the river into our model simulation as a channel gate.

2. Study Area and Cyclones

2.1. Study Area

One of the countries with the most risk of cyclones, Bangladesh, is situated at the northern tip of the Bay of Bengal between longitudes 88° E to 92° E and latitudes 20° N to 26° N.

Nevertheless, the considered area was between the longitudes 85° E to 95° E and latitudes 15° N to 23° N to investigate a storm lasting two or three days in the analysis area before crossing the coast. The study area is bordered on the west, north and east by India, on the southeast by Myanmar and on the south by the Bay of Bengal (see Figure 1). It has long been advocated that the data-poor region, the coast of Bangladesh is the world’s most vulnerable place but has not been sufficiently studied. The reason behind this vulnerability may be the thickly populated offshore islands of different shapes, the curvilinearity of the coastal and island boundaries, long continental shelf, shallow bathymetry, high tidal range with the maximum at Sandwip channel (Figure 1) near Chittagong (about 5–6 m), huge freshwater discharge through the lower Meghna and so forth. Since the Meghna River contributes mainly to surge development in the area of interest, it was also considered in this study. In the river model, the lower part of the Meghna River was considered from mouth (23° N) to the head (23.25° N), wherein at the junction position of the bay and river (mouth of the river), the breadth was considered from longitude 90.40° E to longitude 90.68° E (Figure 2). The area considered was extended from Haiderganj to Moishadi (Chadpur), which is about 40 km in length (red color marked). It is of interest to note here that the sediment discharge from the river Meghna (combined flow of the Ganges–Brahmaputra–Meghna) is the highest and the freshwater discharge is the third highest of all the river systems worldwide [30].

2.2. Cyclones

Bangladesh faces different types of natural cyclones every year. In terms of deaths and economic losses, the cyclones Bhola (1970), BOB1 (1991), SIDR (2007), AILA (2007), Nargis (2008) and MORA (2017) have been the significant cyclones in the last few decades.

In our study, we chose two cyclones for investigation: cyclone 1991 and cyclone MORA.

The 1991 cyclone has been studied enough by some oceanographers and it had a catastrophic nature. The other cyclone is MORA, which struck the Bangladesh coast recently. These cyclones should be investigated properly as one storm proved dangerous for causing human death and the other for economic losses. Furthermore, the parametric characteristics of the storms were different such as wind speed, maximum radius, coast crossing angle and landfall location. Figure 3 shows the path of the cyclone with the time and location that we used in our study.

2.2.1. Cyclone 1991 (BoB 01)

The 1991 cyclone of Bangladesh, known as BOB 01, was the most catastrophic and deadliest tropical cyclone that struck near Chittagong on the southeast side of Bangladesh on 29 April 1991. The approximate losses during the cyclone were 1.72 billion US dollars and the deaths of 138,000 people. The first circulation was formed on 22 April 1991 near the equator in the Indian Ocean and most of the mass cloud surrounded the Bay of Bengal region. The Indian Meteorological Department (IMD) detected the system as a depression on 23 April and it was also detected by the satellite picture taken from NOAA-I1 and GMS-4 satellites [31]. After that, the storm then moved northeastward from its genesis position. From the latitude 11.80° N to 13° N, the storm was a cyclonic storm, having a central pressure of 996 hPa. By this time, on 27 April, the system reserved its intensity until 28 April (UTC 6.00) and changed its nature as a severe cyclonic storm. After that, the system then forwarded to accelerate to the north–northeast side, which made it further intensify with a 938 hPa central pressure. At that time, the astronomical position of the cyclone was 20.60° N latitude and 90.70° E longitude. Finally, the storm crossed the coastal area of Bangladesh with a 6-m surge height on 30 April at 2:00 a.m. local time.

2.2.2. Cyclone MORA

According to the official announcement, on 28 May 2017 (midnight), when the cyclone was known as MORA, the cyclone was located 670 km away from Chittagong Port, 590 km away from the Cox’s Bazar Port, 715 km away from Mongla Port and 650 km away from Payra Port. Before this status, the depression was first observed on the afternoon of 27 May in the central Bay of Bengal. Later on, it became more condensed and crossed Chittagong Port at 6:00 a.m. local time on 30 May 2017. At that time, the wind speed and central pressure was 110 km/h and 978 hPa, respectively. However, the Ministry of Disaster Management and Relief of Bangladesh constituted a strong committee before the storm struck the Bangladesh coast and advised the local council to take necessary action against the cyclone disaster. With adequate alerts and the help of the Government, the number of deaths decreased but the number of economic losses increased greatly. According to the Ministry of Disaster Management and Relief, the death toll consisted of nine people and the economic loss was estimated to be 1.7 billion US dollars. After the storm hit the coast, it became a land deep depression and crossed the Rangamati district and surrounding area. At this time, heavy rainfall occurred in the region, causing floods. According to the Bangladesh Meteorological Department, rainfall was recorded as 225.2 mm in Chittagong, 213 mm in Sandwip, 208 mm in Sitakunda and 187 mm in Rangamati. The detailed statistical information of cyclone MORA is presented in Appendix A.

3. Numerical Model

3.1. Data and Materials

In our numerical model, we used the meteorological, hydrological, geographical and oceanographic data involved with some of the parameters. Meteorological input data like storm track, central pressure, maximum sustained wind radius and maximum sustained wind speed were taken from the Bangladesh Meteorological Department (BMD). The wind distribution data were generated by the empirically based formula of [2] from information from the Bangladesh Meteorological Department. The water depth data were collected from the general bathymetry chart of the ocean (GEBCO) 30 arc-second interval grid data and the data were interpolated to the model grid by using the inverse distance weighted interpolation method. Coastal geometry and island geometry of our study domain were collected from the geophysical data system (GEODAS) coastline shapefile and the shapefile was processed by using coastline extractor software and the MATLAB subroutine. Hydrological data like river discharge Q = 5100.00 m³/s (yearly average) and tidal elevation were collected from the study of [32] and the Bangladesh Inland Water Transport Authority (BIWTA), respectively. The friction coefficient C_f = 0.0026 and the drag coefficient C_d = 0.0028 were taken as uniform from the study of [22]. The drag coefficient could be calculated from the relation $C_d={\rho }_{w}g{n}^2/H^{1/3}$, where ${\rho }_w$ is the seawater density, g is the gravitational acceleration, n is the Manning roughness coefficient, and $H$ represents the water depth. The friction coefficient can be calculated from the relation $C_f=\tau /D_p$, where $\tau$ represents the shear stresses and $D_p$ represents the dynamic pressure. The cyclone track data were 3-h interval data that were collected from the Bangladesh Meteorological Department observed data. Figure 3 represents the track information of cyclones 1991 and MORA.

3.2. Model Description

3.2.1. Parent Model (Bay Model)

A system of rectangular Cartesian coordinates was used, where the curvature of the earth’s surface was assumed to be zero. In this system, the origin, O, is taken in the undisturbed level of the sea surface, OX points towards the south, OY points towards the east and OZ is directed vertically upwards. The displaced position of the free surface is given by $\zeta (x,y,t)$ and the position of the seabed by $h(x,y)$. For any atmospheric or oceanic phenomenon, if the horizontal length scale is large in comparison with the vertical length scale, then the z-component of the momentum equation can be approximated by the hydrostatic equation [19]. By approximating and neglecting the molecular velocity, the flux form of the vertically integrated shallow water equation for the dynamical process in the sea can be written as

$$\frac{\partial \zeta }{\partial t}+\frac{\partial \tilde{u}}{\partial x}+\frac{\partial \tilde{v}}{\partial y}=0,$$

(1)

$$\frac{\partial \tilde{u}}{\partial t}+\frac{\partial (u\tilde{u})}{\partial x}+\frac{\partial (v\tilde{u})}{\partial y}-f\tilde{v}=-g(\zeta +h)\frac{\partial \zeta }{\partial x}+\frac{T_x}{\rho }-\frac{C_f \tilde{u}\sqrt{(u^2+v^2)}}{\zeta +h},$$

(2)

$$\frac{\partial \tilde{v}}{\partial t}+\frac{\partial (u\tilde{v})}{\partial x}+\frac{\partial (v\tilde{v})}{\partial y}+f\tilde{u}=-g(\zeta +h)\frac{\partial \zeta }{\partial y}+\frac{T_y}{\rho }-\frac{C_f \tilde{v}\sqrt{(u^2+v^2)}}{\zeta +h},$$

(3)

where $(\tilde{u},\tilde{v})=(\zeta +h)(u,v)$ and $u,v,$ are the Reynolds averaged components of velocities in the directions of $x$ any $y$, respectively; $\rho$ is the density of the sea water supposed to be homogeneous of the same density; $f=2\mathsf{\Omega }\sin \phi$ is the Coriolis parameter where $\mathsf{\Omega }$ is the angular speed of the earth rotation and $\phi$ is the latitude of the place of interest; and $g$ is the acceleration due to gravity.

In this study, the wind stress components were parameterized by a conventional quadratic law [21] and is given by

$$(T_x,T_y)={\rho }_a{C}_d{(u_a^2+v_a^2)}^{\frac{1}{2}}(u_a,v_a),$$

(4)

where C_d = 2.8 × 10^‒3 is the surface drag coefficient used by [22]; ${\rho }_a$ is the density of the air; and $u_a$, $v_a$ are the $x$ and $y$ components, respectively, of the wind stress. The circulatory wind field based on the available information at the Bangladesh Meteorological Department (BMD) can then be generated by the empirical formulae of [2], being given by

$$V_a=\{\begin{array}{l}{V}_0\sqrt{{\left(r_a/R\right)}^3}\mathrm{for} \mathrm{all} r_a\le R\\ V_0\sqrt{\left(R/r_a\right)}\mathrm{for} \mathrm{all} r_a>R,\end{array}$$

where $V_0$ is the maximum sustained wind at the radial distance $R$ and $r_a$ is the radial distance at which the wind field is desired.

3.2.2. River Model

A separate two-dimensional depth averaged river model was also developed to incorporate the effect of the river from the study of the researcher [29], wherein the modeling, the curvature of the earth was also assumed to be zero and the curvilinear properties of the study region for the river were approximated through proper stair steps. It was previously mentioned that the model domain was considered from 23° N to 23.25° N latitude and 90.40° E to 90.68° E longitude, which turned out to be a length of about 40 km and a breadth of about 16 km. For the developed river model, the origin O′ was taken at the equilibrium level of the free surface and was located at the left corner of the head of the river latitude (23.25° N, 90.40° E). The $x$-axis was taken as positive in the stream direction while the $y$-axis was taken parallel to the coast. The elevation of the free surface above its mean level was denoted by $z={\xi }_{riv}(x,y)$ and the bottom of the river by $z=-h_{riv}(x,y)$. The depth averaged velocity components in the developed river system, $u_{riv}$ and $v_{riv}$, satisfied the following equations

$$\frac{\partial {\xi }_{riv}}{\partial t}+\frac{\partial }{\partial x}\left[H{u}_{riv}\right]+\frac{\partial }{\partial y}\left[H{v}_{riv}\right]=0,$$

(5)

$$\frac{\partial }{\partial t}\left[H{u}_{riv}\right]+\frac{\partial }{\partial x}\left[H{u}_{riv}{}^2\right]+\frac{\partial }{\partial y}\left[H{u}_{riv}{v}_{riv}\right]-H{v}_{riv}=-g\frac{\partial {\xi }_{riv}}{\partial x}-C_f{u}_{riv}{\left(u_{riv}{}^2+v_{riv}{}^2\right)}^{1/2},$$

(6)

$$\frac{\partial }{\partial t}\left[H{v}_{riv}\right]+\frac{\partial }{\partial x}\left[H{u}_{riv}{v}_{riv}\right]+\frac{\partial }{\partial y}\left[H{v}_{riv}{}^2\right]+H{u}_{riv}=-g\frac{\partial {\xi }_{riv}}{\partial y}-C_f{v}_{riv}{\left(u_{riv}{}^2+v_{riv}{}^2\right)}^{1/2},$$

(7)

where $H=[{\xi }_{riv}+h_{riv}](x,y)$ represents the total water column depth where $u_{riv},v_{riv},$ are the velocity component. The bottom stress, following the study [29], was parameterized as $\rho C_f\sqrt{u_{riv}^2+v_{riv}^2}$. To understand the discretization (see Appendix B).

3.3. Boundary Conditions

3.3.1. Boundary Condition of the Parent Model

The normal component of the depth-averaged velocity for open boundaries can be given by various formulas depending on the available data [33]. For the Bay of Bengal region, most of the researchers have used a radiation type of boundary condition, which is able to allow the disturbance generated within the model area to go out through the open boundary as a progressive wave. The following boundary condition was taken from [21], the western (85° E), eastern (95° E) and southern (15° N) are, respectively, given by

The east boundary (85° E):

$$v+{(\frac{g}{h})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\zeta =0$$

(8)

The west boundary (85° E):

$$v-{(\frac{g}{h})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\zeta =0$$

(9)

The south boundary (15° N):

$$u-{(\frac{g}{h})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\zeta =-2{\left(\frac{g}{h}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}aSin\left(\frac{2\pi t}{T}+\phi \right)$$

(10)

For the closed boundaries of the mainland as well as of the islands, the normal components of the mean current were taken as zero. In Equation (7), $a$, $\phi$ and T denote the amplitude, phase and period of the tidal constituent, respectively, under consideration.

3.3.2. Boundary Condition for the River Model

For the eastern and western river boundaries, whenever it is a vertical sidewall, the boundary condition can be set as [25]

$${\xi }_{riv}=0$$

(11)

The extended analysis area of the northern boundary was situated at the head of the river, where the conditions can be put into the following form. If the fresh water discharge from the river is taken into account, then the appropriate boundary condition can be expressed as [25]

$$u_{riv}=\frac{Q}{B_{riv}{h}_{riv}},\text{ and }{u}_{riv}=0$$

(12)

when no fresh water discharge is considered, where $B_{riv}$ is the river breath and Q represents the river discharge.

3.3.3. Matching Boundary Condition

The southern boundary of the river model is influenced by the bay, which is why running a model with a perfect match between the river and bay is important. The continuity of volume flux should be ensured at the mouth of the river.

The previous study [20] explains one dimensional matching equation. By improving their idea, after simplifying the condition for the present study, the two-dimensional matching condition was used in our model by the following equation.

$${\xi }_{riv}={\zeta }_{bay},u_{bay}-u_{riv}=0\mathrm{and}{v}_{bay}-v_{riv}=0$$

(13)

Equation (13) can ensure that the sea surface elevation that influences the river and the continuity of volume flux of the flow. This matching equation can exchange the information of the bay and river model in such a way that the model simulation’s updated value from the river model is again used as the bay model input. The following equation prepares the input for the bay model after the next iteration.

$$u_{bay}=u_{riv}+\frac{Q}{B_{riv}(h_{riv}+{\xi }_{riv})}$$

(14)

This is an appropriate process through the matching to find the appropriate result from the bay and river couple model. In our model simulation, the bay model inner nested grid size and the river model grid size were the same so that the matching could exchange the flow velocity and volume easily with less time and cost. A contraction mapping is needed when the grid size does not match. In this study, the numerical matching computation started from the grid number (40, 7) for the river model and (1, 7) for the bay model.

4. Numerical Procedure

4.1. Nested Model Set Up

In the present study, the study domain of the extended bay model was taken from 15° E to 23° E and 85° E to 95° E (Figure 4), which is referred to as the coarse grid model (CGM) [19]. In this CGM, there were 60 × 61 computation grid points with 15.08 km grid spacing along the north–south direction and 17.52 km along the east–west direction. According to this grid spacing, all major offshore islands along the coast of Bangladesh as well as coastal curvilinearity could not be incorporated properly. Therefore, fine mesh grid is important but the fine resolution for the whole analysis area may result in more computing time with numerical instability.

To reduce computational cost as well as incorporate the coastal complexities of the desired region properly, the nested technique was applied here. In the FGM, the mesh sizes were taken as 2.15 km along the north–south and 3.29 km along the east–west directions and the total number of computational grid points in the model was 92 × 95. Such a grid resolution could incorporate only the major offshore islands. The results obtained from the CGM were interpolated in space with cubic spline interpolation to be able to supply the grid points of the FGM. The parameters $\xi$, $u$ and $v$ were prescribed from those obtained in FGM. To properly incorporate complexities like flat and shallow coastal plains, the nature of the surge conduciveness of the region of interest and the lower Meghna estuarine area, a high-resolution very fine grid model was needed, which was why a fine mesh grid (FGM) was then improved. To make the results more accurate, a very fine high-resolution grid was made, known as VFGM. In this very fine grid model (VFGM), the computational grid points were 190 × 145 (

Table 1. Grid spacing of different nesting models with a number of computational points.

Model	Domain	Grid Spacing along x Axis	Grid Spacing along y Axis	Number of Computational Points
CGM	15° N to 23° N and 85° E to 95° E	15.08 km	17.52 km	60 × 61
FGM	21°15′ N to 23° N and 89° E to 92° E	2.15 km	3.29 km	92 × 95
VFGM	21.77° N to 23° N and 90.40° E to 92° E	720.73 m	1142.39 m	190 × 145
VFGMR	23° N to 23.25° N and 90.40° E to 90.68° E	720.73 m	1142.39 m	40 × 27

) with a grid spacing along the north–south direction of 0.72 km, it is of interest to note here that the resolution of the CGM and FGM were not the same. Such a grid resolution may be able to properly resolve coastal complexities. All of the models (CGM, FGM and VFGM) had the same dynamical Equations (1)–(3) but with different boundary conditions. In the river model, the study domain was extended between 23° N to 23.25° N and 90.40° E to 90.68° E (Figure 2). For the extended region, we considered 40 × 27 computational grid points with a grid spacing of 0.72 km along the north–south direction and 1.14 km in the east–west direction, which was referred to as a very fine grid model for river (VFGMR). The grids of the VFGMR were taken in such a way that along the west–east direction for the y-axis direction, the north–south direction was an x-direction river flow and the VFGM parent grid coincided with the grids of the VFGMR (matching points at the river and bay). The VFGMR depends on VFGM and the parameters $\xi$, $u$ and $v$ of the VFGM can pass through the boundary 23° N to 23.25° N and 90.40° E to 90.68° E.

4.2. Numerical Procedure

Both the bay and river models, the governing Equations (1)–(3) and (5)–(7) as well as boundary conditions Equations (8)–(10) and (11)–(15) except (13) were discretized by a forward in time, central in space finite difference scheme by considering discrete points in the $xy$ plane. If $(i,j)$ represent the position of a grid point in the bay model with $1\le i\le M$ and , then $x_i=(i-1)\Delta x,i=1,2,3,\dots m$ (even), $y_j=(j-1)\Delta y,j=1,2,3,\dots n$ (odd) and the sequence of time instant is defined by $t_k=k\Delta t,k=1,2,3,\dots$. The obtained equations were solved by a conditionally stable semi-implicit manner using a staggered C-grid (Figure 5) where there were three distinct types of computational points, namely $\xi$, $u$ and $v$ for the parent (bay) model and ${\xi }_{riv}$, $u_{riv}$ and $v_{riv}$ for the river model. With $i$ even and $j$ odd, the point is a $\xi$ point at which $\xi$ is computed. If $i$ and $j$ are both odd, the point is a $u$ point at which $u$ is computed and if $i$ and $j$ are both even, the point is $v$, at which $v$ is computed. In a similar fashion to the river model, ${\xi }_{riv}$, $u_{riv}$ and $v_{riv}$ were computed. The $M$ number of computation grid was chosen to be even so that the eastern open sea boundary consisted of both $\xi$-points and $v$-points. The $N$ number of the computation grid was chosen to be odd to ensure that there were only $\xi$-points and $u$-points for all of the grid points of the river mouth, where the continuity of volume flux was ensured. The initial value of $\xi$, $u$ and $v$ were taken as zero to represent a cold start. A 60 second-time step process was used for both the model simulation in order to meet the Courant-Friedrichs-Levy (CFL) stability criteria. A stable tidal regime was then generated with the $M_2$ constituent along the southern open boundary of the CGM. The specified boundary data of FGM were then integrated forward in time. The computed boundary data was exchanged from each nesting model. For the CGM, FGM and VFGM, the computation time step was checked by the CFL criteria for stability. Bay and river matching is important for the river model to work properly. The data obtained from the VFGM were assigned in space to obtain a full set of boundary data by the matching condition Equation (14) and the river bathymetries were interpolated in every grid point from the GECBO data. In the river model, for the solution process, the extrapolated value of the elevation from the bay model at the bay-river junction point provided a forcing for the motion in the river. The use of Equations (5)–(7) demonstrated the updated elevations and currents at the $\xi$ and $u$-points in the river, except for the current at the head of the river where the $u$-velocity was given by one of the boundary conditions specified by Equations (11) and (12).

The updated elevation at the junction point from the river model was chosen to update the elevations in the bay model by using Equations (13) and (14) for the next iteration. The first point of $\xi$ and $u$ in the river model, the elevation and current were computed from the bay model by using the matching phenomenon. This matching procedure guided the updated elevations in the river and, additionally, yielded boundary values of the mouth of the river to be used in the next update of the bay model variables. To find the impact of the river channel, we opened the vertical wall at the northern grid of the bay model where the matching conditions were established. For the next simulation, we set the vertical wall at the mouth of the river to consider no river influence on the bay.

4.3. Numerical Conditions

To perform the model simulation, we set some numerical conditions. Numerical conditions of the model simulation were considered such as the program running 4290 times, the program storing data after every 10 min, 289 was the number of data the program stores and the program started to store data after it was run 2040 times. The calculation period was set from 27/04/1991, 8:00 a.m. local time to 30/04/1991, 8:00 a.m. local time. The time step of our computation was 60 s after checking the stability criteria. The principal lunar semi-diurnal tide M₂ was used in the study with a period of 12.42 h, speed of 28.98 deg/h and the Doodson number of this constituent was 255.55. The tide-surge interaction process in this study was similar to the study by Roy [21]. In our model simulation, the M₂ tide was calculated from Equation (10) and the tidal oscillation was taken as 12.42 h. Commonly, the M₂ tide was dominant our study region, so we selected this as the major constituent in this study. To generate the pure oscillation of the M₂ tide, it was important to calculate the amplitude and phase of the constituent. Therefore, the process of Roy [21] was used to find the precious specifications of the values.

5. Model Validation and Outcomes

5.1. Model Results

Cyclone April 1991 and a recent cyclone MORA were considered to validate our model simulation. The results were computed for 72 h of the storm from 0800 local time (LT) of 27 April to 0800 LT of 30 April for the storm April 1991 and 27 May 2017 UTC 0600 to 30 May 2017 0600 UTC and the results are presented for the last 48 h at some of the stations shown in Figure 1. We ran the model to simulate the surge height at different tidal stations. Figure 6 shows the surge height at some of the tide stations. To understand the investigated surge height clearly in the map locations, we multiplied the simulated surge height with an arbitrary number 100 for the altitude.

The model results that came out in this case fairly agrees with the observed tidal data collected from the Bangladesh Inland Water Transport Authority (BIWTA). In our investigation, we ran the bay and river model for the different cases; the cases were with and without the inclusion of the river. For those cases, we ran the model to consider the tide and surge. Figure 7 shows the water level elevations due to surge at the representative locations along the coast of Bangladesh when there was no tidal influence. It can be seen from Figure 7 that the maximum surge height came out through our calculation (3–6 m). We found that a strong recession could occur at Hiron Point. At Hiron Point, a strong recession occurred after 2:00 a.m. (LT) on 30 April. The reason behind why this recession took place was due to the backwash of water towards the sea. The figure indicates that a recession is delayed if the cyclone went from western to eastern coastal locations. Our calculated maximum water level due to the surge was only 5.98 m at the Cox’s Bazar tide station. This simulation result did not incorporate the influence of the lower Meghna river and the tidal influence. The right part of the same figure shows the simulation result at the same location when there was no river in the model simulation. Due to the river influence, we found that the water level elevation was reduced at some locations. Figure 8 indicates the results of tide, surge and their interaction at the same location presented in Figure 7, however, this figure shows the same condition simulation results but the nonlinear interaction of the M₂ tide was introduced here. From this simulation, the water level was found to be (3–6.8 m).

In this simulation, we used the $M_2$ tide at the southern open sea boundary and nonlinear interaction with the surge height. The results obtained from [22,34,35] fairly agreed well with the computed data. Due to the impact of the tide, the height of the water was very high due to the storm crossing the coast. Finally, we found that the river had some influence on the surge simulation. When we considered the river in our model simulation, we found that the surge height fluctuated at different tidal stations. Near the mouth of the river channel, the surge height decreased by up to 10–20% than at the other tidal stations. Figure 9 shows the investigation results near the mouth of the river (Sona Char Island (Figure 1)). The direct impact of the storm elevation in the river estuarine or nearby was noticed. For this reason, we found that the wave height decreased around the estuarine area. The reason behind the reduction was that the piled-up water penetrated inside the river and flooded the innermost area of the land. However, some areas of the Chandpur district 40 km away from the estuary were also inundated due to storm surge. Figure 9 explains the effects of the river on the water level elevations due to cyclone 1991. The left-hand side of Figure 9 shows the water level height outside of the bay–river junction area and the right side shows the water level height at the junction area. Figure 9 shows that there was an acceptable impact on the storm by the river. The figure explains simulation results of the different conditions mentioned as the tide, without tide, river and without river impact. The Bangladesh Inland Transport Authority (BIWTA) reported that the maximum water level at Cox’s Bazar was six meters and the water level along the coastal region of Bangladesh was computed as 3–6.8 m. Our model computation reported that the water level height at the Cox’s Bazar tide station with both the river and without river influence was 10–30 cm. However, at the junction point and near the estuary, the water level height fluctuated between 1–1.2 m. To test the accuracy of our model simulation results, Figure 10 shows the combined results of our simulation and the observed data.

A direct comparison of our model simulation with the observed data was limited. The reason behind this may be due to the rough weather and the manual process of data collection. The red circle shows the observed data collected from the Bangladesh Inland Water Transport Authority (BIWTA), which were very scarce. The figure legend ‘WRDCH’ represents the surge height considering the river impact at the Cox’s Bazar tide station and ‘WORDCH’ stands for without river discharge at the Chittagong tide station. Similarly, ‘OBCH’ represents the observed data at the Cox’s Bazar tide station. The results ‘WRDJ’, ‘WORDJ’, ‘WRDWOTJ’ and ‘WORDWOTJ’ represent the simulated surge height at the junction point (Sona Char (Figure 1)) with the condition of the inclusion of river (WRDJ), without river (WORDJ), inclusion of the river and tide (WRDWOTJ) and without inclusion of the river and tide (WORDWOTJ).

5.2. Results Discussion

In this study, the impact of river systems has been observed by the two-dimension bay-river couple model to determine the storm surge height. From the above-mentioned figures, it is observed that the river has a visible effect on the water level height due to a storm. Therefore, the river channel and the coastal structure [36,37] of Bangladesh has an inpact on storm surge height.

Table 2. Overall maximum wave peak statistic from the past study and the present study.

Coastal Location	Overall Max. Water Level (m) by [25]	Simulated Overall Max. Water Level (m) by [22]	Simulated Overall Max. Water Level (m) by FDM [38]	Simulated Overall Max. Water Level (m) (without River)	Simulated Overall Max. Water Level (m) (with River)	Observed Overall Max. Water Level (m)
Cox’s Bazar	--	--	6.14	5.98	5.89	6.00
Moheshkhali	--	--	4.12	4.59	4.57	--
Banshkhali	--	--	3.58	--	--	--
Chittagong	6.25	5.45	4.50	4.60	4.61	5.4
Sitakunda	5.78	--	4.48	5.15	5.11	--
Sandwip	5.63	5.33	4.38	5.21	5.20	--
Mirsharai	--	--	5.66	5.05	5.03	--
Companiganj	7.28	--	6.15	5.90	5.89	6.1
Chital Khali	--	--	4.50	--	--	--
Char Jabbar	6.35	5.18	5.69	5.51	5.47	--
Char Changa	5.81	4.31	4.12	4.60	4.58	--
Char Madras	5.81	--	4.32	4.99	4.95	--
Rangabali	4.50	--	4.07	3.56	3.56	--
Kuakata	3.86	--	3.96	3.60	3.58	--
Patharghata	--	--	4.36	3.55	3.57	--
Tiger Point	4.57	--	4.21	4.30	4.30	--
Hiron Point	4.01	0.70	3.80	3.48	3.45	3.5
Monpura Island	--	--	--	5.20	4.88	--
Sona Char	--	--	--	4.51	3.61	--

provides the overall peak of the water level height at some location near the coastal area. The results of Table 2 have a similarity to our study with the results of other research.

To understand more clearly the influence of river on the surge simulation, the water level elevation distribution results are represented by the Figure 11.

It is revealed by our findings that, the surge height influenced by the river system due to the inward inflow through the river channel. Though the freshwater has some influence, the surge height was dominated by the river discharge effect. Figure 11c shows the maximum reduces of surge height (m) near the junction area. This effect could be found due to inland penetration of water flow. In our study, we have used the average discharge of fresh water but there may have the seasonal impact. Therefore, it is important to find the seasonal river discharge impact on surge height near the junction area. However, this study investigates surge height considering the yearly average river discharge. The model can be improved by incorporating the precipitation model in future.

The Figure 12 shows the effect of river discharge on water level elevation due to a storm. For the numerical computation, we set the river discharge at zero and a yearly average value of 5100.00 m³/s. The left figure of Figure 12 shows the time series result of river discharge impact on the water level elevation due to a storm near the junction are Monpura island (see Figure 1). The right figure of Figure 12 shows the discharge impact at some reported location near the coast. Obviously, the surge height fluctuation near the junction area depends upon the characteristics of cyclone behavior like strength and landfall direction of a cyclone. Maybe the reason behind this, we have found a negligible impact of the river in another storm (MORA) surge simulation. We have mostly focused on the cyclone 1991, because of the huge study has been done on this storm by the researcher. We have run our model for the recent cyclone MORA and found the result reported in Figure 13. This figure shows the maximum water level height at Cox’s-Bazar tide station and the Junction area (Sona char). The result represents the water level height due to the tide and surge interaction under the condition of considering the river and without a river.

The results obtained by the different investigations and our computed water levels due to a surge at some stations with other coastal locations compared well. The height computed for both models were compared with and without the inclusion of the rivers for the April 1991 cyclone and MORA. The maximum surge height computed from the bay model at junction point (Monpura Island) was 5.2 m while at the same location, the value computed from the bay–river model was 4.88 m (Table 2). Thus, there was an average of a 10% to 20% reduction in the overall maximum surge height computed the from bay–river model due to the penetration of piled up water inside the Meghna River. In addition, the penetration depended upon the characteristics of the cyclone. Therefore, we found that there was a negligible impact of the river on the water level elevation due to cyclone MORA. Due to the river channel near the head of the river (Chandpur district), the peak surge values computed from the bay–river coupled model were 2.8 m approximately. It can be noted that the effect of the presence of the Meghna River was confined. Thus, the maximum water level elevations computed from the bay–river coupled models along the coast of Bangladesh were different from the bay model simulation where the Meghna River joins with the bay.

6. Conclusions

In this study, a two-dimensional bay–river coupled model was developed to test the effect of the Ganga–Bhramutra–Meghna river system on the water levels associated with a storm along the coast of Bangladesh. The study showed that the effect of a river on the water levels associated with a storm was dependent on the strength of the cyclone. The inclusion of the river in our model showed that the maximum water level due to a storm was reduced by up to 20% at the near location of the river mouth. It was confirmed by this study that the river had some influence on the water level elevation due to a storm. The presence of the river reduced the water level elevation due to a storm but the river discharge influence on the water level elevation increased. Thus, the response of a river on the surge development is appreciable and cannot be ignored for storm surge prediction purposes. However, the model has some limitations but the model can be used for operational predictions of storm surges along the coast of Bangladesh.