Tidal intrusion within a mega delta: An unstructured grid modelling approach

The finite volume community ocean model (FVCOM) has been applied to the Ganges-BrahmaputraMeghna (GBM) delta in the northern part of the Bay of Bengal in order to simulate tidal hydrodynamics and freshwater flow in a complex river system. The delta region is data-poor in observations of both bathymetry and water level; making it a challenge for accurate hydrodynamic models be configured for and validated in this area. This is the first 3D baroclinic model covering the whole GBM delta from deep water beyond the shelf break to 250 km inland, the limit of tidal penetration. This paper examines what controls tidal penetration from the open coast into an intricate system of river channels. A modelling approach is used to improve understanding of the hydrodynamics of the GBM delta system. Tidal penetration is controlled by a combination of bathymetry, channel geometry, bottom friction, and river flow. The simulated tides must be validated before this delta model is used further to investigate baroclinic processes, river salinity and future change in this area. The performance of FVCOM tidal model configuration is evaluated at a range of sites in order to assess its ability to capture water levels which vary over both a tidal and seasonal cycle. FVCOM is seen to capture the leading tidal constituents well at coastal tide gauge stations, with small root-mean-squared errors of 10 cm on average. Inland, the model compares favourably with twice daily observed water levels at thirteen stations where it is able to capture both tidal and annual timescales in the estuarine system. When the river discharge is particularly strong, the tidal range can be reduced as the tide and river are in direct competition. The bathymetry is found to be the most influential control on water levels within the delta, though tidal penetration can be significantly affected by the model's bottom roughness, and the inclusion of large river discharge. We discuss the generic problem of implementing a model in a data-poor region and the challenge of validating a hydrodynamic model from the open coast to narrow river channels. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
With an average freshwater discharge of around 40,000 m 3 s À1 the GBM (Ganges, Brahmaputra and Meghna) river system has the third largest coastal discharge worldwide (after the Amazon and Congo). The GBM river delta is a low-lying fertile area covering over 100,000 km 2 in India and Bangladesh, and is thus classified as a megadelta (Milliman and Meade, 1983). Time varying water levels in the delta are controlled by tide and storm surge at the open ocean in the northern Bay of Bengal, combined with a contribution from river flow through the delta distributaries.
The GBM river delta is part of the lowland Bangladesh Plain: here most land elevations are less than 10 m above sea level, and close to mean sea level (MSL) at the south of the delta. The majority of Bangladesh is comprised of extensive low-lying flood plains (Yu et al., 2010), and this area has been recognised as being highly vulnerable to sea-level rise, (see e.g. Nicholls et al. (1999)). The river system of Bangladesh covers around 7% of the country's surface, and in the delta area there is a hazard from both riverine and coastal flooding. River flooding is responsible for large amounts of river bank erosion, damage to property, and loss of fertile farm land. Coastal flooding can be damaging to crops, bringing saline water inland, while the annual monsoon river flooding brings freshwater to leach dangerous salts out of the soil (Clarke et al., 2015). Salt penetration in-land is linked to tidal intrusion, and the balance between river discharge and tidal range at the mouths.
While the delta is relatively flat and featureless above sea level, it is scoured by deep river channels. These river channels are dynamic -changing shape and position in response to coastal erosion, and sediment deposition from the high river flows during the monsoon. Information about the topography is sparse and as the delta's morphology is dynamic: regular bathymetry surveys are required. Mapping the river channels is particularly challenging due to the seasonal changes in water level and large sediment load. However, a good bathymetric data set is essential for tidal prediction and reproducing the hydrodynamics within the delta.
The high flows seen annually during the wet season have a dramatic effect on the underlying bathymetry. The rivers carry large sediment load, leading to the progradation of the delta. The soft sediments of the delta are very mobile, and the positions and hypsometry of river channels is subject to change year on year with net gains of 14.8 km 2 between 1972 and 1840, and of 4.4 km 2 between 1840 and 1984 (Allison, 1998). As well as altering the morphology of the river channel, erosion can cause the position of the channels to change as islands accrete and erode (Rahman et al., 2011). This is reinforced by Bandyopadhyay (1997) who plots the 140 year evolution of Sagar island at the mouth of the Hooghly, showing gradual erosion over time due to a reduced supply of sediment combined with deforestation. Land erosion and accretion over a 23 year period are also presented by Brammer (2014), showing a net gain of land area but in the context of rapid geomorphological change. It is therefore important that our models are based on contemporary channel maps and bathymetry. Sarker et al. (2014) give a thorough analysis of historic morphological changes in the GBM river system, and Sarker and Thorne (2009) show detail of cross-sectional area changes in the Padma/Jamuna system.
Propagation of tides in the Bay of Bengal is addressed in detail by Sindhu and Unnikrishnan (2013). The semidiurnal tides in the northern Indian Ocean have two distinct regimes e.g. Pugh (1987). The Arabian Sea is broad enough for the standing wave system to develop an amphidrome; however, this amphidrome, being situated close to the equator, cannot be described simply in terms of Kelvin waves. The entrance to the Bay of Bengal is too narrow for an amphidromic system to develop. The wave propagates to the north along the west coast of Sumatra and Thailand and also along the east coast of Sri Lanka, where the range is low and there is a tendency towards a degenerate amphidrome. In the south central Indian Ocean there is an extensive region of large semidiurnal tides over which the phases change only slowly. Dynamically this phenomenon, which corresponds to an antinode, is called an antiamphidrome.
The GBM delta is macro-tidal (mean tidal range > 4 m) with predominantly semi-diurnal tides. Large tides at the coast can penetrate far inland, with tidal variability being observed as far as 200 km inland. Observations of the spring tides on the Indian side of the GBM delta, by Chatterjee et al. (2013), suggest that tides can amplify within major estuaries and tidal creeks, though the rates of amplification are not uniform and follow a complex pattern along the different distributaries, such that the observed spring tidal range for all coastal stations is always larger than 4.5 m and can be as large as 6.7 m. Water levels in the delta are also significantly affected by the large volume of river discharge, which can raise mean water depths by as much as 4 m. Tides, particularly the overall tidal range, are important to the growth of mangrove trees which will not survive if their roots are submerged for long periods (Payo et al.;Auerbach et al., 2015;Gopal and Chauhan, 2006;Barbier, 2012). This paper presents a new 3D baroclinic model of the GBM delta, capable of capturing the hydrodynamic processes, and predicting water levels and currents. Hydrodynamic validation is necessary before examining density changes or making future projections. In order to predict saline intrusion within the river system, any baroclinic model must first be capable of capturing the 2D behaviour of the tides, as water level data is usually the most basic data source available. Here we have maximised the use of limited tidal water level observations to assess the model performance in a data-poor region. We first validate and test the reliability of the model in this region, before discussing the interactions between rivers and tides. By using a numerical model we can perform experiments and sensitivity studies, which provide valuable insight into processes in this data-poor region. The structure of the paper is as follows. Section 2 describes the model configurations and experimental design, with a focus on designing the position of coastline and sources of bathymetry data. Section 3 assesses model performance when compared with available observations. Section 4 discusses the processes controlling tidal penetration into a delta. Results are discussed in section 5 before conclusions drawn in section 6.

Model
An unstructured grid modelling approach was selected for its variable resolution capability in order to get a good representation of narrow river channels while still being able to capture the continental shelf and deep water in the Bay of Bengal with a coarser model. The Finite Volume Community Ocean Model (FVCOM, Chen et al. (2003)) has been applied in a region covering 88.9e93.6 E by 18.9e24.0 N (shown in Fig. 1). FVCOM has previously been successfully applied to small coastal applications (e.g. Yang and Wang (2013) Li et al. (2012), Zhao et al. (2010), Xue et al. (2009), Zheng et al. (2003, and also at a larger scale, simulating shelf tides in the Arctic . For this work, the extent of the model inland was chosen to cover the extent of the tidal penetration, as the tidal signal in places reaches as far as 200 km inland, and storm surges may drive the saline Bay of Bengal waters up to 160 km inland (Snead, 2010). The triangular grid resolution varies from large elements of order (10 km) in the open ocean (maximum 26 km), and small elements of order (100 m) within the narrow river channels (minimum 47 m). The water depth varies from 2 km deep off the shelf, to less than 10 m within the delta. The GBM delta is defined according to Galloway's classification (Galloway, 1975) as a tide-dominated braid delta. The delta's longitudinal elevation profile is such that, in order to capture the full extent of the tidal penetration, we must consider river channels with topographic elevations above mean sea level. This application of a coastal ocean model so far inland is believed to be unique.
Hourly water level was prescribed at the open ocean boundary to the south of the domain. Tidal water levels are taken from a POLCOMs model (Proudman Oceanographic Laboratory Coastal Ocean Model Holt and James (2001)) of the Bay of Bengal (described in Kay et al. (2015). The POLCOMS configuration simulated the 8 leading harmonic constituents(M 2 , S 2 , N 2 , K 2 , K 1 ,O 1 , P 1 ,Q 1 ).
To simulate river discharge into the model, a daily volume flux is applied at the upstream limit. The freshwater volume input as a boundary condition to the delta model comes from a separate river catchment model, INCA. The Integrated Catchment Model (INCA, Wade et al. (2002)) is a process based model which simulates hydrology flow pathways in surface and sub-surface river systems. INCA has been used to simulate freshwater forcing under present day conditions, and a series of possible climate and social futures (Whitehead et al., 2015a, b). Freshwater is introduced into the FVCOM model at a single point at the north of the domain (24.06 N, 89.03 E). The model is started from rest during a flood condition, and the freshwater is then allowed to flow through the distributaries. For this study, a single year's river flow representing the year 2000 is used as freshwater forcing. The lower panel of Fig. 6 shows the daily river volume flux entering the model.

Model grid and bathymetry
Bathymetry is the most fundamental input to any hydrodynamic model, particularly in estuaries (Parker, 1991). There is published topographic data available from The General Bathymetric Chart of the Oceans (GEBCO, 2014), and ETOPO (NGDC, 2011). However, while these data sets are suitable in deeper water offshore, they fail to capture the nearshore and inland bathymetry within the delta. River and tidal flows in the model are critically dependent on good representation of the coastline and bottom bathymetry. In order to produce the most accurate representation of the GBM delta, existing surveys plus new survey work was required. As part of the ESPA Deltas project http://www.espadelta.net/, extra river channel surveys have been performed (during 2013/14), and this up-to-date bathymetry is used here to give the best possible hydrodynamic model. Fig. 2 shows a typical river cross section taken across the Upper Meghna, compared with model bathymetry. Several separate channels can be seen, up to 30e40 m deep in places. Fig. 2 shows the locations where observations were made to improve the bathymetry maps; the left hand panel were already available. The right hand map shows those measured during the ESPA project. Several data sets were then blended to generate the final bathymetry; river cross-sections, digital elevation maps (DEM), polder information, and GEBCO data. The one-dimensional river channel network model HEC-RAS was used to analyse river cross sections, and this also informed bathymetry. Extensive manual work was undertaken at the Bangladesh University of Engineering and Technology (BUET) using ArcMap and Delft-3D Quickin (Haque et al., 2016). Data from several years were combined in order to get the best spatial coverage of the delta. Care is needed when combining bathymetric data spanning several years, due to the soft sediment and dynamic nature of the river delta morphology. However, due to the sparsity of depth measurements available a patchwork of bathymetry data was unavoidable.
Once an acceptable bathymetric data set has been produced the next step is to define a coastline; this defines the extent of the unstructured mesh above which the model will not inundate. In order to develop an unstructured mesh of the delta, the gridded bathymetry generated by Haque et al. (2016) was taken as a starting point. To produce a coastline, first a contour at À0.5 m (relative to MSL) was selected. This contour, while suitably representing the Sundarbans and large river mouths, was not suitable further inland; as the river delta slopes upward, and river channels can have their floors above MSL; for example, the tide can be felt at Pabna (23 56 0 09.2 00 N 89 12 0 17.5 00 E) where 'wet' points can be as much as 5 m above MSL. The gridded bathymetry was read into the surfacewater modelling system (SMS) software (Zundel, 2005). SMS was then used to extend the À0.5 m contour coastline further into the delta. That was done by hand, using the river channel depths as a guide to the location of the boundary.
In order to capture the channels adequately in both wet and dry seasons, the coastline must extend more widely than the central river channel. In the GBM delta, the meandering river channels have deep central sections with wide shallow intertidal mudflats, in some places reinforced with levees and polders at the banks. For example, Fig. 2 shows that the width of the flood plain can be more than twice that of the river channel. It is important for the model to capture these shallow areas, which are often mud flats at the time of tidal low water, especially in the dry season. Brammer (2014) gives more details on the physical geography of the region, focussing on coastal morphology changes, and subsidence. It is critical at this point that the gridded river channels are wide enough and contain sufficient resolution to represent the hydrodynamics of the region. Therefore the model resolution used must vary, based on channel width. For example in the mouth of the Meghna (around 15 km wide) the model resolution across the river section is around 20 points. Further inland the river channels are narrower, for example the Arial Khan is only 300e400 m across. Here, the model resolution has been increased (allowing around 10 points across the channel) to properly represent the shape of the river. This coastline was then used in SMS (http://www.aquaveo.com/ ) to produce a triangular mesh. This mesh is checked for grid quality, including the aspect ratio of model grid cells, and mesh connectivity, using SMS. Finally, the gridded bathymetry was interpolated onto the triangular grid. At this point, the mesh was further checked for regions of steep bathymetry, and the water depths smoothed. At the south of the domain (see Fig. 1) an open boundary was defined in the Bay of Bengal, where tidal and thermohaline forcings can be applied. The elements at the open boundary should be aligned perpendicular to the boundary, to improve model stability. The finalised model domain and a closeup of the river channels are shown in Fig. 1.
As this model must be capable of capturing the peak of the wet season flow, it must also stand to reason that some areas will become inert during the dry season. FVCOM is capable of 'wetting and drying' individual elements and can represent this behaviour. However, we must be confident that individual stretches of river channels remain connected throughout the simulation, to avoid unrealistic ponding as a model artefact. To ensure this connectivity, and numerical stability, the model is initialised from a flood condition during the monsoon season which ensures there is an initial flow throughout the grid.
The bottom friction in FVCOM can be set as a constant drag coefficient (C d ) or a roughness length (z 0 ). In all simulations presented here, a spatially constant z 0 value is used. Roughness length controls bottom friction through: where h is the water depth, and k is the von Karman constant, set to 0.4. s is the number of vertical levels in the model. In shallow waters, below 3 m, C d is calculated based on a minimum depth of 3 m. The minimum value of C d is capped at 0.0025. In river modelling (river channel flow), a Manning coefficient denoted as n, is often used to represent bottom friction (Chanson, 2004). The Manning formula states: where V is the cross-sectional average velocity, R h is the hydraulic radius, and S is the slope of the hydraulic grade line. The Manning coefficient is dependent on many factors, including surface roughness, sinuosity, and the hydraulic depth of the channels. FVCOM is primarily a coastal model, and not designed for river hydraulics and there was also no information available about the river channel type. Therefore we cannot directly implement a Manning's approach in the delta model. However, for comparison in 10 m water depth, a z 0 ¼ 0.015 m is equivalent to C d ¼ 0.0056 and n ¼ 0.035.

Water level data
Bangladesh is notoriously data-poor (Lewis et al., 2013) and as such validation of any model in this area is a challenge. Two data sets have been used in this work: tide tables from the UK Hydrographic Office (UKHO), and twice daily water level observations within the delta. The UK Hydrographic Office (2012) provides tidal phase and amplitude at 19 coastal tide gauges around the Bay of Bengal. The 4 leading constituents (M 2 , S 2 , K 1 , O 1 ) of the tide are generally provided, however the tide is often listed as 'V' for variable. Further inland, within the delta, daily measurements of high and low water level are available at 66 sites, recorded by the Bangladesh Inland Water Transport Authority (BIWTA). The locations where observed data have been collected during the year 2000, the year used for FVCOM model validation in this work, are presented in Table 1 and mapped in Fig. 3. Note that the inland daily observations are difficult to interpret. BIWTA and the Bangladesh Water Development Board (BWDB) measure data twice a day for the daily high and daily low tide level. However, data collection normally occurs during daylight hours from 6am to 6pm. Chatterjee et al. (2013) describe the approach of measuring water levels with the Visual Tide Staff (VTS) or Tide Pole. The confidence level in this manual data is not very high, and although data readers have long experience of collecting these data sets there is little metadata on method and timing of observations, together with an uncertainty in the water level datum.
The tidal type can be classified using (3) using K 1 ¼ K 1 amplitude etc. If F is less than 0.25, the tide is classified as semidiurnal; if the ratio is from 0.25 to 1.5, the tide is classified as mixed, mainly semidiurnal; if the ratio is from 1.5 to 3.0, the tide is classified as mixed, mainly diurnal; if the ratio is greater than 3.0, the tide is classified as diurnal (Defant, 1961;Howarth and Pugh, 1983;Pugh, 1987;Pugh and Woodworth, 2014). The final two columns of Table 2 contain observed and modelled tidal type at the 19 sites. Throughout the model, the tide is classified as semi-diurnal at most sites, except for Sundarikota, and Isla Ghat in mouth of the lower Meghna. Here, the presence of freshwater changes the behaviour of the tide throughout the year and the tidal type is classified as mixed. Where the seasonal cycle is particularly strong the phase or amplitude are marked as 'v' because the behaviour of the tide is annually variable. Tidal range can be classified as micro-(< 2 m spring range), meso-(2e4 m), or macro-(> 4 m) (Davies and Moses, 1964). The western delta around the Sundarbarns is macro-tidal, while the central and eastern section are meso-tidal, classifications derived from the model are largely in agreement with those from the UK Hydrographic Office data, though the highest tidal ranges are slightly underpredicted.

Open coast
To assess the model performance, water levels were compared with UKHO and tide gauge data at several sites throughout the delta. Tables 2 and 3 compare modelled amplitude and phase of the four leading tidal constituents with UKHO chart data at 19 points around the Bay of Bengal. Some summary statistics are presented in Table 4. The statistical measures used are the root-mean-squared error (rmse), coefficient of correlation (R 2 ) and percentage model bias (Pbias), defined as Pbias ¼ 100 where M represents the model prediction, D represents the measured data, and N the total number of data points used. The rmse presents an absolute error for the model data, R 2 is an indicator of how much of the variance is explained by the correlation, and Pbias provides a measure of whether the model is systematically over-or under-predicting the measured data. The modelled tidal amplitudes are in good agreement with the chart data with small rmse of 0.11 m on average, and an average Pbias of below 10%. The dominant semi-diurnal constituents are particularly well represented, while the amplitudes of diurnal constituents tend to be slightly over-predicted and are biased high by around 10%. The diurnal constituents are also less well correlated, however they are of a much smaller magnitude. These components are captured more poorly in the model, with larger rmse and Pbias in their phases. We also compare the FVCOM model performance with that of Rose and Bhaskaran (2015), who report an above 88% correlation coefficient between predicted and observed water levels at the same set of tide gauges, using a harmonic analysis reconstruction approach. Analysing our reconstructed time series from the FVCOM model using the same approach yields comparable correlation coefficients, as listed in Table 5. There is poor agreement at the Hatia Bar site, where there is too much energy in the modelled semidiurnal constituents in FVCOM. It is possible that this is a pinchpoint for the model, where incorrect representation of the bathymetry is limiting the tidal flow. Fig. 2 shows that there is a gap in our knowledge of water depths at the mouth of the Meghna, which  are not covered by Admiralty charts, or river cross sections. A full tidal harmonic analysis including over 100 constituents is also available at five gauges (Hibbert, personal communication). These data were collected over a long period, and some sites represent tidal behaviour up to 120 km inland from the coast: Chandpur (1979e1980), Sundarikota (1980), Ramdaspur (1969), Barisal (1980), and Mongla (1976). The leading 65 constituents were extracted from the model, and a linear regression performed. The amplitudes displayed R 2 ¼ 0.84 and rmse ¼ 0.091 m. The phases exhibited an R 2 correlation of 0.63 and an rmse of 11.59 . Overall this is a very positive validation of the performance of FVCOM at the coastal sites of the GBM delta. Next, the ability of the model to accurately represent water levels within river channels is assessed.

Within the delta
With confidence in the performance of FVCOM near the coast, the model was next evaluated further inland within the complex river channels of the delta. Fig. 4 compares modelled mean tidal range with direct observations made by BIWTA. In the central section of the delta the modelled mean tidal ranges are in good agreement with the observations, varying between 0.7 and 2 m. There are some discrepancies between the observed and modelled ranges however. The model seems to give reduced mean tidal ranges in the entrance to the lower Meghna, and within the Sundarban region to the west. Looking at the map of model results (right hand panel of Fig. 4) it can be seen that large ranges at the open coast are dropping away as the tide enters river channels. This may be due to the bathymetry of restricted river channels at these sites. Fig. 5 shows the model bathymetry on the delta. The two problem sites, where modelled tidal range seems low, coincide with areas of particularly deep bathymetry in excess of 15 m in the model (present in both GEBCO, and newly surveyed cross-sectional data). Wang et al. (2009) discuss the importance of accurate bathymetry for intertidal areas, claiming it only of secondary importance to tidal forcing in terms of inputs to a tidal model.
As there is no absolute datum at the observed sites, direct comparisons of absolute water level must be treated with some reservations. Nonetheless we are able to perform some statistical analysis to evaluate model performance inland. Table 6 shows rmse and R 2 defined as for the open coast sites. Pbias, will not be considered in this case, as the water levels measurements are all relative. The model and observations are well correlated at Babuganj, Barisal Kaitpara, Jhalakati, and Sureswar. These are all river dominated sites, with a large annual cycle driving water levels. Using this measure, the model performs poorly at Rayenda, Bamna, Umedpur, Chalna, and Tajumuddin, however relative measures such as the tidal range, are still well captured. Fig. 6 shows modelled and observed water levels at Rayenda and Sureswar. Both time series show twice daily observations of water levels (black crosses), and hourly model output (continuous red line). As well as a strong diurnal cycle, a secondary 14 day spring/ neap cycle is superimposed, and an annual cycle is also observed: demonstrating the effect of the volume of incoming freshwater on  the surface elevation. Rayenda is seen to be a tidally dominated site, with little change observed in the model when moving from wet to dry season. There is some freshwater signal in the observationsleading to a poor correlation in Table 6. At Sureswar, the freshwater dominates the water levels in both model and observations. The annual cycle is strong, with an amplitude of around 2 m. At Sureswar and other sites (not shown) there is a larger tidal range during dry season. Peaks in discharge are also captured as changes in water level: compare with the lower panel of Fig. 6. Particularly obvious are the two events in October/ November 2000, and to a lesser extent in June. During the wet season, the tidal range is significantly suppressed. Here the tidal wave has to propagate up the estuary in direct competition with the river flow, so the time taken for the tide to propagate is longer. Fig. 7 summarises the relationship between observed and modelled daily high and low tides within the delta; a year's worth of daily data is shown for the 13 sites listed in Table 6. The scatter  plot is coloured by a normalised density. There are two 'hot spots', representing the daily high and low waters. This composite result suggests that FVCOM is overpredicting water levels at low tides (perhaps because the model is drying out, but a negative water level is reported in the observations). The high tides are better correlated, with the majority of high tides being correctly predicted. The slope of the scatter plots suggests that we are generally underpredicting the water levels. Fig. 8 shows the change in tidal range moving from the estuary mouth to the head for 2 sections. The 'eastern' section covers the lower Meghna from Hatiya Bar in the south to Sureswar in the north, and the 'central' section has its mouth at Chardoani, and its head at Babuganj. Both estuaries see the tidal range reduce as the wave penetrates inland. The tidal ranges are larger in the observations, but the dissipation rates (slope of the linear fit) are similar between model and observations. The right hand panel of Fig. 8 shows the phase lag of the tide as it travels up the modelled estuaries. A wave travelling at 160 km in 9 h has a speed of around 4.3 m À1 , closer to the propagation speed of a storm surge simulated by Karim and Mimura (2008). Unfortunately it is not possible to recreate this plot for the observed data, as they are not recorded at a high enough frequency. Fig. 9 shows the co-tidal chart for M2. The tidal amplitude is largest in the North East and North West corners of the Bay of Bengal, consistent with observations e.g. Murty and Henry (1983). The tidal wave propagates inland, with the amplitude decreasing as it does so. The pattern of the S2 tide (not shown) is similar with the phases offset by around 150 .

Friction
In order to determine what controls tidal penetration into a river delta, we can perform a series of model experiments. First, we examine the effect of changing bed roughness in the model. Fig. 10 shows the effect of altering bottom roughness lengths (z 0 ) on tidal range. A suite of tests was run, forced by the same tides and with no river discharge included; z 0 was set to 0.1 cm, 0.15 cm, 0.5 cm, 1.0 cm, 3.0 cm, 5.0 cm, 10 cm, and 50 cm (the default roughness length was 1.5 cm). Their equivalent C d values are calculated in Table 7.
With reduced roughness length, the drag is reduced and tidal range increases. Similarly when the roughness length is increased, Fig. 6. Top panels: Time series plots comparing the annual cycle of water level in BIWTA observations with FVCOM model at Rayenda (top) and Sureswar (below). Red lines are modelled water levels, and black crosses are observations taken at 6am and 6pm daily. The lower panel shows river discharge during the same period. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) more tidal energy is dissipated, and the range reduced. Changing the roughness length by a factor of 10 changes the mean tidal range by 20 cm across all sites. The maximum change in tidal range can be as large as 0.5 m, and the largest changes are seen in areas of high tidal range. Even so, we cannot increase the amplitude of modelled tidal range to that observed merely by altering bed roughness within a reasonable range of lengthscales.

Competition between river and tide
The GBM delta is situated at the interface between large and annually varying river discharge and a meso-to macro-tidal range. However the time scales of these two processes are very different, with the tide dominated by the semi-diurnal (12.4 h) cycle, and the freshwater following a seasonal cycle). Across the delta there are locations where the tide dominates, and others where the river dominates. Some of the most interesting places are where there is competition and the dominance switches during the year.
We can isolate the effect of river flow, by taking an averaged model current speed over a tidal period (or longer). Fig. 11 shows the average current speed from a dry and wet lunar month (29.5 days) respectively. By averaging out the tidal flow, we see small residual current speeds in the open ocean, but much larger values in the river channels. In the dry month, the fast flows are restricted to the narrow and deep sections of the major rivers, most clearly seen in the Upper Meghna. In the wet season, the residual flows are larger, and cover more of the river system. There are fast flows in the Arial Khan and throughout the Meghna. The fast currents are  also visible across more of the width of the river channels, showing that the shoulders of the rivers have been inundated, and allowed to become active 'wet cells' in the model.
The time taken for the tidal wave to propagate into the delta is measured by comparing the arrival time of high tide with the time of high tide at the shelf break (defined at 200 m water depth). The tidal wave travels inland up large river channels (such as the Meghna) at a speed of ffiffiffiffiffi ffi gh p þ U where h is the water depth, and U is the additional velocity provided by the river discharge. An estimation of average speed of the tidal wave up the delta suggests it takes approximately 17 h to penetrate 440 km inland to the head of the river, which equates to a propagation speed of 7 m À1 , or a tidal wave moving in an average water depth of 5.3 m. This speed is quicker than, but comparable to, that seen in the surge propagation by Karim and Mimura (2008). However, in some of the minor river channels a more complex travel pattern is observed, as river channels of different sizes have different hydraulic depths. For example, in a confluence of channels around 90.2 E and 22.5 N we see that one wave is propagating down from the North with a lag of 12 h, which is then meeting the tide approaching from the south whose peak arrive around 5 h earlier (not shown). The hydraulic mean depth (H) (dependent on cross sectional shape) is commonly used for 1D channel flow (where the hydrodynamics cannot be properly captured) and is a measure of channel flow efficiency. For a given channel cross-section, H is defined as the cross section of water flowing through the channel divided by the wetted perimeter of the channel. Thus it can be regarded as the equivalent depth for the same cross-section in a rectangular channel and takes account of varying depths across a cross-section. Our FVCOM model is 3D, so the full river cross-sections are represented, and this simplified 1D approach is not required. This idealised relationship will also break down as we move into tidally controlled estuaries. Nonetheless it is a useful way to compare the efficiency and connectivity of channels. Fig. 12 shows the impact of freshwater on the arrival time of the high water, by running the identical model grid and tidal forcing with, and without river discharge. In the model without rivers, the high tide is seen to propagate all the way up to the head of the river, where the FVCOM domain ends. However when freshwater forcing is enabled, we see that the tide is only able to penetrate to confluence of the Brahmaputra and Padma (23.79 N,89.76 E), and no further inland as the upper Meghna and Gorai rivers are dominated by freshwater. The river also affects the arrival time of the high tide further towards the mouth, by as much as 1 h in places.  . 11. Difference in residual current speed between a dry (top) and wet (below) season month.

Bathymetry and channel width
A consistent model grid and bathymetry were used across all experiments presented here, so it is not possible to run sensitivity tests to changing bathymetry or coastline. However, as the delta model configuration contains a wide range of river channels, the tidal propagation in contrasting channels can be investigated. Consider the arrival time of high tide in the Meghna (around 14 km wide) versus the Rupsha (around 3 km wide). The speed of tidal propagation was calculated at two sections; at 89.6 E through the Rupsha, and at 90.6 E through the Meghna. In the wider channel, the tidal wave travels at around 9.8 m À1 , while in the narrow channel the propagation is much slower: around 4.8 m À1 .
The water depth in the narrow channel (17 m on average) exceeds that of the wide channel (9 m on average). If the wave speed were simply ffiffiffiffiffi ffi gh p , then one would expect faster propagation in the deeper of the two. This is not the case, due to lateral friction at the channel walls. In a wide channel the tidal wave will propagate freely close to the shallow water wave speed. But in a narrow channel the flow and wave speed will be constrained by frictional loss to the walls. This effect has been described by Stoker (1948), and previously Thomas (1937) who present the theory behind wave propagation in open channels. In hydraulics, the frictional effects at the river walls are represented by a "loss of head" term. In turbulent flow this resistance term is proportional to u 2 .
Turning to wave speed Savenije (2005)  proportional to the length scale of the exponential width variation, in other words the distance over which the channel converges. So a sharply funnelling channel will damp the tide much more efficiently than one which narrows over a longer distance. The width profiles of most individual channels can be reasonably approximated with an exponential curve (Davies and Woodroffe, 2010), and the Rupsha and Meghna are no exceptions. The convergence length (L b ) is defined distance over which the channel width narrows by a factor e. For the Meghna L b is approximately 55 km. While for the Rupsha L b z 13 km. The smaller value of L b for the Rupsha compounds the observations that channel width is controlling tidal penetration into this channel.

Discussion
What controls tidal penetration in a river delta? There are four major controlling factors; (1) bathymetry (2) channel width (3) opposing river discharge (4) bed friction. Bathymetry is the most fundamental of these, and is the most important control on tidal penetration. Factors 1,2, and 4 can be considered as static and predetermined by the geometry (and must be captured by the model set-up). River discharge is different as it is externally controlled and varies seasonally due to the monsoon. We can investigate the effects of bathymetry and hydraulic depth by examining maps of tidal range and arrival time within different channels in the delta. The effect of discharge can be assessed by a model experiment with and without river discharge, and sensitivity experiments to subgridscale effects have shown the impact of friction on tidal penetration.
Model bathymetry is one of the most difficult data sets to obtain (Lewis et al., 2013), but also the most important to the hydrodynamic behaviour of estuaries. Several data sources were available in this case, GEBCO and ETOPO offshore, and river channel surveys inland. These data are patchy in both space and time, and also raises the question of "how long are these observations really valid for?" The 'lifetime' of bathymetry data in this dynamic environment is, likely, very short, due to the soft underlying sediment combined with high river and tidal flows. The shape of river cross-sections may be important to flow locally, but when considering tidal penetration, the channel area and river volume will be the controlling factors.
The true position of the coast is further complicated by manmade defences, by poldering and building sluice gates so that stretches of river are made more stable artificially. Outside of the Sundarban forest, many areas are protected by sea defences, as mapped by Islam (2006). However, we are not attempting to predict morphological evolution in this model, so for the purpose of this work it is assumed that the position of the coastline remains constant, with no inundation permitted beyond a predefined channel wall. Land inundation has been addressed in the separate work of Lewis et al. (2013) and Haque et al. (2016). Errors in the model bathymetry are a well known problem in coastal modelling, e.g. French and Clifford (2000) suggests that there should be a feedback loop between the hydrodynamic model run, and bathymetry generation. Cea and French (2012) perform 250 Monte Carlo simulations of the tidal flow in a meso-tidal estuary. Their results indicate that adjusting for bathymetric error may be a more effective and appropriate parameter than bed friction for the calibration of depth-averaged hydrodynamic models. In their case study estuary, tidal current speed is shown to be much more sensitive to small changes in bed elevation than to adjustment of the bed friction coefficient.
The impact of rivers on tidal penetration was considered in section 4. The tidal wave speed is modified by opposing velocity from river discharge. Where the river current speed is strong, the tidal wave is retarded, and the tide cannot penetrate so far inland. This is reflected when examining the tidal range at some river dominated sites (e.g. consider Fig. 6). Where the site has a high river discharge, such as Sureswar the tidal wave is retarded by outflowing river discharge. As the tidal wave dissipates as it propagates up the river channel, this slowing of the tidal wave is manifested as a suppression of the tidal range at the high discharge sites.
Model sensitivity experiments to changing bottom friction showed that the tidal range is sensitive to roughness length (z 0 ). A spatially homogeneous roughness length was used throughout the model. In reality, there is likely to be larger dissipation in river channels, due to the presence of river banks, and vegetation. Li et al. (2012) use FVCOM to simulate changing tides in Darwin Harbour, noting that tidal asymmetry is controlled by an enhanced bottom friction within the mangrove forest areas. This would suggest that to accurately capture the tides in the GBM delta, we would need an area of larger z 0 within the Sundarban forest. Contrary to this, in our model we see that the tidal range is already slightly underpredicted in the delta, suggesting friction is overly strong in the river channels.

Summary
The FVCOM model has been configured and implemented for the Bangladesh delta. Having combined modern river cross-section surveys with GEBCO data and satellite observations, a new bathymetry data set has been produced. This bathymetry was used to generate a grid for the unstructured model, which was then driven with tidal and freshwater forcings.
The modelled sea surface elevations have been compared against tidal data, and water level observations inland. At the open coast, where tides are dominant, the model performs very well, giving us confidence in its abilities. Within the delta FVCOM was able to accurately model the tidal ranges in the central estuarine region, however the model suffered from some under-prediction of amplitude in the Sundarban forest region, and within the lower Meghna. Correct representation of the channels, and bathymetry and bottom friction is crucial, and problems in model bathymetry can be identified by evaluating the behaviour of the tidal wave. Channel width can also control tidal penetration, with rapidly narrowing estuaries seeing more dissipation than wider rivers.
Water levels in the delta are controlled by a balance between river and tidal flow, acting on different timescales. Throughout the year the situation can change; from tides controlling the water levels in the dry season, to a dominance by river flow during the monsoon. Very large flows in the monsoon season are shown to reduced tidal amplitudes by retarding the tidal wave speed. Seasonal water level changes associated with freshwater input can be significantly larger than the total tidal range, and the model is also seen to 'dry-out' where modelled bathymetry is shallower than in reality. A sensitivity to model bottom roughness length z 0 is seen, which can change the tidal range locally by 20 cm on average. However the magnitude of this effect is insufficient to overcome an underprediction of tidal ranges which is thought to be constrained by poor model bathymetry. The delta model developed has been shown to be suitable for predicting tidal behaviour far inland in the delta as well as capturing the impacts of river flow on water levels.