Shoreline Change around a River Delta on the Cox’s Bazar Coast of Bangladesh

Abstract

A recent erosional problem around a river delta on the Cox’s Bazar coast was analyzed in this study. The coastline extends from south to north. Rapid erosion has affected some portions of a 24-km road along the coast, and local authorities have attempted to protect the road via revetment. However, the structure was soon buried with sediment because of a growing sand spit along the river delta, and a new area was eroded. Shoreline positions for a 44-year (1972–2016) period were digitized using Landsat images. From the time stack images, we observed a sand spit growing in a northward direction from 2000 to 2015, and the adjacent erosion area extended in the same direction. We employed a numerical model (MIKE21FM SM) for the computation of wave-driven currents and sediment transport along the coast, and attempted to reproduce recent erosional processes. The numerical result shows that net littoral drift is dominant in the northward direction along the coast, which is the same direction of the spit growth observed in the satellite images. A higher amplitude spit induces higher sediment transport compared to a low amplitude spit because of the difference in local incident wave angles resulting in greater positive gradient of the longshore sediment flux distribution, causing erosion in the downcoast.

1. Introduction

1.1. Aim of the Study

River deltas are extremely dynamic due to the complex interaction between fluvial sediment discharge and nearshore wave climate [1,2]. Deltaic shoreline changes continuously due to natural hazards (like as storm surge, tsunamis, coastal flooding, coastal erosion and sea level rise) and anthropogenic activities (such as river damming). These coastal phenomena occur on the scale of week to decades. Thus, monitoring of these shorelines is essential for understanding the morphological behavior of the nearshore environment.

Cox’s Bazar is the main tourist destination in Bangladesh to enjoy a natural sandy beach, and the government would like to promote this area and build new infrastructure. The overall coast seems stable over the long term. However, recently, the sandy beach has been subject to severe erosion north of the Reju River delta. Rapid erosion has affected some portions of a coastal road along the coast, and local authorities have attempted to protect the road via revetment. However, the structure was soon buried with sediment because of a growing sand spit along the river delta, and a new area was eroded. We speculate that this accumulation was because of sand spit growth north of the Reju River delta, which may have changed longshore sediment transport causing erosion in the down-drift area. Thus, we sought to verify our idea and to understand the morphological processes in this area, including:

• detection of the shoreline position using Landsat images for the period 1972–2016 to understand the morphological behavior of the southeast coast of Bangladesh, focusing on the Reju River mouth delta area;
• identification of the causes of the recent erosion problem north of the Reju River delta using a numerical model.

Remote sensing monitoring, such as Landsat images [3], synthetic aperture radar (SAR) images [4], video images [5], and the global positioning system (GPS) [6], have been widely used for shoreline delineation and change quantification. Among various images, 30 m resolution of Landsat images are cost-effective and suitable for monitoring of sandy areas with particular characteristics (deltas, very large beaches with large variations in time) [7,8]. Thus, we acquired Landsat images for shoreline delineation in the present study.

1.2. Previous Studies

Spits are very dynamic which can be formed in deltas. Sahalin spit [9] and the Trabucador La Banya spit [10] are examples of spits which developed in the Danube and Ebro River deltas, respectively. Rapid growth and migration rates of these spits are mainly due to a large amount of sediment transport in both the along- and cross-shore directions. Kraus [11] explained the governing processes of spit evolution and particularly emphasized the importance of overwash processes. Moreover, he discussed a primary factor for determining the overwash intensity which is dependent on the wave climate, particularly the wave period.

The temporal and spatial variability of the Guadalfeo River mouth was studied by Bergillos et al. [12] with the help of extensive field measurement and numerical models. Their results showed that the river damming led to shoreline retreat and bed-level erosion due to a lack of sediment supply from the river. Recently, they proposed an integrated tool to predict the morphodynamic response over the long-term for the deltaic coasts (Playa Granada, Spain) [13] and also investigated the storm response for this coast under varying wave directions using numerical models [14]. The model results recap the importance of cross-shore and longshore sediment transport in driving the coastal storm response at this location.

Petersen et al. [15] discussed the evolution of spits, suggesting they are mainly driven by gradients in alongshore sediment transport. The study was based on a simple analytical model and a two-dimensional numerical model, and were supported by experimental results. The key findings of the study were that a spit is nearly at equilibrium when constant waves approach the spit at angles greater than 45° and the width of the spit is propositional to the width of the surf zone.

Thomas et al. [16] explained the long-term evolution of a sand spit at Ginst Spit, West Wales, using aerial photographs and topographic surveys. A varied correlation was found between the annually-averaged wave components derived from a wave model and the sand spit rates of the shoreline change. In addition, rainfall and growth of the spit were associated with flood events, which suggested that lowland inundation was caused by the combined contribution of fluvial and coastal processes. Recently, Nienhuis et al. [17] explained that river mouth spit development and spit breaching phenomena with assistance from a numerical model at wave-dominated river mouths and successfully concluded that both wave climate and fluvial discharge play major roles in spit evolution.

From the aforementioned literature review, we can understand spit evolution in different areas of the world and influential factors. Based on this knowledge, we attempted to explain the morphological behavior around the spit at the Reju River delta to assist in developing an efficient countermeasure against erosion, which is presently unknown in this area. In this context, this study describes spit morphology and how it affects downcoast areas with the aid of remote-sensing data and a numerical model.

2. Study Area

The study area is on the southeastern coast of Bangladesh, Cox’s Bazar (Figure 1a), which is the longest uninterrupted natural sandy beach in the world [18,19]. It is bounded on the west by the Bay of Bengal and the north by the Moheskhali Channel. There are number of seasonal and intermittent streams that serve as the natural drainage in this area [20]. The Reju River is among the important channels and is at approximately x = 50 km, where x is the longshore axis defined in Figure 1a. The major commercial part of the beach is between x = 64 km and x = 67 km, and is the country’s most attractive tourist destination. The overall shoreline seems stable; however, erosion has caused problems maintaining a coastal road near the Reju River mouth area, which is the focus of this study (Figure 1b).

A 24-km-long road, Marine Drive Road (Figure 1a between x = 43 km and x = 63.7 km), was built along the coast to promote tourism opportunities, facilitate the fishing industry, exploit more natural resources, and establish a portion of a regional highway. Road construction began in 1993–1994 and was completed in 2008. Soon after completion, the road has been subject to severe erosion near the Reju River delta. The fluctuation in the shoreline in front of the road section was minimal and the area was not vulnerable to erosion prior to construction: the road was placed with a setback of approximately 500 m from the shore. Soon erosion became apparent, and the shore width off the road decreased to 50 m. The road construction authority has attempted to protect the road via temporary protection work (sand geo-bags) as well as a hard structure (a revetment). The revetment (400 m in length) was built between x = 53.7 km to 54.1 km during 2009–2010 (Figure 1b). However, soon after the construction of the revetment, sand accumulated in front of this structure and the shore widened again to approximately 400 m, and another eroded area appeared further to the north.

The tidal characteristic of this area is semi-diurnal, and the tidal range is approximately 3 m during spring and 1 m during neap tide. The climate of Cox’s Bazar includes four seasons: dry (December–February), pre-monsoon (March–May), monsoon (June–September), and post-monsoon (October–November) [21,22]. Most of the cyclones approach the area during the pre-monsoon and post-monsoon seasons. Waves are dominant mainly during the monsoon season, originating from the southwest direction.

3. Data

3.1. Shoreline Data

In the study area, 26 cloud-free scenes of Landsat imagery (

Table 1. List of the Landsat images used in this study.

Serial No.	Satellite	Sensor	Date
1	LANDSAT_1	MSS	1972/12/27
2	LANDSAT_2	MSS	1976/01/26
3	LANDSAT_2	MSS	1978/01/15
4	LANDSAT_3	MSS	1980/01/14
5	LANDSAT_5	TM	1988/02/28
6	LANDSAT_5	TM	1989/01/13
7	LANDSAT_5	TM	1990/03/05
8	LANDSAT_5	TM	1991/11/19
9	LANDSAT_5	TM	1993/11/08
10	LANDSAT_5	TM	1994/12/13
11	LANDSAT_5	TM	1996/02/18
12	LANDSAT_5	TM	1997/11/03
13	LANDSAT_5	TM	1998/12/24
14	LANDSAT_7	ETM	1999/12/19
15	LANDSAT_7	ETM	2000/12/21
16	LANDSAT_7	ETM	2002/12/27
17	LANDSAT_7	ETM	2004/12/08
18	LANDSAT_7	ETM	2005/11/25
19	LANDSAT_7	ETM	2006/12/30
20	LANDSAT_7	ETM	2008/11/17
21	LANDSAT_7	ETM	2010/01/23
22	LANDSAT_5	TM	2011/02/11
23	LANDSAT_7	ETM	2012/12/06
24	LANDSAT_8	OLI_TIRS	2013/12/17
25	LANDSAT_8	OLI_TIRS	2015/03/01
26	LANDSAT_8	OLI_TIRS	2016/10/22

) were acquired for the period of 1972–2016 from the Earth Explorer database of the U.S. Geological Survey [23]. The time of acquisition of the images was mainly during the dry season from November to February as there were no clear images during other seasons. There are several methods for shoreline detection using Landsat imagery. Kuleli [24] showed that the mid- and near-infrared spectral bands help to separate land from water, which distinguishes the shoreline position using satellite images. Benny [25] successfully extracted the shoreline using band 7 of Landsat Multispectral Scanner (MSS). In addition, Alesheikh et al. [26], Frazier and Page [27], and Sarwar and Woodroffe [28] have used band 5 of Landsat Thematic Mapper (TM) and Enhanced Thematic Mapper (ETM) for shoreline delineation. Recently, Josep et al. [29] extracted the shoreline from the infrared bands of Landsat Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS). Thus, both the near-infrared and infrared bands were adopted in this study. The separation of the land surface from the water bodies can be easily identified via visual inspection. The entire shoreline was digitized manually using the images. Moreover, a tidal correction was completed for the entire shoreline assuming the slope of the beach area was 1:80. Shore positions relative to a fixed shore-parallel baseline (x-axis in Figure 1a) at 700 transects at a 100-m interval were extracted. All the extracted shoreline is presented in Figure 2, in which the red stripe indicates the erosion at Marine Drive Road.

3.2. Wave and Rainfall Data

Wave hindcast data for the period between 1988 and 2014 were downloaded at a grid point (latitude 92°; longitude 21.25°) from the European Center for Medium-Range Weather Forecasts (ECMWF) ERA-Interim archive [30]. The wave data includes significant wave height, mean wave period, and direction at 6-h intervals as shown in Figure 3. As previously mentioned, higher waves were found mainly during the monsoon season mostly originating from the southwest direction (Figure 4). Furthermore, daily rainfall data was collected from the Bangladesh Metrological Department (BMD) during 1988–2013. Figure 5 shows the daily rainfall at Cox’s Bazar station; the annual average rainfall was 3832 mm during 1988–2013. Most of the rainfall occurs during the monsoon season.

4. Morphological Analysis

4.1. Shoreline and Image Interpretation

We calculated the statistical properties of the shoreline variation including the mean range and standard deviation in which the range was defined as the envelope of the shoreline positions occupied at each profile location over the duration of the study. Small standard deviations and narrow ranges indicate stable regions, while higher standard deviations and wide envelopes indicate high variability regions. The mean shoreline profile with standard deviation is presented in Figure 6, showing relatively large shoreline changes at the northern end (x = 70 km), as well as around the river mouth of the Reju River (x = 50 km). The range around the Reju River area is approximately 500 m, and the standard deviation is approximately 120 m. Variations in the shoreline for the remainder of the region (x < 40 km) are relatively small, indicating stable beach conditions. Moreover, we also calculated the percentiles of the shoreline, which is an efficient statistical property for check variability of the datasets [31]. The 10th and 90th percentiles with respect to mean shoreline are shown in Figure 6. Variability is relatively large at the northern end (x = 70 km), as well as around the river mouth of the Reju River (x = 50 km).

Figure 7 shows the time stack of Landsat images near the Reju River mouth during 1972–2015. As it can be seen, there is a significant morphological change at the Reju River mouth during 1996–2000: by 2000, remarkable sand deposition is observed at the river mouth (red circle). Following the period 2000–2013, the sand deposition extends northward and is transformed into a sand spit. The coast in front of Marine Drive Road (53.7 < x < 55.8 km) was severely eroded during 2006–2008. Then, sediment accumulated in this area until 2011 because of the growing sand spit. This accumulation also covered the revetment which was constructed during 2009–2010. However, another eroded area emerged further north of the sand spit.

The shore width in front of Marine Drive Road was measured from the shoreline to better understand the progress of erosion in the roadside area. Here we chose a section (x = 54 km) where the road erosion started; its variation is shown in Figure 8. When the road construction started in 1993, the shore width was approximately 500 m. The road construction authority may have considered the width sufficient to safely maintain the road at this time. However, following 2002, the beach width narrowed considerably and reached a minimum (approximately 49 m) during 2008. A revetment was built during 2008–2009 for protection of the road. However, a growing spit covered the revetment in 2010, and the shore width recovered to approximately 500 m. The revetment work in this area for road protection was not necessary, or excessive; low-cost temporary protection would have been sufficient. Moreover, the shore width at a northern point (x = 56.8 km) was not affected by the erosion in 2008; however, this area was eroded in 2010 (Figure 8). The shore width of this area was approximately 160 m during road construction in 1993; however, the width decreased and reached a minimum (approximately 30 m) in 2010. The revetment had no influence on the erosion at this section because it was buried by a growing spit soon after construction.

4.2. Variation in the Area of the Delta

We calculated the area of the delta around the Reju River mouth to trace the long-term morphological change, which may have a connection to the erosional problem of Marine Drive Road, and attempted to identify influential factors that affect the growth and decay of the delta.

The area of the delta was calculated between the distances of x = 49.7 km and x = 55.2 km, which covers the most dynamic region given the greater standard deviation along this stretch (Figure 6). Northern (x = 49.7 km) and southern (x = 55.2 km) boundaries were chosen considering these are stable (or less variable) locations over a long-time period. The eastern boundary was set as Marine Drive Road (Figure 9). As the Reju River flows through the delta, we separated the delta into two portions: a northern delta and a southern delta.

Figure 10 shows the variation in the area of the delta. The variation in the total area of the delta (the sum of the area of the northern and southern deltas) shows some periodic change over a long-time period. The minimum area of the northern delta was observed in 2008. However, the area of this delta dramatically increased during 2008–2010, perhaps caused by sedimentation in front of the newly constructed revetment. The variation in the area of the southern delta showed an opposite pattern compared to the northern delta, indicating the area of the southern delta increases with a decrease in the area of the northern delta. Thus, a negative correlation (R = −0.72) was found between the area of the northern and southern deltas. The area of the deltas (northern, southern, and total) changes with time because of sediment supply from the river, longshore sediment transport, and other factors (storms, waves) as discussed in the following sections:

4.2.1. Influence of the River Mouth Position on the Delta

The Reju River mouth position may influence the area of both deltas (north and south delta). We identified the river mouth position from the images during 1972–2016. Figure 11 shows the variation in the river mouth position. The river mouth moved northward for most of the periods, but abruptly changed to southward in 1976, 1991, 1996, and 2010.

From the Landsat images, we observed several breachings during 1972–2016. Breaching can be defined as a shift in the position of the river mouth as a result of erosion, as shown in Figure 12. Breaching can occur during catastrophic storms, and a new inlet that provides a shorter route for tidal exchange normally stays open while the less efficient old inlet gradually closes [32]. From the historical storm surge record, several cyclones struck this coast in 1991, 1995, 1996, and 1997. Breachings in 1991 and 1996 may have occurred because of these cyclones. The river mouth moved to southward mainly because of breaching through the southern delta; part of this delta attached to the northern delta and likely increased its area.

The areas of both deltas (northern and southern) change in accordance with the position of the Reju River mouth. The correlation between the area of the two deltas and the river mouth position is shown in Figure 13. There is a good correlation (R = 0.95) between the area of the southern delta and the river mouth position; the southern area increases if the river mouth position moves northward. However, the reverse phenomenon was found for the northern delta with relatively low correlation (R = −0.74); the area of the northern delta increases when the river mouth position shifts southward.

4.2.2. Influence of Rainfall on the Delta

As the amount of river sediment flux was not measured, rainfall data was taken as a proxy of river sediment input. Twenty images were available for the period of 1988–2013 that matched the rainfall collection period. The change rates of the area of the delta were calculated between two consecutive periods, along with the average rainfall for the same period.

The change rate of the area for total, northern, and southern delta with average rainfall is shown in Figure 14. There is a negative correlation (R = −0.4) between rainfall and the total delta area change rate which indicates higher rainfall reduces the area of the delta. We found a similar negative correlation (R = −0.36) for the area of the northern delta; however, there was nearly no correlation (R = 0.07) with rainfall for the southern delta. Prior to the analysis, we thought that higher rainfall would result in a greater sediment load that might increase the area of the delta. However, this was not the case: higher rainfall produces a higher flow volume in the river flushing sediment further offshore area and not contributing immediately to the surrounding delta. Alam et al. [33] explained some sediment characteristics of different geomorphic units in the northern portion of the Cox’s Bazar area very close to our study area. They summarized that the sediment size of this area is very fine sand to silt. If the river sediment is very fine, it is suspended in the water column and is flushed offshore via high river flow.

4.2.3. Influence of Offshore Waves on the Delta

Among the simplest and still useful deep-water wave parameters that includes the wave period is the wave steepness, H/L [34], where H and L are the deep-water wave height and length, respectively. In the early 1930s, it was known that steeper waves, usually observed during storms, were responsible for beach erosion, while mild conditions are favorable to beach accretion. The offshore wave height and period were used to calculate the dimensionless wave steepness parameter during 1988–2014. Furthermore, the wave steepness was averaged per day using the same period as the area of the delta calculation.

The averaged deep-water wave steepness plotted with the change rate of the area of the delta is shown in Figure 15. We found a negative correlation (R = −0.20) between average wave steepness and the change rate of the area of the total delta, and similar negative correlation (R = −0.22) was found for the northern delta, indicating steeper waves reduced the area of the delta. Steeper waves during high energy events may have eroded the delta, and the area of the total delta reduced with the increasing wave steepness and vice versa. However, there was nearly no correlation (R = 0.10) for the southern delta.

A discussion of the variation in area of the delta is summarized as follows: The area of the delta changes in response to the position of the river mouth, rainfall, and waves. Their correlations are summarized in

Table 2. Summary of correlation analysis between area of deltas and different influential factors.

Deltas	Position of River Mouth	Rainfall	Wave Steepness
Northern	−0.74	−0.36	−0.22
Southern	0.95	0.07	0.10
Total	−0.05	−0.40	−0.20

. Moreover, the variation in the total area of the delta shows some periodic change over the long study period, but it was difficult to detect any trend over the long term. Thus, it was not possible to connect the shoreline change downcoast of the delta with the area of the delta. However, we observed a sand spit growing northward in the northern delta, which may have a relationship with downcoast erosion. Thus, we defined the amplitude of the sand spit and analyzed its variation in connection with the shoreline change downcoast, as described in the next section.

4.3. Variation in the Amplitude of the Sand Spit

To measure the amplitude of the sand spit, we defined the tip and toe of the spit. The tip is the peak (maximum offshore position) of the sand spit, and the toe is the minimum cross-shore position of the sand spit that is connected to the adjacent coastline. The toe position may change between longshore and cross-shore depending on the tip position of the spit. Then the amplitude (A) was defined as the perpendicular distance from the toe to the tip, as shown in Figure 16.

The maximum erosion north of the sand spit was observed in 2008, and the erosion extended approximately 2 km from the toe as measured along the x axis shown in Figure 1a. This result suggests the growth of the sand spit may have affected 2 km of shoreline from the toe downcoast. The spatial mean of shoreline change was calculated between the subsequent years from the toe to 2 km northward and compared to the amplitude of the spit.

The variation in the amplitude of the spit is shown in Figure 17. The amplitude starts to increase in 2000 and reaches a maximum in 2006, which is the same period as the northward sand spit growth. Figure 18 shows the variation in the spatial mean of the shoreline change. The spatial mean of the shoreline change shows the erosive nature starting in 1996 and the maximum erosion was found in 2008. There is a weak correlation (R = −0.53) between the amplitude of the spit and the spatial mean of shoreline change during 1972–2016. This may be because of the smaller amplitude of the spit which continued until 2000 (Figure 19). However, a higher correlation (R = −0.82) was found for the period of 2000–2016, when the sand spit was growing northward: the higher amplitude of the spit may have caused erosion to the north (Figure 19).

A summary of the recent morphological processes around the sand spit is illustrated schematically in Figure 20. According to image interpretation, the growing sand spit was accompanied by an adjacent erosion area observed in the northern delta. The amplitude of the sand spit seems to have a correlation with the amount of erosion downcoast. However, we do not know how the sand spit caused erosion. Generally, shoreline changes in a sandy beach are mainly because of longshore sediment transport. Thus, a numerical model was needed to evaluate the longshore sediment transport around the spit area as discussed in the following section.

5. Numerical Computation

5.1. Model Description

It is widely accepted that longshore sediment transport, or littoral drift, is one of the key factors in shoreline change for a sandy beach. On the other hand, cross-shore sediment transport is also vital because the coastal profile is formed by the erosion/deposition associated with shore-normal transport. However, cross-shore transport was not included in the present analysis as it is a seasonal process and can be neglected over the long-term. Therefore, the main goal of this part of the study was to estimate the longshore sediment transport numerically to explain the erosion/deposition pattern near the Reju River mouth area. To achieve this goal, the MIKE21FM Shoreline numerical model was used.

The basic concept (Figure 21) and a detailed description of the MIKE21FM Shoreline model is in Kristensen et al. [35,36] and Danish Hydraulic Institute [37]. This model consists of four models: a spectral wave model, a hydrodynamic model, a sand transport model, and a shoreline morphology model. The wave model computes the wave field considering linear wave refraction, linear wave shoaling, and wave breaking using the model of Battjes and Janssen [38]. The computed wave field from the wave model simulation is used as a forcing for the hydrodynamic model. The hydrodynamic model solves the depth-integrated Navier-Stokes equations. The computed flow field and the wave field are used as the input for the sediment transport model, which computes the sediment transport field. The computation of the sediment transport rates is determined using an intra-wave force balance description where the time evolution of the wave boundary layer is solved using the integrated momentum approach by Fredsoe [39]. The shoreline morphology model is based on a one-line formulation in terms of coordinates which follow the shape of the shoreline, instead of the more common approach where the two orthogonal horizontal directions are used. For long-term simulations, the current field, wave field, and sediment transport field from each module were coupled with the Shoreline Morphology Module to update the morphological bed level according to the sediment continuity equation. The bed level changes were continuously updated by the flow and wave fields until the simulation was complete.

5.2. Model Setup

The model domain covers the longshore extent from x = 49.4 km to x = 56.6 km, which is around the Reju River mouth area. The domain was gridded with an unstructured fine triangle volume mesh. We maintained a high-resolution mesh near the nearshore zone area, whereas a lower-resolution mesh was selected for further offshore. The bathymetric data used for this model originated from field survey data of the Institute of Water Modelling (IWM) from 2012. Offshore wave boundary data were collected from the ECMWF ERA-Interim archive. Sediment with a mean grain size d₅₀ = 0.1 mm, a grading coefficient = 1.5, and a porosity of n = 0.4 were used for the entire domain in all the presented simulations. The sediment data was provided by the IWM. Reju River input was not considered as there is no measurement of the flow rates of this river. We separately simulated the wave field and longshore sediment for (1) the road erosion period (2006–2008), and (2) after the revetment construction period (2008–2010). The initial shoreline position for every simulation was taken from the Landsat images, and the coastal profile of each simulation remained constant. We ignored cross-transport and tidal effects in the simulation. The time step used for simulations is 3600 s.

5.3. Model Results

Figure 22 shows the spatial distribution of the wave field (wave height and direction), wave ray, sediment transport field, and the distribution of cross-shore integrated alongshore sediment transport flux around a high-amplitude sand spit during a high energy event in 2007. Figure 23 illustrates the results for a low-amplitude sand spit during a high-energy event in 2009. The steep incident wave relative to the shoreline induced higher longshore sediment transport around the tip of the spit because the wave ray was concentrated around the tip of the spit in both cases. However, the longshore sediment transport was higher around the higher-amplitude spit than around the lower-amplitude spit because of the distribution of the alongshore sediment transport flux.

Since the littoral drift movement is parallel to the shoreline, there are two possible directions of motion, right or left, relative to the observer on the shore looking out to sea. The net longshore transport/littoral drift is the difference between the amount of littoral drift transported to the right and that to the left past a point on the shoreline. A numerical model was able to calculate the longshore sediment transport for every time step, and the results from the numerical simulation were used to compute the annual net littoral drift. The distribution of the net annual drift along the shore during 2006–2008 is shown in Figure 24. The positive net littoral drift indicates a predominantly northward longshore transport along the coast, which was expected because the dominant waves originate from the southwest direction. The net drift around the spit was larger compared to that of the straight shoreline up coast and down coast. Similar characteristics of net annual drift were observed during 2008–2010 as shown in Figure 25. However, the quantitative amount of annual drift was higher during 2006–2008 around the tip of the spit than during the later period.

According to one-line theory, the longshore sediment transport gradient is of interest because accretion is expected where the transport rate decreases alongshore, while erosion is expected where it increases. We attempted to verify the model result using the relation between the calculated gradient of net drift and the observed shoreline change rates. Figure 26 and Figure 27 show the relation between the gradient of the net drift and shoreline change rates during 2006–2008 and 2008–2010, respectively. Both results show a negative correlation between these two items which qualitatively supports the validity of the computation of the one-line theory. Moreover, we found a higher positive gradient of longshore sediment flux downcoast of the spit during 2006–2008 compared to that of 2008–2010, resulting in greater erosion downcoast during 2006–2008.

The net littoral drift from the up-drift coast and along the spit is mainly dependent on the dominant wave height and the orientation of the shoreline relative to the wave direction [15]. The sediment transport capacity becomes zero for waves approaching normal to the coast, while it reaches a maximum at a wave incidence of approximately 45°. We estimated the angle between the local incident wave angle to the local shore normal (α). The local incident wave angles were extracted at the closest grid point near the shoreline, and then averaged over the simulation period. Figure 24 shows the distribution of the mean incident wave angle ā with net littoral drift during 2006–2008. The variation in ā along the shore is no greater than 45°. We found higher transport around the tip of the spit compared to the straight shoreline up-coast and down-coast because of the high angle ā. A similar pattern of the mean incident wave angle ā was found during 2008–2010 as shown in Figure 25. However, the value of ā was greater during 2006–2008 around the tip of the spit than during the later period. Thus, the net annual drift around the tip of the spit was higher during 2006–2008 than the net annual drift during 2008–2010.

6. Discussion

The area of the delta changes in response to the position of the river mouth, rainfall, and waves. The position of the river mouth changes due to several breachings. This type of morphology is also found at the Kiawah River mouth, South Carolina [32]. The negative correlation between rainfall and the total delta area change indicates higher rainfall reduces the area of the delta. Thomas et al. [16] performed the correlation between rainfall patterns and Ginst spit growth, and they summarized that the higher rainfall results in a reduction in spit growth. Miller and Dean [34] established the correlation between deep-water wave steepness and shoreline changes in three different sites, and they found a negative correlation in most of the cases, which indicates the higher value of this wave property erodes the shoreline, and vice versa. In our study area, we also observed a similar relationship between wave steepness and the area of the total delta.

From the schematic diagram (Figure 20), we found a higher-amplitude spit (2006) induces greater erosion downcoast compared to a lower-amplitude spit (2008). According to the numerical results, the longshore sediment transport was greater around the spit in 2006 (a high-amplitude spit) compared to the spit in 2008 (a low-amplitude spit) because of the difference in the local incident wave angles (Figure 24 and Figure 25). Moreover, we found a higher positive gradient of longshore sediment flux downcoast of the high-amplitude spit resulting in greater erosion downcoast compared to that of a low-amplitude spit (Figure 26 and Figure 27).

Numerous studies have focused on the influence of sea level rise on the evaluation of different deltas around the globe, such as the Ganges [40], Mississippi [41], Song Hong [42], Nile [43], and Ebro [44]. Several studies [45,46] estimated sea level rise for Cox’s Bazar using 20-year tide data (1979–2000). The rate of sea level rise is about 1.4 mm/year. Smith [47] stated that global sea-level rise had increased from its 130-year average rate of 1.7 mm/year to about 3 mm/year over past 20 years. The rate of sea level rise for Cox’s Bazar is almost half compared to the global rate of sea level rise. Moreover, the impact of sea level rise could be devastating when land subsidence occurs. To our best knowledge, there is no alarming information about the land subsidence for this area.

In this study, the river sediment input is not measured, and it may have an important role in delta development. The smaller the delta, the stronger the influence of waves on overall delta evolution is likely to be [48]. The Reju River delta is relatively small and wave-induced longshore sediment may have a greater influence on this delta evolution than river sediment input. In addition, we used a single value for sediment size for the numerical simulation due to a lack of measurements. Thus, we assume the sediment distribution would not be drastically changed in the area.

7. Conclusions

In this study, we analyzed recent erosion downcoast of the Reju River delta on the Cox’s Bazar coast. Rapid erosion has affected portions of the 24-km-long Marine Drive Road along the coast, and local authorities have attempted to protect the road via a revetment. However, the structure was soon buried with sediment because of a growing sand spit along the northern delta. From image interpretation, we observed that a sand spit grew northward direction from 2000 to 2015, and the adjacent eroded region also moved in the same direction. The area of the delta changed in response to the position of the river mouth, rainfall, and waves. The amplitude of the sand spit seems to be correlated with the erosional amount downcoast.

The numerical results show that net littoral drift is dominant in the northward direction along the coast. The higher-amplitude spit induced higher sediment transport than that of the lower-amplitude spit because of the larger incident wave angle, and the higher positive gradient of the drift has caused erosion downcoast. Thus, a higher amplitude results in erosion downcoast.

Protection work near the delta is challenging because of the high variability in the shoreline position in this area. The revetment work in this area was not necessary; low-cost temporary protection would have sufficed. This study aids our understanding of the morphological behavior around the delta using limited data.

In this study, the seasonal change in the shoreline/bathymetry could not be verified because of the lack of data, which is among the key limitations. Moreover, the river sediment input is unknown and it may have an important role in delta area development. Long-term monitoring (field surveys) can aid in the accurate prediction of the future morphological behavior of this area that may benefit implementation of appropriate protection (hard/soft structures) measures.