INTRODUCTION
A storm surge is an abnormal rise of the water level above the predicted water
level of the astronomical water caused by a combination of wind and low atmospheric
pressure. Storm surges generated by tropical cyclones are most common in the
Bay of Bengal and the Gulf of Mexico. Tropical storms along with tidal surges
are the most devastating among all natural calamities befalling Bangladesh.
The costal belt of Bangladesh is frequently lashed by tropical storms and associated
surges. Specifically, the Meghna estuarine area is the most affected region.
Extreme bending of the coast line, shallowness of water, huge discharge through
the Meghna and other rivers, low lying big and small offshore islands and coastal
regions, complex tidal phenomena etc. make the coastal belt of Bangladesh favorable
for high surge. Thus, a proper understanding of the facts affecting surge development
and its accurate prediction in the coastal region is highly desirable. Many
analyses on prediction of surge due to tropical storms have been made all over
the world. Most of the study have been done for the Atlantic, Pacific, North
Sea and Bay of Bengal. For Atlantic region, the pioneering works are due to
Jelesnianski (1965, 1966, 1967)
and Thacker (1977, 1979). Extensive
studies have been done for East Chaina Sea, North Sea and NorthWest European
Continental Shelf. Particular mention is made here of the study of Yu
and Zhang (2002), Flather (1981), Flather
and Davies (1976), Heaps (1973), Heaps
and Jones (1981), Nihoul (1977, 1982).
Also, several numerical models Haque et al. (2005),
Das et al. (1974), Jhons
and Ali (1980), Jhons et al. (1981, 1982,
1985), Dube et al. (1985,
1986, 2004), Sinha
et al. (1986) and Roy (1985), Rao
et al. (1997), Murty and Flather (1994),
have been developed for the prediction of storm surges in the Bay of Bengal.
Though considerable number of studies have been done for the coast of Bangladesh,
none of them, except Talukder et al. (1992) and
Roy (1995, 1999), considered
the existence of Islands along the head Bay region. But it is known that there
is influence of offshore islands over the surge intensity along the coastal
belt (Talukder et al., 1992; Roy,
1999). Moreover, since the islands are thickly populated, it is necessary
to estimate the water levels due to tropical storms along these Islands also.
Therefore, the inclusion of the islands in a model is essential in order to
incorporate the real solution of the head bay region. To estimate maximum flood
levels due to both tide and surge Roy (1995), nested
a fine resolution scheme for Meghna estuary into a coarse mesh scheme extending
up to 15° N latitude. This model is similar to that of Jhons
and Ali (1980), for the east coast of India and the coast of Bangladesh.
But until now no attempt has been made to develop an efficient model suitable
for operational forecasting purpose for the coast of Bangladesh. One of the
major problems for numerical modeling of storm surge lies with the fact that
fine resolution is required in the finite difference scheme near the coast in
order to properly incorporate the bending of the coast line and offshore islands
whereas this is unnecessary away from the coast.
The present study attempts to develop a surge forecasting model based on the Cartesian coordinate system that incorporates the whole coastal belt and offshore islands very accurately. For that purpose, a Fine Mesh Scheme (FMS) for the coastal belt of Bangladesh is nested into the Coarse Mesh Scheme (CMS) covering up to 15° N latitude. As the region between Barisal and Chittagong is full of so many small and big islands and also the bending of the coastline is very high, considering those facts, a Very Fine Mesh Scheme (VFMS) for this region is nested into FMS. Vertically integrated shallow water equations in Cartesian coordinates are solved with the above mentioned doubly nested schemes. The specialty of the nesting is that, the CMS is completely independent, whereas along the open boundaries of the FMS the elevations are prescribed from those obtained from the CMS at each time step of the solution process. Similarly, along the open boundaries of the VFMS the elevations are prescribed from those obtained from the FMS. The model is verified using waterlevel data from a number of previous storms hitting the coast of Bangladesh during the period 19701991. For analysis and verification of the results, 10 locations along the coastal belt are considered (Fig. 1). The locations are Hiron Point, Kuakata, Char Madras (South Bhola), Char Chenga (South Hatiya), Char Jabbar, Companigonj, Sandwip, Chittagong, Bashkhali and Cox's Bazar.
MATHEMATICAL FORMULATION
Vertically integrated shallow water equations: A system of rectangular
Cartesian coordinates is used in which the origin, O, is in the undisturbed
level of the sea surface as the xyplane and Oz is directed vertically upwards.
The displaced position of the free surface is considered as z = ζ(x, y,
t) and position of the sea floor as z = h(x,y). Then the vertically integrated
shallow water equations given by Roy (1995) are:
Where:
u, v 
= 
Velocity components of sea water in x and y directions (m
sec^{1}) 
f 
= 
Coriolis parameter 
g 
= 
Gravitational acceleration (m sec^{2}) 
T_{x}, T_{y} 
= 
Components of wind stress (N m^{2}) 
h 
= 
Ocean depth from the mean sea level (m) 
C_{f} 
= 
The friction coefficient 
ρ 
= 
Water density (kg m^{3}) 
Also, in the above equations u and v are the vertically integrated
components given by:
where,
and
are x and y components of the Reynolds averaged velocity.
Using Eq. 1 we may express the Eq. 2 and 3 in the flux form and thus, Eq. 13 may be written as:
where, .
Here, u and v in the bottom stress terms of Eq. 2 and 3
have been replaced by
and
in Eq. 6 and 7 in order to solve the equations
numerically in a semiimplicit manner.
Boundary conditions: At the closed boundaries of the mainland as well as of the islands the normal components of the mean current are taken as zero. The following radiation type of boundary conditions are taken for open sea boundaries:
Generation of wind stress: The wind stress is, in general, parameterized in terms of the wind field associated with the storm. This is frequently done by the conventional quadratic law:
where C_{D}, ρ_{a} are the drag coefficient and air density, respectively.
The circulatory wind field can be generated by various empirical formulae.
For the Bay of Bengal region most frequently used formula is due to Jelesnianski
(1965), which is given by:
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. The x and y components (u_{a}, v_{a}) of the wind field are derived from V_{a} given by the above empirical formula.
NUMERICAL METHOD
Grid generation: In order to estimate the water level due to surge in the study area, it is necessary that the study area should be considerably big so that, a storm can move over the area at least for, say 3 days before crossing the coast. This is because; the surge response along the coast becomes significant well before a storm reaching the coast. On the other hand, in order to include the major islands in the estuary the mesh size (the distance between two consecutive grid points) should be considerably smaller whereas, this is unnecessary away from the estuary. Considering the above facts, a highresolution numerical scheme (FMS) is nested into a Coarse Mesh Scheme (CMS). Further, a Very Fine Mesh Scheme (VFMS) is nested in the FMS. The CMS covers the area between 15° N to 23° N latitude and 85° E to 95° E longitudes. The mesh size along NorthSouth direction (along xaxis) is Δx = 15.08 km and EastWest direction (along yaxis) is Δy = 17.52 km. There are 60x61 grid points in the computational xyplane. The FMS covers the area between 21°15’ N to 23° N latitude and 89° E to 92° E longitudes. This scheme area includes the study area and almost all the offshore islands. The mesh size along northsouth direction (along xaxis) is Δx = 2.15 km and along eastwest direction (along yaxis) is Δy = 3.29 km. There are 92x95 grid points in the computational xyplane. The VFMS covers the area between 21.77° N to 23° N latitude and 90.40° E to 92° E longitudes. The mesh size along northsouth direction (along xaxis) is Δx = 720.73 m and along EastWest direction (along yaxis) is Δy = 1142.39 m. There are 190x145 grid points in the computational xyplane.
As the computational method uses a rectangular grid, the coastline and the island boundaries have been approximated following the grid lines using stair step representation. It may be mentioned that, it is not possible to show clearly the actual grid specification (60x61, 92x95 and 190x145) on the real boundary line, as this needs a very big space. A typical grid specification on the physical domain is shown in the Fig. 1.
Numerical scheme and model data setup: All of the schemes have the
same dynamical Eq. 57 with different boundary
conditions. An important feature of this doubly nesting is that the CMS is completely
independent. On the other hand, along the open boundaries of the FMS the parameter
æ is prescribed from those obtained in CMS in each time step of
the solution process. Similarly for the VFMS, the parameter æ is
prescribed from those obtained in FMS in each time step of the solution process.
The governing equations given by Eq. 57
and the boundary conditions given by Eq. 810
are discretised by finitedifference (forward in time and central in space)
and are solved by conditionally stable semiimplicit method using staggered
grid system. For numerical stability, the x and y components of the momentum
equations are modeled in a semiimplicit manner. For example, from Eq.
6 the term:
is discretised as:
where, the subscript k+1 indicate that
is to be evaluated in advanced time level. Moreover, the CFL criterion has been
followed in order to ensure the stability of the numerical scheme. Along the
closed boundary, the normal component of the velocity is considered as zero
and this is easily achieved through appropriate stair step representation as
mentioned earlier.
Table 1: 
History of the chosen storms 

*Maximum wind speed: 62.00 m sec^{1}. **Maximum wind
speed: 65.00 m sec^{1} 
The initial values of ζ, u and v are taken as zero. The time step is taken
as 60 sec that ensures stability of the numerical scheme. In the solution process,
the values of the friction coefficient C_{f} and the drag coefficient
C_{D} are taken as uniform throughout the physical domain, which are
0.0026 and 0.0028, respectively.
Analysis of results and verification: For the purpose of analysis, the computed results are presented at ten locations along the coast of Bangladesh for some major storms hitting the coast of Bangladesh during the period 19701991. Records show that November 1970 and April 1991 storms are the most severe of this century having maximum wind speeds 223 km h^{1} (62 m sec^{1}) and 235 km h^{1} (65 m sec^{1}), respectively. Histories about the abovementioned storms are shown in Table 1, the data are received from the Bangladesh Meteorological Department (BMD).
Analysis for the storm of 1970: The storm of 1970 is one of the severe storms of this century that hit the coast of Bangladesh and is also favorable for high surge due to both wind intensity and path of its movement. So, we give more stress on the 1970 storm data for computation of results.
The storm of November 1970 moved approximately northward and on November 12
it gradually turned towards right and crossed the coast (land fall) of Bangladesh
near Chittagong at early morning (approximately 0030 h) of November 13 (Table
1). Figure 2 depicts the computed surge levels associated
with November 1970 storm at Hiron Point, Char Madras and Chittagong. It may
be observed that, the maximum surge level is increasing with time as the storm
approaches towards the coast and finally there is recession. At Hiron Point
a strong recession is occurred after 17 h of 12th November, earlier than in
any other locations and about 6 h before landfall of the storm (Fig.
2).

Fig. 2: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with November 1970 storm 

Fig. 3: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with November 1970 storm 
The recession takes place due to backwash of water from the shore towards the
sea. In fact, Hiron Point is situated far left (West) of the storm path and
so the direction of the anticlock wise circulatory wind becomes northerly (i.e.,
towards the sea) at Hiron Point long before the storm reaches the coast and
thus driving the water towards the sea. The recession reaches up to 2.3 m at
2200 h of 12th November. It may be noticed that the recession at Kuakata, Char
Jabbar and Bashkhali began approximately at 1800, 2200 h of 12th November and
0030 hrs of 13th November respectively (Fig. 3). Figure
4 shows that the recession at Char Chenga, Companigonj, Sandwip and Cox’s
Bazar began approximately at 2100, 2200, 2300 and 2300 h, respectively. Thus,
the beginning of recession delays as we proceed towards east as is expected.
At every location, the peak surge is attaining before the land falls time of
the storm. This is expected, as the circulatory wind intensity is highest along
the coast when the storm reaches near the coast. The maximum elevation varies
between 2.46 m (at Hiron Point) to 9 m (at Companigonj). At Chittagong the computed
water level increases up to 6.85 m before recession starts at 2330 h of 12th
November (Fig. 2).
Figure 9 shows the contours of the recommend water level associated with November 1970 storm in the Meghna estuary region (in the VFMS region of our computed domain).

Fig. 4: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with November 1970 storm 

Fig. 5: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with April 1991 storm 
We could not compare our computed results for November 1970 storm to the observed data due to non availability of observed time series data. But it is reported by the Bangladesh Meteorological Department (BMD) that there was 6 to 9 m surge near Chittagong. It is seen that the model results are in good agreement with observed data.
Analysis for the storm of 1991: Figure 57
depict the surge levels associated with April 1991 storm at ten locations from
Hiron Point to Cox’s Bazar. It is clear from our calculated results that
the maximum surge level is increasing along the locations from west to east.
Early recession of surge levels at the western coastal stations is evident from
the figures. The storm hits the coast at early morning of 30th April near Chittagong.
The intensity of the storm of 1991 is highest among the considered storms of
1970, 1981, 1985 and 1991. So, the surge is high in many locations; but surge
intensity is, in general, less than that of November 1970. This is because of
the fact that the path of 1970 storm is very favorable for very high surge along
the coastal belt of Bangladesh. Figure 10 shows the contours
of the recommend water level associated with April 1991 storm in the Meghna
estuary region (in the VFMS region of our computed domain).
According to the Banglapedia: National Encyclopedia of Bangladesh, the maximum storm surge height during April 1991 was estimated to be about 5 to 8 m. This also justify the accuracy of our model results.

Fig. 6: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with April 1991 storm 

Fig. 7: 
Computed water level (w.r.t. MSL) due to surge at different
locations associated with April 1991 storm 

Fig. 8: 
Computed maximum surge levels associated with four chosen
storms along the coastal belt between Hiron Point and Cox’s Bazar 
Analysis of computed maximum surge response along the coastal belt:
Finally, the maximum surge responses associated with four chosen storms (the
storm of 1970, 1981, 1985 and 1991) at ten representative locations from West
to East along the coastal belt of Bangladesh are shown in Figure
8.

Fig. 9: 
Contours of the recommend water level (in the VFMS region
of the computed domain) associated with November 1970 storm 

Fig. 10: 
Contours of the recommend water level (in the VFMS region
of the computed domain) associated with April 1991 storm 
For the November 1970 storm, the maximum surge increases from west to east
up to Companigonj and then gradually decrease reaching minimum at Cox’s
Bazar. The surge profile of the April 1991 storm is similar to that of the November
1970 storm. The path of the April 1991 storm is almost parallel and towards
the right of the November 1970 storm and it crossed the coast just north of
Chittagong keeping Sandwip to its left. Although the intensity of the 1991 storm
is higher than that of 1970, the surge responses of 1991 are in general less.
This is attributed to the path of the storm and the landfall position as mentioned
earlier. Thus the path of a storm as well as the landfall position is the important
deciding factor of the over all intensity of surge response through out the
coastal belt. Similar analysis can be made for December 1981 and May 1985 storms,
the landfall positions being near Hiron Point for the former and near Chittagong
for the latter.
CONCLUSION
In this study we have developed a numerical model based on the Cartesian coordinate system to compute surge levels along the coast of Bangladesh. The model is capable of incorporating bending of the coastline and island boundaries with a considerable accuracy. The nested numerical scheme ensures very fine mesh near the coast and coarse mesh away from the coast. Although, we could not incorporate the ocean depth data properly, the model may be efficient to compute surge levels in the head Bay of Bengal.