The vertically integrated SWEs for the dynamical process in the sea can be given by (Paul and Ismail, 2013; Paul et al., 2014, 2018a, b)

where $x$ and $y$ represent the coordinate axes directed towards the south and the east, respectively, with the reference point $({23^{\circ }}{\rm{N}}, \, \, {85^{\circ }}{\rm{E}})$ on the mean sea level at the north-west corner of the study domain; $u$ and $v$ symbolize the Reynold’s averaged velocity components in x- and y-directions, respectively; $h{{(x, y)}} $ denotes the undisturbed water depth and $\zeta {{(x, y, t)}}$ stands for denoting the elevation of the free surface at any point over time $t$ about the undisturbed surface level; $f = 2\varOmega \sin \phi $ is the Coriolis parameter, where $\varOmega $ represents the angular speed of the earth rotation and $\phi $ is the latitude of the place of interest; $g$ is the acceleration due to local gravity; $\rho $ denotes the density of sea water; ${C_f}$ the friction coefficient; ${T_x}$ and ${T_y}$ stand for the components of wind stress acting on the water surface along $x$ and $y$ directions, respectively.

In this study, the surface stresses are parameterized using a conventional quadratic law (Jelesnianski, 1965):

where ${C_D}$ is the surface drag coefficient, which varies with the wind speed and the stability of the atmosphere (Das et al., 1974). In our study, we have used ${C_D}{\rm{ = 0}}{\rm{.002\;8}}$, which is found to be used by most modelers (Das, 1972); ${\rho _a}$ is the density of the air; and ${u_a}, \, \, {v_a}$ are the velocity components of the surface wind in x- and y-directions, respectively. Though there are various empirical formulae for generating wind stress depending upon the meteorological information, we have adopted the formula due to Jelesnianski (1965), being given by

where ${V_0}$ is the maximum sustained wind speed and $R$ is the corresponding radial distance from the eye of the storm; ${r_a}$ represents the radial distance at which the wind field is to be generated. It is to be noted here that the meteorological information that are required in generating the wind field led by the formula mentioned above are available in BMD. The generation of wind field led by the above formula over the region of interest can be found to be supported by various studies (Roy, 1995, 1999a; Paul and Ismail, 2012a; Paul et al., 2014, 2016, 2018a).

There are two types of boundaries, namely closed and open boundaries. The coastal belt (northern boundary) and island boundaries are treated as the closed boundaries, where the normal component of the Reynold’s depth-averaged velocity is assumed to be zero. The other three boundaries are treated as open boundaries, where the effects of different factors are required to be included depending upon the availability of the data and specific modelling application (Mills, 1985). Due to unavailability of sufficient information, in most of the studies conducted over the BOB region, a radiation type boundary condition is found to use for the open boundaries (Roy, 1995), which allows internally generated disturbances to propagate outwards normally across the open boundaries without disturbing them (Sinha et al., 1986). In our study, the following open boundary conditions are used (Paul and Ismail, 2012b, 2013):

the west boundary:

the east boundary:

the south boundary:

where $a$ and $\varphi $ denote the prescribed amplitude and phase, respectively, of the tidal constituent of interest, and $T$ denotes its period.

Along the northern boundary of the innermost child scheme, as we will see later, Meghna River is considered between longitudes ${90.46^{\circ}}{\rm{E}}\, \, {\rm{and}}\, \, \, {90.61^{\circ}}{\rm{E}}$ along latitude ${23^{\circ}}{\rm{N}}$ and following Roy (1995), the river fresh water discharge is taken into account by the formula:

where $Q$ is the fresh water discharge through the river, and $B$ is the breath of the river, which is approximately 16 km in width at the mentioned position.

The study domain (Fig. 2) of interest is extended from latitudes ${15^\circ}\, \, {\rm{to}}\, \, \, {23^\circ}{\rm{N}}$ and longitudes $\, {85^{\circ}}$ to $\, {95^{\circ}}{\rm{E}}$. The reason behind the consideration of such a big domain is that it is required for keeping the storm at least for three days on the domain, which is prerequisite to make people aware of the situation along the coastal region (Rahman et al., 2017). To incorporate the complex land-sea interface and bottom topographic details, a high resolution of grid is required. But the use of such an improved grid resolution over such a big domain increases the computational cost tremendously. In order to solve this fact, following Paul et al. (2016), one-way nested grid technique is used in this study. The coarse mesh scheme (CMS) is the parent one extended over the whole domain of interest with relatively lower resolution (Fig. 3). It is to be noted here that in the deep sea, the elevation of sea surface is rarely affected by the factors mentioned above (Paul et al., 2016). But they are highly effective on surge towards the coast (Paul and Ismail, 2013; Paul et al., 2016, 2017). Therefore, a fine mesh scheme (FMS) with a relatively high resolution (Fig. 3) is nested into the CMS covering the whole coastal region of the country. The FMS is designed in such a way that it can be able to incorporate the coastal complexities and the effect of islands properly except that of the Meghna deltaic region. Here the mouth of the Meghna River falls into the sea and the sedimentation and fresh water discharge through the river are changing the topography of the sea bed continuously. So in this area, an extra care is needed to be taken into account, which is made sure with the nesting of another very fine mesh scheme (VFMS) into the FMS (Fig. 3). A detailed description of the three schemes is presented in Table 2, whereas Fig. 2 may help to make understand the study domain with nested schemes. The coupling of the schemes is made by the process used in Paul et al. (2016).

The governing equations specified by Eqs (3)–(5) along with the boundary conditions prescribed by Eqs (8)–(10) are discretized with the aid of Eq. (2) providing the discrete points in the computational xy-plane defined by ${x_i} = (i - 1)\Delta x, \, i = 1, \, \,$ $ 2, \, \, 3, \,..., M$ (even) and ${y_j} = (j - 1)\Delta y, \, $ $j = 1, \, \, 2, \, \, 3, \, \, ..., N$ (odd) using the staggered C-grid in which there are three distinct types of computational points. With $i$ even and $\,j$ odd, the point is a $\zeta $-point at which $\zeta $ is computed. If $i$ and $\,j$ both are odd, the point is a $u$-point at which $u$ is computed and if $i$ and $\,j$ both are even, the point is a $v$-point at which $v$ is computed. $M$ and $N$ are chosen to be even and odd, respectively, to ensure $\zeta $ and $v$ points along the southern open sea boundary and $\zeta $ and $u$ points along the eastern and western open sea boundaries. It is to be noted here that to discretize model equations by finite difference approximation method, the domain should be rectangularized. Stair step method can be of help in this regard. A fictitious coastal boundary and its stair step representation on an Arakawa C-grid are shown in Fig. 4 for better understanding. It is pertinent to point out here that in approximating the complex land-sea interface on the Arakawa C-grid, following the rule of stair step representation, we used our developed MATLAB routine. Our approximated coastal and island boundaries after imposing stair step method on the different schemes are shown in Fig. 3, which show incorporation of the coastal complexities differently in different schemes. However, after the discretizations, the continuity equation and the $x$ and $y$ components of the momentum equation can be written in the following forms:

For every grid point $({x_i}, \, {y_j})$, where $\, i = \, 2, \, 4, \, 6, \, ..., M - 2$, $\, j = 1, \, 3, \, 5,..., N - 2$, and Eq. (3) can be written as

where $l = i - 3, \, \, i - 2, \, \, i - 1, \, \, i, \, \, i + 1, \, \, i + 2, \, \, i + 3;\, \, m = j - 3,$$ \, \, j - 2, \, \, j - 1, \, \, j, \, \, j + 1, \, \, j + 2, \, \, j + 3.$

Similarly, for every grid point $({x_i}, \, {y_j})$, where $ i = \, 3,\,\, 5, \,$$ 7, ..., M - 1$, $\, j = \, 3, \, 5, \, 7, ..., N - 2$, and the $x$-component of momentum equation specified by Eq. (4) can be set to the following form:

where $l = i - 3, \, \, i - 2, \, \, i - 1, \, \, \, i, \, \, i + 1, \, \, i + 2, \, \, i + 3;\,{\rm and} \, m = j - $$3, \, j - 2, \, \, j - 1, \, \, j, \, \, j + 1, \, \, j + 2, \, \, j + 3.$

Finally, for the grid point $({x_i}, \, {y_j})$, where $\, i = \, 2, \, 4, \, 6, ...,$$ M - 2$, and $\, j = \, 2, \, 4, \, 6,..., N - 1$, $y$-component of momentum equation specified by Eq. (5) can be written as

where $l = i - 3, \, \, i - 2, \, \, i - 1, \, \, \, i, \, \, i + 1, \, \, i + 2, \, \, i + 3;\,{\rm and} \, m = j - $$3, \,\, j - 2, \, \, j - 1, \, \, j, \, \, j + 1, \, \, j + 2, \, \, j + 3.$

In the above equations, the prefix $k$ means that the values of the variables are at the current time level and $k + 1$ means that the values of the variables are at the advanced time level.

In order to run the model, some data are to be allocated as input data including oceanographical, meteorological, hydrological, and geographical. To specify the water depth data accurately to the grids discussed in Section 4.1, the area covering our domain is cropped from Fig. 3, which was used in the study of Johns et al. (1985). Then by a MATLAB routine, a matrix with water depth information specified by the given contours is produced. Depth information at the remaining elements of the matrix is generated using Shepard interpolation. A contour map of our interpolated bathymetry is shown in Fig. 5. The bathymetric data at each of the longitudinal and latitudinal positions of the computational grid points of the CMS, FMS and VFMS are then compiled from the generated matrix using again the Shepard interpolation. The meteorological inputs namely, storm track, maximum sustained wind speed and maximum wind radius are collected from the BMD and in this regard, Table 3 and Fig. 6 are included for better understanding. The rotation speed of the earth and density of the sea water are taken as $\varOmega = 1\;040\, \, {\rm{m}}/{{\rm{h}}}$, and $\rho = 1\;025\, \, {\rm{kg}}/ {{\rm{m}^3}}$(Paul et al., 2014), respectively. The value of the friction coefficient ${C_f}$ is taken as constant for computational simplicity, which is 0.002 6. The tidal period is taken as $T = 12.4\, \, {\rm{h}}$ (Roy, 1995). Initial values of the momentum and scalar variables are taken as zero to represent an initial state of static equilibrium, which is supported by the studies due to Loy et al. (2014) and Paul et al. (2016, 2018a, b). It is to be mentioned here that the mentioned initial state is the best assumption; otherwise, the real values of $\zeta $, $u$ and $v$ have to be made available at all grid points in the region of analysis (Loy et al., 2014). The time step is taken as $\Delta t = 60\, \, {\rm{s}}$ that ensured CFL stability criterion.

The obtained discretized equations given by Eqs (12)–(14) are solved with the RK(4,4) method.

First, Eq. (12) with the initial conditions and the defined values of the parameters involved (Appendix, Eq. (A1)) is solved by the RK(4,4) method for the values of $\zeta $. A detailed description of the solution procedure can be found in Paul et al. (2014). After getting the values of $\zeta $ at the points specified by the Arakawa C-grid, its values are calculated at the specified grids of the boundaries through Eqs (A3)–(A5) (Appendix) and then averaging procedure is taken into account to obtain $\zeta $ at each of the grid points representing water, coastal and island boundaries. Path of the storm is then generated with the obtained data from the BMD (Table 3) and the wind field is generated subsequently with the values of the required parameters from the BMD with the help of Eqs (6) and (7). Equation (13) is then solved with the help of the RK(4,4) method. Finally, Eq. (14) is solved for the y-component of momentum equation. The procedure mentioned above is made for all the schemes. The only difference is in the boundaries as is mentioned earlier. The Meghna River fresh water discharge is incorporated in the northern open boundary of the VFMS using Eq. (A6) at the points $(1, \, j)$ where $j = 7, 9, 11, 13, 15,$ 17, 19. This procedure is repeated over time supplying the updated values as initial ones for calculating water levels due to tide, surge, and their interaction. It is pertinent to point out here that for generating a stable tidal regime over the model domain, the model was run setting the initial conditions to zero, known as "cold start" in the absence of atmospheric pressure gradient force and wind stress. First, the values of amplitude (a) and phase ($\varphi $) related to the ${\rm M_2}$ tidal constituent are prescribed in Eq. (10) along the southern open boundary of the CMS from the study of McCammon and Wunsch (1977). The period of the tidal oscillation T is taken as 12.4 h, as the mean period is always found to be approximately 12.4 h, which is of ${\rm M_2}$ tide. Then from the cold start in the absence of wind field, a stable tidal regime is achieved after four cycles of integration. But the problem to generate a pure oscillatory response in the BOB region corresponding to the tidal constituent with period T lies with adjustment of the exact values of a and $\varphi $. Thus some techniques are needed to be adopted for precise specification of the values of the constants and in this regards, we use the technique of Roy (1995) to adjust the values of a and $\varphi $. Generation of a stable tidal condition by the process along the area of interest is supported by several studies (Johns et al., 1985; Roy, 1995, 1999b; Paul and Ismail, 2012b, 2013; Paul et al., 2016). For the pure surge (in the absence of astronomical tide), the model was also run from the cold start. It pertinent to pinpoint out here that the zero initial values will not affect the calculated results as their effect on the results over several hours will almost disappear (Zhang et al., 2007). For estimating water levels due to the interaction of tide and surge, the generated pure tidal regime provided the initial condition at the model time $t = 0$.

The numerical experimentations were carried out with the cyclone April 1991. The reasons behind the choice of the cyclone are that it was a severe one and it passed nearby the Meghna River mouth, the region of our interest. Also, the investigated information about the cyclone is relatively more. The cyclone was so ferocious that about 150 000 people and 70 000 cattle were died having an economical loss of about Tk 80 billion (Khalil, 1992). For discussions of our model simulated results, the time history of the storm is utmost urgent. A detailed description of the time history of the storm can be found in Paul et al. (2016, 2017, 2018a), therefore, a gist of it is presented here for convenience. The cyclone was first detected as a depression on 23 April 1991 from the satellite pictures taken at Space Research and Remote Sensing Organization (SPARRSO) from National Oceanic and Atmospheric Administration (NOAA)-II and Geostationary Meteorological Satellite-4 (Khalil, 1993). It was gradually propagated towards the coast gaining energy. At 0300 universal time coordinate (UTC) of 25 April, it was identified as a deep depression near $({10^{\circ}}{\rm{N}}, \, \, {89^{\circ}}{\rm{E}})$. The depression was gradually intensified and at 0300 UTC of 28 April, it was detected at $({15.5^{\circ}}{\rm{N}}, \, \, {87.5^{\circ}}{\rm{E}})$ having central speed of 30 ${\rm{m}}/ {{\rm{s}}}$. Suddenly it changed its direction and started moving easterly towards the estuary of Meghna River. Finally, it hit the coast of Bangladesh near Chittagong along Noakhali coast at 2000 UTC of 29 April. The maximum sustained wind speed 65 m/s of the storm was examined at Sandwip.