Shoreline Change Analysis of Hooghly Estuary Using Multi-Temporal Landsat Data and Digital Shoreline Analysis System

Abstract

Long-term (1973–2021) shoreline displacement, rate of change, and temporal pattern were examined using multi-date Landsat data and Digital Shoreline Analysis System (DSAS) along the 200 km coast of Hooghly estuary. Orthogonal transects of 100 m apart were casted for calculation of End Point Rate and Weighted Linear Regression rate on different temporal scales for seven analysis zones. The shoreline change pattern was established using Hierarchical agglomerative clustering. The study reveals that almost 43.45% of the beachfront has eroded and 56.55% has accreted during the past four decades. The average erosion rate varies between −0.01 and −13.71 m yr.−1 and accretion of −0.01 to −22.30 m yr.−1. The littoral drift resulted in a maximum seaward aggression by 1096.89 m in the zone 1. Landward movement was maximum (−602.96 m) in the zone 4. Although west bank is prograding @ 3.47 m yr.−1 (±5.83), the east bank is eroding @ 1.30 m yr.−1 (±4.08). Based on the cluster analysis about –1.87% of the shoreline exhibits consistent erosion over all the intervals, whereas trend was evident in 4.73% of the coastline. The portions of coastlines, which exhibit high erosion rate and consistent erosion need immediate attention and policy intervention.

1. Introduction

The shoreline is the physical interface or intertidal margin between land and sea and constitutes one of the 27 global “Geo-indicators” referred by the International Union of Geological Science [1] and International Geographic Data Committee (IGDC). Shoreline change is a dynamic natural process in the coastal areas induced by erosion/accretion that occurs over a range of temporal scales. The morphological evolution of the Hooghly estuary and its coastline is the result of two counteracting transport processes of sediment supply versus removal. When both the processes are balanced an equilibrium is reached. However, most often this balance is disturbed due to the influence of episodic and/or long-term natural forcing and anthropogenic interventions. As a consequence, the shoreline keeps changing its position [2, 3, 4, 5, 6] over a wide temporal scale, from geologic age to short-lived, extreme weather events such as storms and tsunamis. The long-term processes that shape the shoreline include sea-level rise (SLR), altered wind patterns [7], frequency and intensity of storms [7], offshore bathymetric changes [8], high energy swells [9] and supply of fluvial sediment input. In addition, anthropogenic activities viz., landcover changes in the river catchment, port and harbor and dam construction, dredging for maintaining navigation channels, aquaculture, protective embankments, beach nourishment, economic and tourist activities also exacerbate the coastline change on a short temporal scale. Engineering structures change the estuarine circulation patterns and may change the freshwater flow along with sediment and nutrient supply. In several instances, engineering modifications to the beach creates discontinuities in the historical shoreline position and mask underlying long-term behavior [10]. Another less reported phenomena are land subsidence which may occur naturally due to compaction of sediments or triggered by the excessive withdrawals of ground water. In general, the coastal landform establishes a morphodynamic equilibrium after episodic short-term perturbations. However, many times the combination of natural and manmade activities exacerbates the shoreline change and exhibits non-linear morphological responses to change [11].

According to Williams [12], the study of shoreline variation and forecast plays an important role in coastal zone management and it becomes more crucial in the context of anticipated climate change and sea-level rise [13]. In this context, one of the key requirements for effective coastal zone management is the availability of accurate position of the shorelines for analysis of changes in the past and future trends. Traditional methods of shoreline delineation include terrestrial surveys using landmarks, aerial photos [14, 15], Global Positioning Systems (GPS), terrestrial Light Detection and Ranging (LiDAR) or 3D scanners. But they are time-consuming, labour intensive and costly. In contrast the remote sensing data form space platform is more convenient, easy to process and above all freely available in the public domain. Remote sensing data has been extensively used in shoreline change studies because of their synoptic and repetitive coverage, multispectral capabilities enabling contrast between land and water in the infrared portion, and cost-effectiveness [14, 16]. Advanced image processing techniques can be employed on satellite data for precise extraction of the shoreline. Some of the methods used by different researchers include threshold level slicing and image classification technique [17], density slicing of TM band 5 [18], canny edge detection using DN threshold [19], mean shift segmentation [20], pixel-based segmentation using DN threshold [21], neural network [22], fuzzy logic [23, 24], texture analysis [25], machine learning [26] and incorporation of ancillary spatial data in the classification scheme [27, 28, 29]. Quantitative assessment of the spatio-temporal variation of shoreline at global scale has been carried out by several authors [30, 31, 32]. In this endeavor the twin technologies of Remote Sensing and Geographic Information System has been recognized as the most useful tools for quantifying the historic shoreline change [33, 34]. To avoid the discrepancy which might be introduced due to fluctuation of water level Yu et al. [35] have used satellite images obtained at similar tidal heights. Chen and Chang [36] have done the tidal correction using high spatial resolution satellite images and real-time data of tidal level to reduce the impact of tidal level variability on the estimation of coastline change. In India also several studies have been carried out for shoreline change analysis using remote sensing data [37, 38, 39]. Most of the studies have used Digital Shoreline Analysis System [40], a software extension within the ArcGIS tool for measuring, quantifying, calculating and estimating of rate of change from multiple historic shoreline positions at different temporal scales [41, 42, 43, 44]. The change metrics of DSAS are Net Shoreline Movement (NSM), Shoreline Change Envelope (SCE), End Point Rate (EPR), Linear Regression Rate (LRR) and Weighted Linear Regression Rate (WLR) among others. LRR and WLR enable multiple historic shorelines to be used to determine the rate of change by fitting a least-square regression line to all shoreline points for particular transects.

In the present study, Landsat satellite data of 8 temporal intervals between 1973 and 2021 were used for land-water discrimination, generation of shorelines and long-term change rate along with change pattern along the Hooghly estuary. The instantaneous land-water boundary was used as coastline which is relatively simple and can easily be identified using image transformation. The main objectives of the study are i) medium- and long-term changes in the shoreline at high spatial resolution using DSAS ii) to identify the erosion/accretion pattern and iii) to examine the role of change drivers.

The findings of the study will be useful for the managers and engineers to make scientific and rational policies for land use planning, to develop effective coastal protection strategies, predicting capacity for future coastal change due to climate and other drivers and improving impact and vulnerability assessments that include natural human sub-system interactions.

2. Study area

The Hooghly estuary is located in the southernmost part of Indo-Gangetic plain, flanked between East Midnapur (in the West) and South 24 Parganas district (in the East), extending between 21^o33′10′′N to 22^o13′16′′N latitude and 87^o45′00′′ to 88^o18′22′′E longitude (Figure 1). The head Bay is a unique deltaic environment comprising a wide continental shelf, complex coastal geometry and high tidal range. Tide domination is indicated by exponentially tapering channels, with funnel-shaped mouths [45]. The region has formed, sculptured and modified due to continuous fluvial action of the Ganga and the Brahmaputra systems, intense tidal hydrodynamic behavior, climatic disturbances and anthropogenic activities [46]. The funnel-shaped estuary has a width of 6 km at its head and 25 km at the mouth, responsible for tidal asymmetry and flow variation leading to bank erosion [47]. The average depth of the water column is only 6 m [48]. The estuary receives 4 tributaries viz. Damodar and Rupnarayan rivers at its head, and Haldi and Rasulpur rivers at the middle at its west bank. In contrast, the east coast is punctuated by several closely spaced inlets. Based upon the tidal amplitude, the coastal region of West Bengal can be sub-divided into i) macro tidal (tidal range > 4 m) from Sagar to Bangladesh border and ii) meso-tidal (tidal range of 2–4 m) mostly Medinipur (Digha-Sankarpur-Junut) coastal plains to the west of the Hooghly estuary.

Figure 1.

Index map of the Hooghly estuary.

Geologically the basement of the Bengal basin is a part of the eastern edge of the Indian plate, which is being subducted beneath the China plate along the Sunda subduction zone and Naga-Lushai orogenic belt. The tectonic and depositional history of the Bengal basin has been controlled by several movements during Cretaceous-Tertiary periods. Due to the tectonic activity the Bengal basin has been tilted towards east resulting in successive changes in the course of the Ganga River towards east from the historical past. Due to this shifting, the deltaic region suffers from the paucity of fresh water discharge and sedimentation. Auto compaction of loosely attached sediments and gradual land subsidence is also another prominent geomorphic event occurring in this region [49, 50, 51, 52, 53, 54] which mostly remains unnoticed. Morphometrically the Hooghly estuary is the product of continuous fluvial sedimentation in a series of para-deltaic lobe progradation systems developed on the western shelf margin areas and eastern troughs of the Bengal basin caused by the eustatic, isostatic and tectonic forces. The coastline presents various landforms such as tidal/mud flats, sandy beaches (located near Digha, Duttapur, Shyampur, Dadanpatra, Baguranjalpai, Dariapur and Nij Kasba), salt marshes (near Khejuri and at the mouth of Rasulpur river near Nij Kasba) and mangrove marsh (south of Patibunia). A vast extension of the muddy beach is found in South 24 Parganas, especially to the east of Bakkhali. The most striking feature is the development of successive rows of dunes (both Palaeo and Neo dunes) with intervening clayey tidal flats in the south of East Midnapur district between the stretches of Subarnarekha and Hooghly estuary is due to punctuations in the regression of the sea during Holocene [55]. Banerjee and Sen [56] opined that the regression of sea along this coastal tract is around 6000-year BP which resulted in seaward shifting of shoreline and formation of Paleo-dunes. Accordingly, to Niyogi [57], six regular cycles of beach ridges alternating with a variable number of bars are visible in the area, which is indicative of the shifting of shorelines. According to Gaur and Vora [58], the shoreline position was 5–15 km inland from the present shoreline around 6000-year BP. The erosion and accretion patterns clearly show a continued geomorphic sculpturing of the Hooghly coast.

To capture the micro-level variability, alongshore is divided into 7 analysis zones (Figure 2) covering both the west and east bank. The zones in the west bank are delimited by the main inlets which are the freshwater sources, eventually draining into Bay of Bengal. The area delimitation of various zones, constituting transects and shoreline distances is given in Table 1. The west bank is divided into 3 zones whereas the east bank into 4 zones (Table 1). The total length of the coastline studied is 200 km of which 90 km on the western side and 110 km on the eastern side of the estuary. The studied coastline was divided into 1924 number of transects (T_n) separated by 100 m. The number of transects increases from west to east bank in the clockwise direction.

Figure 2.

Different analysis zones.

Zone	No of transects; (from-to)	Location	Distance (km)
West bank
Zone–1	187 (T25-T211)	Pichhabani outlet to Rashulpur river	19.37
Zone–2	307 (T218-T524)	Rashulpur river to Haldi river outlet	31.04
Zone–3	384 (T534-T917)	Haldi river to the confluence of Rupnarayan and Hooghly River	39.52
East bank
Zone–4	256 (T919-T1174)	Confluence of Rupnarayan and Hooghly river to Kulpi	26.23
Zone–5	280 (T1175-T1454)	Kulpi to Kakdwip	27.75
Zone–6	120 (T1455-T1574)	Kakwip to Namkhana	13.41
Zone–7	390 (T1575-T1967)	Namkhana to Henry Island	43.10

Table 1.

Salient description of different analysis zones.

3. Materials and methodology

3.1 Data used

3.1.1 Army and survey of India Topomaps

The historic shorelines were digitized from Army Map Series (NSS&H, Edition-1, AMS) in 1:250,000 scale surveyed during 1942–1943) number NF-45: 7 (north of study area) and 11 (south of the study area) were used for the coastline change analysis. Besides Survey of India topomaps of 73 N-16, 73O -13,14; 79B - 4; 79C-1,2,6 surveyed during 1967 were also used for generation of high-water level (HWL) coastlines.

3.1.2 Satellite data

Landsat satellite data of 1973 to 2021 have been used for decadal and long-term trend analyses. The data has been selected based on clear sky condition, high tide date and time as well as season. For discrimination of land-water boundary shortwave infrared bands 5 (1.55–1.75 μm) and 7 (2.08–2.35 μm) of Landsat - 4, 5, 7 and bands 6 (1.566–1.651 μm) and 7 (2.107–2.294 μm) of Landsat – 8 (OLI) were used. The details of the satellite data used in the study are given in Table 2.

Satellite/sensor	Path/Row	Date of overpass	Spatial resolution; (m)	Overpass time; (local time)	Time of high tide; (local time)	Tide height; (m)
MSS1	149/45	17.01.73	60	NA	09:40	4.05
MSS3	149/45	17.01.80	60	03:52	10:33	4.59
TM4	138/45	19.01.89	30	04:03	9:08	3.45
TM5	138/45	28.01.95	30	03:43	8:54	3.85
ETM + 7	138/45	06.03.00	30	04:23	11:15	4:73
TM5	138/45	07.01.05	30	04:17	8:14	3.81
ETM + 7	138/45	29.01.10	30	04:22	10:17	4.38
OLI 8	138/45	09.02.17	30	04:31	9:48	4.04
OLI 8	138/45	24.03.21	30	04:30	7:42	3.19

Table 2.

Details of the satellite data used.

The tide information is pertaining to the Diamond Harbor station.

For the of satellite data tide and current prediction programme viz. WXTide32 was used which allows knowing the time of high-tide and low-tide, as well as the tide height (m). It is supported by more than 9500 stations worldwide with the capability to predict tides from 1970 through 2037. One of the nearest stations of the study area, e.g., Diamond Harbor was selected to know the high tide date, time and magnitude. The high tide timings were compared with the satellite overpass dates and time for the selection of images representative of coastlines on high-tide date. This enabled comparison of the shorelines under identical tide conditions by minimizing the variability due to the tidal cycle.

3.2 Methodology

3.2.1 Land-water discrimination and shoreline extraction

There are seven types of coastline indicators viz. geomorphological reference lines, vegetation limits, instant tidal levels and wetting limits, tidal data, beach contours and storm lines. In the present study, the high-water levels (HWL) during spring tide have been considered for historical coastline change analysis. Several water indices have been used by many authors to extract the coastlines [59, 60]. Two of the most popular water indices are the Normalized Difference Water Index (NDWI) and Modified Normalized Difference Water Index (MNDWI). NDWI was primarily developed for Landsat MSS whereas MNDWI was developed for TM, ETM+ and OLI sensors [61]. The bands used to generate these indices mostly consist of green, near-infrared, middle-infra-red and shortwave infrared bands. The objective of all these indices is to enhance the contrast between land-water interfaces. Band 2 (0.52–0.60 μm) and band 5 (1.55–1.75 μm) of Landsat-7 and Band 3 (053–0.59 μm) and band 6 (1.57–1.65 μm) of Landsat-8 OLI were used for computing NDWI. The formulae of NDWI and MNDWI are given below:

$${\text{NDWI}}_{\mathrm{MSS}}=\left({\mathrm{R}}_{\text{green}}-{\mathrm{R}}_{\mathrm{nir}}\right)/\left({\mathrm{R}}_{\text{green}}+{\mathrm{R}}_{\mathrm{nir}}\right)=\left({\mathrm{B}}_1-{\mathrm{B}}_4\right)/\left({\mathrm{B}}_1+{\mathrm{B}}_4\right)$$E1

$$\mathrm{M}{\text{NDWI}}_{\mathrm{TM},\mathrm{ETM}+}=\left({\mathrm{R}}_{\text{green}}-{\mathrm{R}}_{\text{swir}}\right)/\left({\mathrm{R}}_{\text{green}}+{\mathrm{R}}_{\text{swir}}\right)=\left({\mathrm{B}}_2-{\mathrm{B}}_5\right)/\left({\mathrm{B}}_2+{\mathrm{B}}_5\right)$$E2

$$\mathrm{M}{\text{NDWI}}_{\mathrm{OLI}}=\left({\mathrm{R}}_{\text{green}}-{\mathrm{R}}_{\text{swir}}\right)/\left({\mathrm{R}}_{\text{green}}+{\mathrm{R}}_{\text{swir}}\right)=\left({\mathrm{B}}_3-{\mathrm{B}}_6\right)/\left({\mathrm{B}}_3+{\mathrm{B}}_6\right)$$E3

Where, R_green = spectral reflectance of the green band, R_nir = spectral reflectance of near-infrared band and R_swir = spectral reflectance of the shortwave infrared band.

Before applying the water index on Landsat MSS data of 1973, the image was resampled to 30 m spatial resolution to make the resolution comparable with the rest of the datasets. A Boolean approach was used on the NDWI/MNDWI images to create two classes viz. land and water. The threshold for land-water boundary was kept >0.15 for OLI and > 0.1 for TM/ETM+. The resulting image had only two classes viz. land and water. However, some of the inter-tidal zones could not be demarcated due to very high sediment loading. Hence, a hybrid approach was followed wherein NDWI/MNDWI, SWIR and topomaps were used for precise demarcation of the shorelines. Vectorization of the land-water boundary was done using the region growing tool of ERDAS/Imagine (ver. 9.1). The spectral Euclidean distance was set interactively to accurately capture the land water boundary. Some manual editing was also done on the shoreline vector.

3.2.2 Shoreline change analyses

Historical shoreline behavior was examined using Digital Shoreline Analysis System (DSAS, ver. 5.0), an extension tool of ArcGIS software (developed by the US Geological Survey) which calculates several change statistics viz. Net Shoreline Movement (NSM), Shoreline Change Envelope (SCE), End Point Rate (EPR), Linear Regression Rate (LRR) and Weighted Linear Regression Rate (WLR). It can analyze the time series of multiple shoreline positions [34] by using linear regression fit. Provision is there to include uncertainty of the input data in terms of assigned weights. Functionally DSAS performs 3 major activities viz., i) defining a baseline ii) generation of orthogonal transects and iii) computation of rates of changes. The tool enables the calculation of scales and rates of change statistics from multiple historic shoreline positions and sources. DSAS is a freely downloadable tool and is available at Woodshole [62] http://woodshole.er.usgs.gov/project-pages/dsas/. The details are available in Thieler and Danforth [63, 64] and Thieler et al. [34]. A brief description of the change statistics used in the study is given below.

Shoreline Change Envelope (SCE): A measure of the overall change in shoreline (m) at each transect considering the farthest and nearest position of the shoreline for the baseline location [65] irrespective of the dates.

Net Shoreline Movement (NSM): It is the distance (m) between the oldest and the youngest shorelines [66].

End Point Rate (EPR): It is derived by dividing the distance of shoreline change between two time periods by the time interval and expressed as m yr.⁻¹ [65, 66, 67, 68, 69]. This method provides the net rate of change over the long term. It has both advantages and disadvantages; the advantage is that only two shorelines are required for computation of change rate but unable to use more than two date shoreline data.

Linear Regression Rate (LRR): It determines a rate-of-change statistic by fitting a least square regression line using all the intersection points all shorelines and individual transect [66, 68, 70]. The slope of the line is the rate of shoreline change. The advantages of linear regression include i) all the time-series data are used and ii) can reduce the impact of spurious values on the overall accuracy of change rate [71]. To ensure meaningful results from the regression model the temporal intervals of shorelines were kept well distributed over the analysis period.

Weighted Linear Regression Rate (WLR): This method takes into account the positional and measurement uncertainty of the shoreline positions [63, 72, 73, 74] especially when the shorelines were digitized from various sources and scales. The uncertainty values are incorporated in the DSAS as weights. The slope of this regression line is the shoreline change rate (m yr.⁻¹). Using the weighing factor, WLR, standard error of the estimate (WSE), standard error of slope with user selected confidence interval (WCI) and R-squared value (WR2) are obtained [34]. The results of this method are controlled by the points with smaller positional uncertainty on the best-fit regression line [65]. If no values are provided by the user, DSAS uses the default uncertainty value. The weight (w) is defined as a function of the variance in the uncertainty of the measurement (e) [34]:

$$\mathrm{w}=1/\left({\mathrm{e}}^2\right)$$E4

Where e = shoreline uncertainty value.

Coefficient of Determination (R²): It is the percentage of variance in the data that is explained by a regression [65]. It is used to verify the quality of the best-fit line regression.

3.2.3 Calculation of data uncertainty

The errors or uncertainties that arise due to different data sources, time of data acquisition, and the type of shoreline indicator were quantified based on several studies [73, 75]. According to Fletcher et al. [75] and Romine and Fletcher [76] there are two types of uncertainty: positional (seasonal and tidal fluctuations) and measurement (digitizing, pixel and rectification error). The uncertainty for each dataset was worked out considering the data product with due weightage of the quality of each data. The total uncertainty is used to calculate the weight and further working in the DSAS. Different uncertainties are explained below.

Seasonal error (Es): It is the error that arises due to seasonal changes of shoreline positions under the action of the waves and storms [75]. In the present study, all the scenes are of the winter season and hence this error was neglected.

Tidal fluctuation error (Et): It is the error associated with horizontal variability in shoreline position due to tides [75]. All the images in the present study correspond to high-tide values. Based upon the values the tidal range was considered as 2.12 m.

Digitizing error (Ed): It is the error related to shoreline digitization [75]. The digitizing error was kept within half a pixel of Landsat data (15 m).

Pixel error (Ep): It relates to image precision (resolution). In the present study except for the Landsat MSS data all the TM, ETM+ and OLI data have a spatial resolution of 30 m. To make the pixel uniform the 60 m MSS pixels were resampled to 30 m. Thus, the average pixel error was neglected.

Rectification error (Er): It is the square root of the mean error of the image rectification process [75, 76]. The rectification error in the level-2 Landsat products was found to be one fourth of a pixel, i.e., in the present case was decided as 7.5 m.

Total Positional Uncertainty: The total positional uncertainty (U_t) is the result of all errors that were previously estimated. It is defined as the square root of the sum of the squares of the sources of different errors [75, 76]. The formula of U_t is given as follows:

$$U_t=\pm \sqrt{\left(E_s^2+E_t^2+E_d^2+E_p^2+E_r^2\right)}$$E5

Where E_s is the seasonal error, E_t = tidal error, E_d = digitizing error, E_p = pixel error, and E_r = rectification error. The annualized uncertainty (U_a) was calculated using the square root of the sum of the squares of total positional uncertainty for each shoreline divided by the analysis period [75] as is given below.

$$U_a=\pm \frac{\sqrt{{\sum }_1^n{U}_{\mathit{ti}}^2}}{T}$$E6

Various uncertainties in the historical shoreline position between 1948 and 2021 is given in Table 3.

Uncertainty	Positional uncertainty		Measurement uncertainty			Total positional uncertainty
	Es (m)	Et (m)	Ed (m)	Ep (m)	Er (m)	Ut (m)
Landsat images	0	±2.12	±10	0	±5	11.37
SOI topo map	0	0	±15	0	±15	21.21
Army topo map	0	0	±15	±15	±30	33.54

Table 3.

Uncertainties associated with shorelines obtained from different sources.

The weight (w) is defined as a function of the variance in the uncertainty of the measurement (e). Weighted Linear Regression Rate (WLR) was computed using the total positional uncertainty values.

3.2.4 DSAS configuration

It consists of four main steps as is given below.

1. Baseline definition: An offshore baseline, at 200 m distance, almost parallel to the shoreline, was created in shape file format (.shp) with required attributes. The baseline is required to calculate the distance from shoreline to it at each orthogonal transect. The movement of all shorelines over the 48-year period (1973–2021) has been computed about the baseline. Due to the large uncertainty of the shoreline generated from Army topo maps, it was not considered for DSAS but used for future reference. When the shoreline moves landward to the baseline it is considered as erosion (denoted as negative values). In contrast, when the shoreline moves seaward, it is considered as accretion (denoted as positive values).

2. Collection of shorelines: Each shoreline vector represents a specific position in time and space, hence each of them is assigned a date (in dd-mmm-yy format) in the shoreline feature-class attribute table.

3. Generation of transects: Orthogonal transects were generated at 100 m intervals alongshore by using DSAS. At each of the representative locations, shore normal transects were generated against which relative changes in shoreline position were determined. A total of 1924 effective orthogonal transects were casted along the baseline from the south-west corner of the west bank to the south-east corner of the east bank in a clockwise direction and numbered (as Transects ID). The measurement transects that are casted by DSAS from the baseline, intersect the shoreline vectors. The points of intersection stores location and time information which are subsequently used to calculate the rate of change. The distances from the baseline to each intersection point along a transect are used to compute the distance and rate statistics. Some of the transects were removed from the analyses as they were falling on the confluence of the tributaries of the Hooghly River.

4. Calculation of change in shoreline position and other statistics: Two types of statistics are generated by DSAS viz. distance (NSM) and rate (EPR, LRR, and WLR). The regional change rate is calculated by averaging the rates of changes from all the transects. The average coefficient of determination (R²) and uncertainties of the annual rate-of-change (m yr.⁻¹) are computed at a 95% confidence interval (LCI95 or WCI95).

3.2.5 Prediction of tide condition

The date, time and height of tide were calculated using WXTide32 package. The height of tide is governed by the following harmonic equation given in the Manual of Harmonic Analysis and Prediction of Tides, special publication no. 98, US Department of commerce [77].

$$h=H_0+\sum _n^N{f}_n{H}_n\mathit{cos}\left[a_{n}t+{\left(V_0+u\right)}_n-K_{n^{\prime }}\right]$$E7

Where, h is the height of tide at any time t.

H₀ = the mean height of water level above datum used for prediction.

H_n = the mean amplitude of any constituent A_n.

f_n = the factor for reducing mean amplitude to year of prediction.

a_n = the hourly speed of constituent A_n.

t = the time, in hours, reckoned from beginning of year or prediction.

(V₀ + u)_n = the Greenwich equilibrium argument of constituent A_n when t = 0.

K_n’ = the modified epoch of constituent A_n.

N = the number or constituents used for the particular station.

In this equation except h and t, all other parameters are considered as the harmonic constant for any particular year and the place. Using these Harmonic constant, the successive value of tide height can be generated at any point of time. WXTide32 data pertaining to Diamond Harbor was considered to be representative of the present study. Table 2 provides the satellite overpass time versus the low and high tide timings.

3.2.6 Cluster analysis

Cluster analysis is a technique used to classify cases into groups that are relatively homogeneous within themselves and heterogeneous between each other, based on a defined set of variables [78, 79]. Hierarchical agglomerative clustering using the Ward linkage method was followed in the present study. In this method, clusters are merged to reduce the variability within the cluster. At every stage the average similarity of the cluster is measured. A case is selected to enter the cluster if the inclusion in the cluster produces the least increase in the error. The number of the cluster centres was determined from ‘Scree diagram’ in which ‘distance coefficients’ are plotted against the ‘stages’. The point at which there is a significant jump in the distance values was considered as the ‘elbow’ of the ‘Scree plot’. The numbers of clusters were decided as the number of cases minus the step of the elbow. Once the clustering is done, K-mean classification is performed for all the transects using the number of cluster centres from ‘Scree plot’. K-mean classification assign cluster membership and distance from the cluster centre to each case. Distance of the cluster centres are determined by using Euclidean distance as is given below:

$$d_{\mathit{ij}}={\left[\sum _{l=1}^q{\left(x_{\mathit{il}}-x_{\mathit{jl}}\right)}^2\right]}^{\frac{1}{2}}$$E8

Where dij = ED for two individuals i and j, each measured on q variables, Xij, Xji, i = 1, …q. ED (dij) is calculated as the sum of squared differences between relative cross-shore positions at each transect (T₁, T₂, T_3, etc.) during each epoch (a) and for all epochs (N). Smaller ED values indicate the cases are more similar. To evaluate the robustness of the clusters Kruskal-Wallis one-way ANOVA test [80] was carried out. A detailed description of cluster analysis can be found in Everitt et al. [81] and Hennig et al. [78].

4. Results and discussion

4.1 Shoreline configuration

The 200 km stretch of the study region has varied beach types including wide sandy beaches to mudflat, the mixture of sand and mud, mangrove wetlands as well as open mixed jungle at the backdrop of sandy/muddy beaches. The considerable length of the shorelines has embankments (Table 4). The western bank consists mainly of sandy and muddy beaches whereas the east bank predominantly consists of a muddy and mangrove systems with intermittent gap areas where the beach is absent. Zone-wise brief description of the beach configuration is given below.

Statistics	1973–1980	1980–1989	1989–1995	1995–2000	2000–2005	2005–2010	2010–2017	2017–2021
Mean	−14.32	−15.43	1.69	−43.00	39.79	−7.64	42.06	70.67
sd	200.03	86.59	77.90	110.99	122.09	100.35	135.17	162.57
Max	1386.98	566.43	681.41	629.64	1141.27	908.46	1000.23	1057.30
Min	−404.98	−527.71	−293.49	−1061.84	−313.23	−478.47	−519.50	−180.76

Table 4.

Mean shoreline change (m) over different time intervals.

Zone-1: This zone constitutes the western bank of Hooghly estuary and falls between the outlet of Pichhaboni Khal and Rashulpur river, confined between the transects T₂₅ and T₂₁₁. The foreshore is mostly sandy. Several sandy beaches (near Haripur, Boral, Bankiput and Digene) and the large number of aquaculture ponds are present in this zone. A large difference between high and low water lines (maximum of about 850 m) is observed between Pichhaboni Khal to the south of Gopalpur based on Survey of India topomaps. From satellite imagery, significant accretion and expansion of forest land can be seen between the Pichhaboni outlet and the north-east of Junput.

Zone-2: This zone is represented by the transect number T₂₁₈ to T₅₂₄ and falls between the confluence of the Rashulpur river and Haldi river on the west bank of Hooghly. The foreshore is mostly muddy with small sandy beaches near Hijli and Khejuri. Thick forest vegetation is observed in the east of Thanaberia along the coast. Large numbers of aquaculture ponds are present between Kandlamari and Khejuri. There is no embankment present in the northern half of this zone. The difference between high and low water lines is comparatively less (∼ 80 m) between Talpati Khal and Haldi river.

Zone-3: This zone extends from the Haldi river to the Rupnarayan river confluence, demarcated by the transects T₅₃₄ to T₉₁₇. This zone is highly convex towards the bay and the foreshore is muddy. The difference between high and low water line increases from Gilabaria (T₆₀₀) till Jhikarkhali (T₆₈₅) which further increases near Horkhali. Series of brick kilns and few aquaculture ponds can be seen in this zone. There is no protective embankment present is this zone.

Zone-4: This zone falls in the east bank of Hooghly estuary between the transect T₉₁₉ and T₁₁₇₄. Very narrow muddy foreshore in present here, however, in some places it is absent. Most of the river bank of this zone is embanked especially upto Kantabaria, however, from here till Kulpikata Khal there is no embankment. The narrow difference exists between high and low water lines except near Diamond Harbor. Series of brick kilns can be observed in the north (Simulbaria to Diamond Harbor) and south (Kantabaria to Kulpi) of this zone.

Zone-5: The zone covers the coastline between Kulpi and Kakdwip, defined by the transects T₁₁₇₅-T₁₄₅₄. The difference between high and low water lines is very small and sometimes absent. Almost the entire coastline has an embankment. Few brick kilns and aquaculture ponds are present in this zone.

Zone-6: This zone falls between Kakwip and Namkhana (T₁₄₅₅-T₁₅₇₄), totally embanked and shortest among all. The difference between high and low water lines is very small especially at south of Kakdwip and near Nadabhanga Khal. No brick kiln or aquaculture pond is present here.

Zone-7: This is the last and longest zone between Namkhana and Henry Island in the east bank of Hooghly, represented by the transect T₁₅₇₅-T₁₉₆₇. The foreshore consists mostly of a muddy area (Dakshin Durgapur) but a narrow sandy beach is present in the west of Namkhana and a wide sandy beach in the extreme south near Lakshmipur till the outlet of the Bakkhali river. In the south of Patibunia considerable area is under a mangrove swamp and open mixed jungle. Moderate to good forest exists in the extreme south of the region confined between T₁₉₁₂-T₁₉₆₇. Near Dakshin Durgapur, the difference between high and low water line is more (∼300 m). Most of the coastline is embanked. No brick kilns and aquaculture ponds are present in this zone.

4.2 Spatiotemporal change in shoreline

The large difference in the shoreline position was observed within each time interval and among different intervals. The dynamics of the shoreline are mainly due to disequilibrium in the morphological state and northward tapering nature of the estuary coupled with plausible subsidence due to auto-compaction of the Holocene sediments. One commonality among all the time intervals is the large variation in the seaward end of both the banks (Figure 3). During 1995–2000 and 2005–2010, the overall variation in the shoreline position is minimum. In comparison to the east bank west bank has more variation except for 1973–1980. Considering all the temporal intervals between 1973 and 2021 average recession is maximum in 1995–2000 (− 43 m ± 110.99) especially due to erosion in the southern part of the east bank of the estuary. In contrast, there is an increasing trend in the seaward extension of the shoreline since 2010. Between 2017 and 2021 the average accretion is 70.67 m (± 162.57). The maximum accretion length was 1386 m at T₁₄₈₆ (south-west of Kalinagar) in 1973–1980 whereas maximum erosion was −1062 m at T₁₈₆₅ (west of Fraserganj) during 1995–2000 (Table 4).

Figure 3.

Shoreline changes recorded at different transects over different temporal intervals.

The percentages of transects recorded aggradation or recession is given in Figure 4. From the figure, it is apparent that the proportion of aggradation and erosion does not match over the time intervals. The percentage of the transects exhibiting erosion was comparable during 2000–2005 (29%), 2010–2017 (30.20%) and 2017–2021 (27.81%). There was an abrupt increase in the erosion by 69.91% in 2010–2017. In general, there is a decreasing trend of erosion, especially after 2000 (Figure 4).

Figure 4.

Percentage of transcets showing erosion at different time intervals.

Figure 5 depicts how each zone contributes to the total shoreline change. Between 1973 and 2021, zone 5 contributed maximum towards erosion. Other zones that contributed marginally to erosion include zone 7 and zone 6. Zone 6 showed consistent erosion in all the intervals except for 1973–1980. Very high annualized aggradation of 69.17 m and 29.93 m was recorded in zone 1 and 2 respectively over the entire period of 1973–2021.

Figure 5.

Contribution of each zone towards erosion / accretion at different time intervals (Z represents the zones).

It is interesting to note that while comparing the coastline of 2021 with respect to 1948 (not used in the DSAS), there is a significant recession (∼900 m) in the zone 2 (between Talpati Khal and Kaldalmari) and in the zone 3 (near Horkhali) by about 600 m. In the east bank, most significant erosion is noticeable in zone 5, between Jadabnagar and Tilakmandal chak. The maximum landward retreat recorded was 2700 m near Uttar Chandannagar. On the other hand, accretion was observed in the south of zone 1 and 2 as well as in the north of zone 6. Quantitative analysis of the coastline change in this region has been carried out by Bandyopadhyay et al. [82], Raju et al. [83], Jana et al. [84], Rudra [85], Chakraborty [49] and Das et al. [86] along with their underlying mechanism. They have opined that beach erosion is attributed to various causes such as decrease of sediment supply from rivers, land subsidence, and interruption of longshore sediment transport by man-made structures. As the sea level rises, it causes waves to act on higher parts of the beach profile, resulting in enhanced erosion. If the sandy beaches disappear as a result sea-level rise, waves and storm surges, it will impact higher areas along the coastline [87].

Jana and Bhattacharya [88] used multi-resolution Landsat satellite imagery of 1972–2010 for shoreline change study along the 65 km long coastal stretch located between Rashulpur (Purba Medinipur) and Subarnarekha (Balasore) estuarine complex. The authors revealed that about 23 km of coastline recorded accretion, which was observed on several beaches such as at Talsari, Udaipur and Haripur, which were not affected by anthropogenic activities.

4.3 Shoreline change rate

The shoreline change rates were computed by linear regression and end point rate method at a lateral spatial interval of 100 m along the coast. The rates of changes of shoreline at different transect points estimated by EPR and LRR methods are given in Figure 6. Large variation in net shoreline movement and change rates were observed in the study region among various analysis zones (Table 5). Considering long term change between 1973 and 2021 four of the zones viz., 1, 2, 3 and 4 showed positive change (aggradation) by WRR method, the range of which varies between 0.24 m yr.⁻¹ (zone 3) to as high as 9.45 m yr.⁻¹ (zone 1). The very high recession was found in east bank at zone 6 (− 4.35 m yr.⁻¹), followed by zone 5 (−3.02 m yr.⁻¹). Overall, the rate of aggradation superseded the rate of erosion in the 48 years span. The zones which experienced maximum net seaward movement include zone 1 (553.34 m) whereas maximum net landward movement (erosion) of shoreline was found in zone 5 (137.22 m) (Table 5). There is a good agreement between both the methods (EPR and WRR) in respect of zone 2, 3, 4, 5 and 7. In zone 1 large difference in the shoreline change rate calculated by both the methods was recorded between Junput and Jagannathpur (T₁₀₂ to T₁₆₈), whereas in zone 6, the differences are significant in the region south of Kakdwip to Budhakhali along the Kakdwip river (T₁₄₅₅ to T₁₅₁₁).

Figure 6.

The rate of change of shoreline by WRR and EPR method.

Zones	SCE	NSM	EPR	WLR
Zone-1	740.83 ± 359.22	553.34 ± 298.84	11.48 ± 6.20	9.45 ± 6.22
Zone-2	391.0 ± 263.97	239.42 ± 307.19	4.97 ± 6.38	3.47 ± 4.89
Zone-3	170.61 ± 128.57	−16.81 ± 137.25	−0.35 ± 2.85	0.24 ± 3.28
Zone-4	81.31 ± 41.47	36.91 ± 54.96	0.77 ± 1.14	0.36 ± 0.93
Zone-5	249.19 ± 150.22	−137.22 ± 206.24	−2.85 ± 4.28	−3.02 ± 3.54
Zone-6	668.09 ± 351.75	−5.18 ± 456.16	−0.11 ± 9.47	−4.35 ± 5.36
Zone-7	303.21 ± 259.99	−2.30 ± 210.13	−0.05 ± 4.36	−0.38 ± 4.47

Table 5.

Zone-wise average shoreline change envelope (SCE), net shoreline movement and change rate by EPR and WRR method.

Although, the net shoreline movement (NSM) values are less in zones 3, 6 and 7 but the shoreline change envelope records large variation which indicates that the inter-annual fluctuation is very high in these zones and morphodynamic processes are very active.

Based upon the rate of erosion/accretion by WRR method, the transects were grouped into 7 classes (Table 6). From the table, it is evident that most of the shoreline (more than 73.33% by WRR and 69.95% by EPR) exhibit erosion/accretion rate between −5 and + 5 m yr.⁻¹. Low erosion rate (< 1.0 m/yr) was exhibited by 13.46% and 11.43% of the shoreline in WRR and EPR method respectively (not presented in the table). The proportion of very high erosion (<−10 m yr.⁻¹) and aggradation (>20 m yr.⁻¹) is limited to less than 2% of the shoreline. The spatial distribution of different change classes by WRR method is given in Figure 7a. It can be seen from the figure that in the west bank only one segment exhibits high erosion (−10 to −5 m yr.⁻¹) whereas in the east bank at least 6 segments (east of Kharibaria) show high erosion. This area exhibits has a large difference between low and high tide lines. While comparing with the Army Series map of 1948, it was found that there is a significant landward movement of shoreline between 1948 and 1973. In the east bank, there is no area under high erosion in zone 4, however, in zone 5, 6, and 7 considerable area along the shoreline is under high to very high erosion state. There are 3 distinct stretches near Uttar Chandannagar, Ramganunagar, Madhusudanpur and Lakshimipur. Close observation with the Army toposheet of 1948 reveals that there is an extensive recession in this area. The Rangatala island which used to be an integral part of the east bank has almost reduced to half between Kulpi and Madhusudanpur. The southern half of zone 6 has a high to very high rate of erosion between Budhakhali and north or Namkhana. The zone 7 is punctuated by two major areas of high erosion i) in the west of Edward creek, dominated by mangrove swamp and open mixed jungle and ii) in the east of Henrys island. In contrast to erosion, high to very high aggradation (> 20 m yr.⁻¹) is recorded between south of Gopalpur to Junput dominated by a wide sandy beach and inter-tidal difference. High aggradation is also observed in the south of the Rashulpur river confluence. In zone 2 high rates of accretion is observed in the north of Rashulpur river and east of Nij Kasba. In the east bank, there is no area of high accretion except in zone 7, near Lakshmipur dominated by mangrove swamps. This observation is in good agreement while comparing with the Army topo map of 1948.

Class	Range	WRR method		EPR method
		No of transects	% of total transects	No of transects	% of total transects
1	< −10	37	1.92	28	1.45
2	−5 to −10	187	9.71	204	10.60
3	−5 to +5	1411	73.33	1346	69.95
4	5 to 10	166	8.62	121	6.28
5	10 to 15	77	4.00	130	6.75
6	15 to 20	29	1.50	71	3.69
7	> 20	17	0.88	24	1.24

Table 6.

Different classes of erosion/accretion rates and their contribution to the shoreline.

Figure 7.

Shoreline changes a) rate of erosion/accretion (m yr.⁻¹) and b) change pattern.

4.4 Temporal pattern of erosion/accretion

To understand the temporal pattern of change direction, transects were grouped into two categories viz., eroding and aggrading type based upon displacement direction in each time interval. From 8 different temporal intervals 256 unique combinations were generated which were further grouped into 8 categories viz. i) consistent erosion (when in all the 8 temporal intervals the changes are negative) ii) mostly erosion (when erosion is recorded at least in 5-time intervals) iii) recent erosion (when last 3 or more consecutive intervals erosion is dominant) iv) mostly accretion (when accretion is recorded at least in 5-time intervals) v) recent accretion (when last 3 or more consecutive intervals accretion is dominant) vi) alternate (when erosion and accretion takes place alternatively) vii) trend reversal (changes from erosion to accretion over the years in a consecutive manner) and viii) others (when no definite trend is observed). Table 7 provides change patterns and their contribution to the entire shoreline under study. At a long-temporal scale, 1.87% of the shoreline shows consistent erosion which is alarming and another 36.69% are mostly eroded. We could not find any transect recording accretion in all the time intervals. All the 3 accretion types (MA, RA and TREA) together account for 38.93% of the shorelines (Table 7). Reversal of trend towards aggradation during last 3 or more consecutive time intervals was found from 4.73% of the shoreline, mostly located in zone 1 and 2. Only a small proportion of the shoreline (0.42%) exhibits alternate erosion and accretion over the years and does not yield a definite trend. Figure 7b shows the changing pattern along the shorelines.

Type	Description	Change direction	No of transects	% of total shoreline
CE	Consistent erosion	-ve	36	1.87
ME	Mostly erosion	-ve	706	36.69
RE	Recent erosion	-ve	10	0.52
MA	Mostly accretion	+ve	571	29.68
RA	Recent accretion	+ve	87	4.52
ALT	Alternate	mixed	8	0.42
TREA	Trend reversal (erosion to accretion)	+ve	91	4.73
OTH	Others	mixed	415	21.57

Table 7.

Temporal change pattern of shoreline behavior and their contribution.

Some of the transects that recorded both high erosion rate (more than 5 m yr.⁻¹) and consistent erosion are located in the north of Sibkalinagar (T₁₃₇₂-T₁₃₇₄), south of Budhakhali near Ghiya Khal (T₁₅₂₀-T₁₅₃₃), south of Nadabhanga Khal (T₁₅₅₂-T₁₅₅₇) and north of Duaragra Gang in zone 6 (T₁₅₆₉-T₁₅₇₄).

4.5 Hierarchical agglomerative clustering

Although, shoreline change analysis quantifies rates and directions of change, further analyses are needed to resolve distinct modes of coastal system behavior. Traditional shoreline changes analyses quantify the rate and direction of change by analyzing multi-date/historical data. However, there are some commonalities in terms of coastal system behavior. The Hierarchical agglomerative clustering was performed using the change matrix of all 8 temporal intervals to define the distinct coastal change behavior. Clustering was done using the Ward method which computes the sum of squared distance within the clusters and aggregates the clusters with the minimum increase in the overall sum of squares. The distance coefficients were plotted against the stage to generate a ‘Scree diagram’ (Figure 8). The number of clusters in the present study was 5 which was used for K-mean clustering. The cluster centres and the distances between cluster centres are given in Tables 8 and 9 respectively.

Figure 8.

Scree diagram defining the optimum number of clusters using elbow rule.

Temporal intervals	Clusters
	1	2	3	4	5
1973–1980	−44.37	−52.95	−64.31	926.40	193.82
1980–1989	−6.50	−64.44	−51.74	−76.24	226.13
1989–1995	−.32	24.22	1.62	−57.19	82.35
1995–2000	−25.80	−89.62	−35.18	−50.58	−931.35
2000–2005	27.71	97.21	−14.32	−59.69	941.45
2005–2010	−21.90	85.34	48.19	−83.01	−3.82
2010–2017	32.63	−28.49	464.74	−87.87	22.23
2017–2021	12.53	341.60	445.03	−54.35	109.87

Table 8.

Various cluster centres.

Cluster	1	2	3	4	5
1	—	369	619	991	1335
2	369	—	521	1087	1278
3	619	521	—	1249	1476
4	991	1087	1249	—	1572
5	1335	1278	1476	1572	—
6	—	369	619	991	1335

Table 9.

Distance between cluster centres.

The clusters captured a unique pattern of change at a temporal scale (Table 8). Among all the transects, 79.15% are represented by cluster 1 and only 0.94% by cluster 5. In clusters 1, 2 and 3 most of the transects show a balancing act of aggradation and erosion at different temporal intervals. The transects that recorded consistent erosion (Figure 7) were found in cluster 1 only. In cluster 4, erosion is dominant, while in cluster 5 accretion is dominant in most of the time span. The mean displacement of the shoreline in cluster 1 is −3.25 m and the maximum is 99.25 m in cluster 3, constituting only 4.80% of the total transects. All the clusters show aggradation in terms of their mean displacement values except cluster 1.

5. Discussion

Beach profile morphology and coastline, change over a range of time and spatial scales. The short-term variability occurs over a period of days to a month as a result of i) episodic events (storms) ii) medium-term variability over several months (e.g., winter summer wave change) to several years (e.g., due to regional climate variability, engineering intervention and prevailing sedimentary processes) and iii) long term variability that occurs over a period of a decade to a century, associated mainly with climate change impact; and very long term millennial-scale evolution as a result of quaternary sea-level changes [89]. Broad-scale analysis of changes in shoreline position has the potential to highlight the role of regional forcing on large-scale coastal behavior, e.g., long-term tidal cycles [90] or sea level rise [4]. Shoreline change analysis is also useful to identify notable ‘hotspots’ of contrasting behavior [91, 92]. The Hooghly estuarine shoreline analyses studied here comprehend synthesis of historical shoreline change over 48 years supported by limited ground observations. The data has been analyzed at high spatial resolution (100 m, alongshore interval) along the entirety of a 200 km shoreline. In the area evidence for strong met-ocean forcing is ostensibly compelling. The phenomena of erosion and accretion are largely regulated by littoral current patterns and sediment influx from different rivers and the adjacent Bay of Bengal. The shoreline of this 200 km stretch has different configurations from the sandy beach to muddy swamp punctuated by anthropogenic footprints including brick kilns, aquaculture ponds, protective embankments and beach nourishment treatments. Beach nourishment projects and coastline protection structures can result in an artificial accretion of coastline in a short period [93]. Large variations exist in shoreline position within the same year and also among different years indicate the disequilibrium in the morphological state. There could be several external factors responsible for shoreline change including sea level rise, changes in the wave climatology and storm intensity as well as changes in the catchment characteristics due to deforestation and land degradation which results in higher sediment load in the terrestrial run-off. In contrast to surface runoff, engineering intervention through the construction of dams and barrages also makes the estuary sediment starved. In long-term perspective, temporal data of PSMSL (Diamond Harbor and Haldia) reveals that the sea level is rising at the rate of 2.41 and 3.02 mm yr.⁻¹ respectively. The sea surface temperature induced El Niño-Southern Oscillation (ENSO) has a significant role in global atmospheric circulation influencing the temperature and precipitation. The irregular pattern of El Nino and La Nina triggers rainfall variability over the Indian sub-continent. In recent years strong La Nina and very strong El Nino have been witnessed in 2010–2011 and 2015–2015 respectively. The monsoon rainfall variability has a direct relation with terrestrial run-off and estuarine water level. Since 1951 there were 8 strong to very strong El Nino and 7 strong types of La Nina years. The storm surges are another strong forcing factor in a short temporal scale that can change the shoreline configuration. Although, the frequency of cyclonic storms is declining over the Bay of Bengal but the intensity is increasing. Extremely severe cyclonic storms of 2019 and 2020 are the best examples causing extensive damage to the coastline embankments. Karunarathna et al. [89] found single storms or storm clusters predominantly change the supra tidal and inter-tidal part of the beach profile and that beach erosion volumes are strongly correlated to the power of the storm. Once the astronomical tides coincide with storms, extreme sea level occurs resulting in large-scale inundation and damage to the coastal structures. Besides warming of sea surface relative, sea level change can also happen due to vertical land motion that can result from glacial isostatic adjustment, tectonic processes, coastal subsidence and uplift caused by anthropogenic factors. High-frequency and short temporal scale sea level variability due to seiches, meteotsunamis are frequently under-represented in sea level studies and yet contribute to the extreme sea levels which are of great research interest and importance to coastal dwellers [94]. In general, coastal landforms affected due to short-term perturbations viz. cyclone generally attains a morphodynamic equilibrium often by adopting different ‘states’ in response to varying wave energy and sediment supply [95]. Nevertheless, elevated surge water levels are known to be important drivers of longer-term dynamics on sedimentary shorelines [96]. Besides sea level rise and storms alongshore sediment transport can also have an impact on the coastline change, in particular, it is likely to result in coastline accretion.

Most of the west bank of Hooghly estuary is prograding at the rate of 0.24 m yr.⁻¹ in zone 3 to as high as 9.45 m yr.⁻¹ in zone 1. Whereas recession is pre-dominant in the east bank, especially in zone 5, 6 and 7 accounting −0.38 to −4.35 m yr.⁻¹. In general, aggradation dominates over erosion. Large variation in the shoreline change envelope in zone 3, 6 and 7 reveals an active morphodynamic process. The different suite of behaviors in recent intra-decadal scale suggests that forcing of coastal change can be interpreted as a form of the time-dependent complex response of the kind envisaged by Schumm and Lichty [97] whereby changes over shorter time scale, are inherently associated with tighter cause-effect linkages at smaller spatial scales, and broader trends emerge over longer time-scales. Additionally, the phenomena of erosion and accretion are largely regulated by littoral current patterns and sediment influx from different rivers and the adjacent Bay of Bengal. The west bank of the estuary having sandy inter-tidal plain is aggrading over longer time scale whereas several areas in the east bank of muddy beaches record the high rate of erosion. The temporal pattern of erosion/accretion has been captured using the direction of change in each time interval. Some portions of the shorelines especially north of Kakdwip and Namkhana recorded a consistent high rate of erosion (> −5 m yr.⁻¹) over each interval. Although, only 1.87% of the area of the shoreline showed consistent erosion for all the time intervals but together with ‘mostly erosion’ type it constitutes 38.56% which is alarming. These areas need to be protected from anthropogenic intervention and to be stabilized by rejuvenating protective embankments or vegetative barrier. Contrasting modes of prograding stretch adjacent to retreating stretch can be found in close proximity, particularly in zone 1 and 2, which suggests that local influences may be particularly important. Both these transitions in behavior suggest localized net littoral fluxes of sand and gravel from the north of estuary to the south-west. These localized instances of coupled behavior have led to a distinct net change in regional shoreline planform over a longer time scale. Some of the stretches of the shoreline exhibit distinct change of cuspate foreland from rounded to sharp apex especially north of Jhikarkhali and Madhusudanpur, north of Kakdwip. Erosion at the north and progradation at the west and south-west, illustrates south-west transport of sediments over the studied time scale whereas diffusive behavior dominated decadal-scale shoreline change.

The inter-temporal analysis using spatial smoothing windows of 1000 m showed that there is no consistent association between convexities/concavities and the erosion/accretion. Some concave stretches of shoreline exhibit erosional signatures, whilst others are accretional. The convexity of the shoreline near Horkhali (in the west bank) increased over time but decreased near Madhusudanpur on the east coast, however near Kakdwip and Patibunia the convexity remained almost unchanged over the years. Some of the concave stretches of the shoreline showed seaward accretion in the west bank, e.g., at Nij Kasba. The eroding sediments move parallel to the coast by alongshore currents from north to south direction and are expected to deposit around the concave coast owing to the lower current velocity [93]. As a result, the coastline can advance to the ocean around these regions. Several studies claim that a concave-shaped coastline tends to exhibit accretion while a convex-shaped coastline tends to exhibit erosion [93]. However, in the present study, several concave stretches of the east coast exhibited landward retreat of coastline typically along the Rangafala channel near Lakshmipur, between Ghiya Khal and Duraragra Gang (north of Namkhana) and small patches in Patibunia island. Presumably, both diffusive and anti-diffusive (unstable) behavior is operational [98] which are likely to change as the shoreline planform adjusts in response to the consequent patterns of erosion and deposition.

With the anticipated increase in global mean temperature by about 0.5°C, the thermal expansion and melting of ice caps and glaciers are inevitable [13] but this effect may be masked by inlet dynamics and coastal engineering projects even over extended time periods. However, the implication is that sea level rise is a secondary but inexorable cause of beach erosion in such areas which may lead to high-energy swells to reach further up the beach and redistribute sand offshore. Apart from the external and natural forces there are alarming uncontrolled anthropogenic activities which have imposed excessive pressure on the coastal landuse and exacerbating beach erosion problems along the Hooghly estuary. This will have ominous implications for ever-increasing coastal population and associated livelihood [99]. There is a need for decoupling the long-term forces from the anthropogenic effects and projecting the future scenario of coastal changes for effect coastal planning and enforcement.

6. Conclusion

The study of historical evolution and sculpturing of the coastal areas of Hooghly estuary in terms of short and longer time scale has significant importance in evaluating the criticality in shoreline change. The findings of the present study revealed that geospatial techniques are very useful for analyzing and predicting shoreline dynamics. The short- and long-term changes have been estimated using the DSAS extension tool of ArcGIS. The tool enables the calculation of several change metrics and also the rate of changes from time-series shoreline positions and helps in determining the zones of erosion and accretion. The variation is higher in the west bank than east bank except for 1973–1980. Considering the entire study period average recession is maximum in 1995–2000 (−43 m ± 110.99) especially due to erosion in the southern part of the east bank of the estuary. Zone 5 contributed maximum towards erosion, however, in general, there is a decreasing trend of erosion, especially after 2000. While comparing with 1948-topomaps there is a significant recession (∼900 m) in zone 2 (between Talpati Khal and Kaldalmari) and in zone 3 (near Horkhali) by about 600 m. On the east bank, the most significant erosion is noticeable in zone 5, between Jadabnagar and Tilakmandal chak. The maximum landward retreat recorded was 2700 m near Uttar Chandannagar. The shoreline erosion is attributed to various causes such as decrease of sediment supply from rivers after construction of barrages in the upstream, land subsidence due to natural compaction or extraction of ground water, interruption of longshore sediment transport by man-made structures and dredging operation to maintain the navigation channel. In contrast, there is an increasing trend in the seaward extension of the shoreline since 2010. Between 2017 and 2021 the average accretion is 70.67 m (±162.57). Very high annualized aggradation of 69.17 m and 29.93 m was recorded in zone 1 and 2 respectively over the study period. The shoreline change rate computed using WLR method reveals that zone 1, 2, 3 and 4 show the positive change (aggradation) which varies between 0.24 m yr.⁻¹ (zone 3) to as high as 9.45 m yr.⁻¹ (zone 1). The very high recession was found in the east bank in zone 6 (− 4.35 m yr.⁻¹), followed by zone 5 (−3.02 m yr.⁻¹). More than 73% area of the shoreline exhibits erosion/accretion between −5 and 5 m yr.⁻¹. The proportion of very high erosion (< −10 m yr.⁻¹) and aggradation (> 20 m yr.⁻¹) is limited to less than 2% of the shoreline. The temporal change pattern was examined using change direction in each time interval. About 1.87% of the shoreline shows consistent erosion in which at all the time-interval the direction of change was negative and an additional 36.69% constitutes of mostly eroded, characterized by erosion in at least 5 epochs. There is no area where consistent accretion was observed. In about 4.37% of the shoreline trend reversal from erosion to accretion has been observed. The change rate and pattern maps generated in the study will be helpful for policy makers to prepare a strategic coastal management plan and for future policy intervention. It is suggested that there should have a regular monitoring mechanism of this estuarine region to keep watch on the shoreline change and triggering factors and regulatory purpose.

Acknowledgments

The authors are thankful to the Chief General Manager, RRSCs (NRSC) for his keep interest and sustained support to carry out this study. Thanks, are also due to earthexplorer.usgs.gov for providing satellite data freely to the user community.