Time-series analysis of landcover dynamics and their relation with coastline migration along Kuakata coast, Bangladesh using remote sensing techniques

ABSTRACT Time-series analysis of satellite imageries is useful in studying changes in coastlines and the nature of the landcover dynamics in coastal environments. This study aimed to investigate the relationship between coastal erosion-accretion and landcover changes in Kuakata. Landsat TM and Landsat OLI/TIRS satellite imageries at a nearly 5-year interval between the years 1989 and 2020 were used to compare changes within five major landcover classes in the study area – 1) mangrove vegetation, 2) settlements, 3) agricultural land, 4) waterbody and 5) beach. Net land loss over the past 31 years was estimated to be 1.73 km2 within the 93.05 km2 study area. Linear regression rates were calculated to identify the area most prone to erosion. The average erosion rate along the coastline was estimated to be 2.09 m/year. Most erosion occurred along the western part of the coast while the highest accretion was limited to an area in the east. Study findings also suggested that changes in beach and mangrove vegetation classes have a significant spatial and statistical correlation with coastal erosion-accretion processes. These findings can help the policymakers implement coastal zonation, and preventive and rehabilitative measures to save the tourism industry and agriculture in Kuakata.


Introduction
The erosion-accretion processes of the channels and islands of the Meghna estuarine delta make the coastline of Bangladesh highly dynamic.Most studies carried out on the changes in the coastline of Bangladesh indicate higher accretion than the erosion of the coasts ruled by estuarine conditions (Brammer, 2014;Islam et al., 2011).But this newly-formed land is less suitable for settlement and agriculture than older eroded land due to the raw state and alluvium salinity and is being continuously exposed to tidal and storm surge flooding.That's why the land use and landcover of the coastal regions are unstable as well (Mondal et al., 2018).In a recent study on the positional change of the Ganges deltaic shoreline for the period of 1977-2017, it was found that the shoreline is shifting landward at an average linear regression rate of 0.27 m/ year (Mullick et al., 2019).Examining the erosional processes can help the decision-makers, and administrators involved in coastal planning to understand the causes, develop forecasts, and plan how to implement preventive and rehabilitative measures (S. A. Hossain, Mondal, Thakur, Fadhil Al-Quraishi et al., 2022;S. A. Hossain, Mondal, Thakur, Linh et al., 2022).
Remote sensing techniques are extensively used to study changes in coastlines from time-series data (Mondal et al., 2021;White & El Asmar, 1999).Most studies suggest that the calculation of a single transform like the Normalized Difference Vegetation Index (NDVI) may be suitable for water delineation (M. A. Islam et al., 2014;Mondal et al., 2020).But the chaotic nature of the wave environment in the nearshore surf zone may result in high variability in the reflectance values in the red and infrared bands which is why they alone are unable to adequately separate the beach from the surf zone (Daniels, 2012).A hybrid algorithm for coastline delineation is likely to be preferable (Bushra, 2013;Daniels, 2012;Emran et al., 2017;Mullick et al., 2019).With the application of GIS, the accuracy of coastline delineation can be enhanced.Digitization is the most popular coastline delineation option (Rahman et al., 2013).Using hightide coastlines is deemed to be the most convenient technique in case of unavailability of time and/or data on bathymetry and beach slope (A.Z. Z. Mondal et al., 2020;A. Z. Z. Islam et al., 2013).On the other hand, analyzing landcover changes from medium resolution images of coastal areas using supervised classifiers like SVM, ANN, fuzzy analysis and segmentation comprises complexity and bias (Yuan et al., 2009).This results in lower overall accuracy compared to unsupervised classification schemes in classifying homogeneous (50% pixels classified as one landcover class) study areas (Kaliraj et al., 2017;Peacock, 2014).
Kuakata, located at the southernmost tip of the coastline of Bangladesh is a panoramic beach that offers a full view of sunset and sunrise.In 2016, the total gross recreational benefit of this area was estimated to be approximately 29.55 million per year in Bangladeshi Taka (M. S. Hossain & Islam, 2016).Geographic demarcation of the area is 90° 4' 12.0648 " E, 21° 47' 32.3628" N and 90° 16' 52.6296" E, 21° 53' 11.0256" N (Figure 1).This 93.05 km 2 coastal area at the north of the Bay of Bengal belongs to the Lata Chapli and Dhulasar unions of Kalapara Upazila in Patuakhali district within the Barishal division.Being part of the greater Ganges-Brahmaputra-Meghna (GBM) basin and Meghna Estuary, this area is susceptible to the estuarine processes of the river Galachipa and tidal actions of the Bay.It has a tropical climate with an average annual temperature of 26°C and about 2748 mm of precipitation occurs throughout the year (Climate-data.org,2022).Its climate makes Kuakata prone to tropical cyclones (i.e., cyclone Sidr in 2007, cyclone Aila in 2009, cyclone Mahasen in 2013 etc.) and storm surges.Apart from that, coastal erosion and salinity intrusion are severe problems in the area, threatening tourism and agriculture there.The high and low-velocity longshore currents combined with the convex-shaped coastline makes Kuakata more susceptible to marine processes reflecting seasonal variations, and cause coastal erosion in some places and accretion in others (Bushra, 2013;Rashid & Mahmud, 2011).Previous studies on this coast suggest that more erosional activities occur along the coastline compared to coastal accretion  (Bushra, 2013;A. Z. Z. A. Z. Z. Islam et al., 2013;Rahman et al., 2013).However, these studies on the Kuakata coast neither focused on landcover change analysis nor tried concluding the interrelationship between coastline shifting and subsequent landcover changes.But to mitigate the effects of coastal erosion and save the tourism in Kuakata, changes in landcover associated with erosion and accretion need to be explored.The present study aimed to determine whether any distinctive pattern exists relating coastal erosion and landcover dynamics of Kuakata.

Data set and pre-processing
For this study, satellite imageries (Landsat TM and Landsat OLI/TIRS) from 1989 to 2020 were used (Table 1).These images are from the same season of the respective years and were acquired around the same time of the day.All Landsat image bands have 30 m spatial resolution (bands used are listed in Appendix A).
The Landsat images were downloaded in GeoTIFF format (USGS, 2019).These Landsat images had been radiometrically calibrated and orthorectified using ground control points (GCPs) and digital elevation model (DEM) data to correct for relief displacement (USGS, 2019).

Landcover classification
Experimentation was done using both supervised and unsupervised pixel-based classification algorithms but an unsupervised clustering algorithm was chosen for this research over a supervised algorithm.Because these clustering algorithms seek to find classes that are most different while supervised algorithms based on training data were causing a higher level of lumping between classes.Among various clustering algorithms, ISODATA was used in this study.The images were initially classified into 100 landcover classes from which they were recoded to 7 broad classes -1) mangrove vegetation, 2) settlements, 3) dry agricultural land, 4) wet agricultural land, 5) inland waterbody, 6) sea, and 7) beach based on ground data and visual interpretation.Dry and wet agricultural lands were later merged as "agricultural land" while inland and sea water were merged to become a single class called "waterbody" to reduce complexity in change detection.

Validation of landcover classifications and post-classification comparison
A four-day field campaign was carried out in the study area focusing especially along the coastline to collect ground data on high waterlines and landcover types, as well as GPS coordinates using a handheld device.One hundred and forty-two ground data points were collected which were then used for calculating four accuracy metrics to evaluate classification accuracyoverall accuracy, kappa coefficient, user's accuracy, and producer's accuracy.
Change detection was performed between successive pairs of classified Landsat images to show the sequence of change and also between the first date and last date to show the resultant changes over the study years.

Coastline extraction and erosion rate calculation
For coastline delineation, a hybrid algorithm using NDVI (Normalized Difference Vegetation Index) and Tasseled Cap Analysis followed by ESRI's unsupervised ISO classification was used in this study (Daniels, 2012;Mondal et al., 2020).In the case of NDVI, vegetative "spectral signatures" are calculated based on the principle that "actively growing green plants strongly absorb radiation in the visible region of the spectrum and strongly reflect radiation in the NIR region" (Sultana et al., 2014).
Tasseled Cap Analysis is based on a linear transformation of data from the original image into three new axes (brightness, greenness, and wetness) which become features of the transformation and allows the creation of other useful raster by-product files like polylines (coastlines; Scott et al., 2003;Vorovencii, 2007).These three bands and the NDVI band were used for ESRI's unsupervised ISO classification algorithm to obtain a binary classification of land and sea for all seven Landsat images (Daniels, 2012).Then the coastlines were extracted using GIS functions -majority filtering, contour, smooth line, and finally, manual editing to correct for overlapping of other features.Additionally, the 7-class unsupervised ISODATA classified Landsat images were reclassified to land and sea to extract erosional and accretional polygons of each study year.These polygons were used for calculating total coastal erosion and accretion in 5-year windows (1989-1994, 1994-2001, 2001-2006, 2006-2010, 2010-2015, 2015-2020, and 1989-2020), as well as between the study years of 1989-2020.Erosion rates for the extracted coastlines were calculated automatically in Digital Shoreline Analysis System (DSAS) version 5.0 at 30-m interval transects (Himmelstoss et al., 2018).As collecting data for the total uncertainty term is time-consuming and expensive, it was ignored.Instead, the coastlines were collected from the same season to reduce the magnitude of uncertainty (Ruggiero & List, 2009).

Landcover classification, results validation and post-classification comparison
The total area of each landcover class for the years 1989, 1994, 2001, 2006, 2010, 2015, and 2020 in km 2 (Figure 3) was estimated from the unsupervised ISODATA classification of Landsat images (Figure 2).It is seen that even though the share of agricultural lands had been reducing over the years, it is still the most dominant landcover class in the study area.Despite the afforestation attempts of forest division, mangrove vegetation has been reduced to some extent.However, the area coverage of this class has increased in 2020 compared to 2015.The most increase in area is seen in the settlements class.
The total beach area has been fluctuating over the years due to erosional and accretional activities.The highest accretion occurred between 1994 and 2010 after which the beach area started decreasing, and reaching the lowest value of this class in 2020.It indicates that the overall erosion rates were highest from 2010 to 2020.The drastic changes in landcover classes between 1989 and 2020 are shown in Figure 4.
Validation of image classification was done by calculating accuracy metrics (Table 2) from the confusion matrix (Appendix C) generated using the ground truth data against the pixel data from the Landsat OLI image of 2020.
Change detection results (Appendix B) suggest that mangrove vegetation was mostly changed to water class (Figure 4).We can assume this occurs due to wave actions causing coastal erosion.Over the years, 1.51 km 2 of mangrove forests were consumed by the sea.From 1989 to 2020, 0.72 km 2 of mangrove vegetation was changed to agricultural land, and the 0.59 km 2 changed to beach were most likely to occur due to the erosion of mangrove areas and subsequent coastal accretion.In the case of settlements, it is seen that between 1989 and 2020 about 1.22 km 2 were converted to agricultural land.It indicates that the settlement vegetation was cleared for agricultural purposes.But the highest change of 3.04 km 2 occurred between 2010 and 2015.0.67 km 2 of settlements were changed to inland water while 0.15 km 2 of settlements were changed to beach features during the study years.This can be explained as the combined result of coastal erosion-accretion processes and forestation attempts.The same can be said for the changes in agricultural lands to beach and mangrove class.In fact, the most dramatic change has been observed in agricultural use.From 1989 to 2020, 11.68 km 2 of agricultural land was changed to settlements, and 4.95 km 2 were lost to waterbodies (Figure 5).1.84 km 2 of beach area were lost to the sea and 1.21 km 2 were changed for agricultural purposes during the study years (Figure 6).
Change detection results also suggest that mostly waterbodies were changed to agricultural lands, the highest being 3.50 km 2 between 2015 and 2020 (Figure 6).This indicates sediment deposition was highest during this period that allowed the newly formed land to be used for agriculture.

Coastline extraction and erosion rate calculation
In this study, linear regression rates (LRR) were chosen to present the computational results of DSAS.LRR is determined by fitting a least-squares regression line to all coastline points for a transect so that the sum of the squared residuals is minimized.LRR is the slope of the line whose positive and negative values correspond to coastal accretion and erosion (Figure 7(b) and Table 3), respectively (Sutikno, 2016).
It can be seen that erosion rates are highest along the western coast facing the bay while accretion rates are higher at estuary mouths on both sides of the coastline.
Total eroded and accreted areas extracted from unsupervised ISODATA classification of 1989 and 2020 images (Figure 7a) seem to correspond to the LRR calculated by DSAS.Areas having negative LRR are found to be eroded while areas having positive LRR are attributed to accretional areas.Based on these results, it can be said that most erosion occurred along the western part of the coast while accretion was highest in a limited area in the eastern part.From 1989 to 2020, 4.76 km 2 of coastal land was eroded while 3.03 km 2 was accreted in Kuakata.So, the net land loss over the past 31 years is 1.73 km 2 (Figure 8).

Relationship between coastal erosion and landcover dynamics
The graphical representation of the relationship between changes in major landcover classes and coastline movement in the study area is shown in Figure 9.
It can be realized that erosion and accretion are a continuous process along the coastline.In the case of accretion, the trend was alternating increase and decrease every 5 years.However, the trend in erosion was not similar.From 1994 to 2010, the total eroded area per 5 years was on the decline and the lowest value being 0.79 km 2 between 2006 and 2010.In fact, 2006-2010 was the only period when the total eroded area was a little lower than the total accreted area.It means during this period some overall gain may have occurred along the coast.But after 2010, the total eroded area had been increasing per 5 years that    means continuous erosion had been occurring along the coast in the past decade.In the case of landcover changes, mangrove vegetation was shown to be on the decline until 2015.Such a scenario was also evidenced by field observation (Figure 10).Agricultural land share showed a decreasing trend, the exception being during 2010-2015 when a significant increase was noticed (Figure 9).That is also the period when the highest coastal accretion occurred.On the other hand, the settlements class shows an increasing trend in area coverage that is only natural due to population rise.But an exception was noticed between 2010 and 2015 when it showed a decrease.The loss due to erosion will not be reflected in the statistics of this class as the change will be shared by other classes.To truly understand the change, we have to take a look at Figure 11 which shows the settlement areas that were lost to coastal erosion.
It can be seen that major classes -mangrove forests, beaches and agricultural lands were lost to erosion and these classes are now changed to seawater.On the western part of the coast, roads, embankments, and beachside settlements were lost due to erosion.Beach class, however, shows a very high correlation with coastal erosion (Figure 9).From 1994 to 2010, erosion per 5 years was on the decline while the percent change of beach area was on the rise.Again during 2010-2020, the beach area had been declining while the total eroded area has increased per 5-year interval.

Discussion
This study aims to explore the relationship between coastal erosion and landcover changes in Kuakata.It was evident that as long as erosion continues, so will the mangrove deforestation.If afforestation efforts succeed then there may be an overall gain in mangrove forest share which is the case occurring between 2015 and 2020 (Figure 9).According to the Gongamati Reserved Forest Range officer, afforestation efforts on accretional lands of the Dhulasar coast started in 2007.The second phase of the plantation project was carried out in 2011.These efforts resulted in an increase of 0.85 km 2 in mangrove vegetation between the study years.The 1.72 km 2 of waterbodies changed to settlements (Appendix B) attributes to river sedimentation allowing settlement build-up.The decrease in settlement area between 2010 and 2015 (Figure 9) could be due to the impact of Cyclone Mahasen in 2013 which had devastating impacts on infrastructures and settlements in the Patuakhali district.As the population grows, the settlement area is bound to expand.When one region is flooded or eroded, families relocate landwards and compensate for the loss of housing and infrastructure by clearing forests or using up agricultural lands for settlement purposes.
1.28 km 2 of the beach area was accreted between 1989 and 2020 (Appendix B).But accretion rates are much lower than erosion rates as sediment deposition and settlement is a slow process.Even though there is no net land gain in the study area, coastal accretion does occur.Initially, the newly accreted lands are shown to be beach feature class.Later on, they are either changed to mangrove vegetation or agricultural land.Similar study can be conducted on the entire coast of Bangladesh since previous studies on the coastline of Bangladesh focused only on the erosion rates, not the implications in landcover changes.In addition, considering total uncertainty terms should provide more accurate erosion rates.

Conclusion
Kuakata beach, an important tourist destination in Bangladesh is subjected to continuous coastal erosionaccretion processes over the years and currently, this scenic landscape is on the brink of disappearance.Interaction between two strong opposing factors -    fluvial delta-building activities and marine processes results in this dynamism.Along the coast, accretion and erosion are taking place at different rates and the natural environments are being changed in various ways by human interventions.The findings from this study suggest that beach and mangrove vegetation are most affected by coastal erosion-accretion processes and changes in these classes have a significant correlation with coastal processes.In order to save the tourism of this economically and aesthetically significant area, appropriate beach protection measures should be taken in the areas shown to be most prone to coastal erosion.Newly accreted areas can be considered for relocation of the tourism industry and future developments in Kuakata.

Figure 1 .
Figure 1.Location map of the study area.

Figure 2 .
Figure 2. Landcover classification results of Landsat images from 1989 to 2020.

Figure 3 .
Figure 3.Total area (km 2 ) occupied by landcover types from the year 1989 to 2020.

Figure 8 .
Figure 8. Eroded and accreted areas in Kuakata at 5-year intervals from 1989 to 2020.

Figure 9 .
Figure 9. Landcover changes shown as percentages vs. coastline movement in km 2 .Negative values indicate the percentage by which the share of each landcover class has been decreased during the 5-year interval while positive values indicate an increase.The total area of coastal erosion and accretion during the time interval are shown on the secondary axis.

Figure 10 .
Figure 10.(a) and (b) denuded forests on the coast of Dhulasar; (c) denuded mangrove vegetation on the western coast of Lata Chapli; (d) seawater encroachment into forests near Dhulasar coast.

Figure 11 .
Figure 11.Landcover changes specific to the accretional and erosional zones between 1989 to 2020.

Table 1 .
Overview of data used.

Table 2 .
Image classification accuracy assessment results.

Table 3 .
Summary statistics of erosion rate calculation.