Remote Sensing of Narrowing Barrier Islands along the Coast of Pakistan over Past 30 Years

Barrier islands (BIs) are the first line of defense against the sea/wave actions in coastal areas, and assessing their stability is crucial in the context of effective coastal planning. Therefore, this study evaluates the spatial–temporal shoreline changes of the BIs in Pakistan over the past three decades (1989–2018). Satellite data from Landsat missions are used to delineate the shorelines of 19 BIs in Pakistan. After delineating the shorelines from satellite observations, two well-known statistical methods (i.e., end point rate (EPR) and linear regression rate (LRR)) are used to capture the localized changes in the BIs. The results ascertain that nearly all of the BIs have experienced noteworthy erosion during the past three decades. While the mean erosion over all the BIs during the study period is estimated to be >6 m/y, significant spatial heterogeneities among the individual BIs exist. The interdecadal comparison indicates that the highest mean erosion of the BIs occurred during the period 1989–1999 (13.03 ± 0.62 m/y), which gradually reduced over the preceding decades (i.e., 7.76 ± 0.62 m/y during the period 1999–2009 and 3.8 ± 0.7 m/y during the period 2009–2018). Nevertheless, ~65% of the total BIs experienced high erosion (>2 m/y), ~15% experienced moderate (<2 m/y), and ~20% experienced low erosion (<1 m/y) during the period 1989-2018. This situation implies that while ~65% of these BIs need immediate interventions from the concerned authorities, the 15% BIs with moderate erosion might experience high erosion in the wake of rising sea levels and decreasing sediment influx in the near future without proper measures. This depletion of the BIs might not only affect Pakistan but also have regional consequences due to their various services.


Introduction
Due to the significance of barrier islands (BIs) as the first line of natural defense to protect coastal communities against storm surges and provisioning of multiple services, considerable interest surrounds their health/stability for a better regional environment [1][2][3]. BIs provide essential economic services to coastal populations (e.g., supporting tourism) as well as serve valuable ecological functions (e.g., providing natural habitats to several species, such as migratory birds) supporting regional coastal sustainability. However, at the same time, these fragile and dynamic coastal features are also vulnerable to degradation from both natural processes and human interventions [4][5][6]. In the wake of climate change-induced sea-level rise and increasing anthropogenic pressures on coastal systems, it is essential to track the changes in the geometry, area, and position of BIs, particularly those bordering the developing and low-lying deltaic coasts. This information on coastal changes could help in taking appropriate measures for coastal resilience and sustainability [7].
Shoreline position determines the protection provided by various coastal barriers and landforms. The spatial-temporal shoreline change can be used as a proxy to monitor the coastal environmental change, as the assessment of long-term trends in shoreline changes is essential to improve the understanding regarding coastal responses to rising sea levels [8]. Several methods have been used for studying the shorelines, such as unmanned aerial systems (UAS), aerial photogrammetry, topographic and global positioning system (GPS) surveys [9][10][11][12]. Despite having a higher spatial resolution, they are money-and labor-intensive, particularly for the long-term monitoring of large-scale assessments for developing coasts. In the wake of climate change, it is imperative to track the shoreline changes [5], particularly those bordering the developing and low-lying deltaic coasts, to take appropriate measures for their stability. The information on spatial-temporal shoreline changes is of particular importance to coastal planners and decision makers [12,13]. To support effective policy production for coastal management, it is imperative to have baseline scientific information on the spatial-temporal trends and rates of shoreline changes [4]. Previous studies conducted in different parts of the world, such as Bangladesh [14], Egypt [15], the Mediterranean sea [16], the coast of North America [17,18], Mississippi River Deltaic Plain [19], and deltaic regions of China [20,21], highlight the need as well as the importance of monitoring the changes in the dynamic coastal areas for effective coastal policy and natural resource management. Similarly, some researchers map the perimeters of BIs [22][23][24] and quantify their geomorphological changes to inform national/sub-national coastal planning. On the contrary, few studies [25][26][27][28] have focused on studying shoreline degradation and morphological changes, particularly erosion, in the Indus River Delta, Pakistan-the world's 6th largest delta-on a large scale. Most of the existing studies are limited in terms of spatial-temporal scales and have primarily focused on the mainland shoreline [25,27,[29][30][31][32], neglecting the BIs in the Indus Deltaic Region (IDR). One seminal study by Waqas et al. [33] analyzed the changes in the area and perimeter of the BIs along the shoreline of Pakistan using the pixel frequency method. Relatively little is known about the rate at which these BIs have evolved and eroded. Furthermore, the approach adopted by Waqas et al. [33] is not adequate to provide localized changes regarding coastal erosion, which is a very important aspect of coastal adaptation. This lack of BI shoreline change assessment in Pakistan not only hinders the efficient coastal spatial planning and management but limits the coastal adaptation actions and opportunities to gain benefits in the form of many services provided by them. Furthermore, this lack of information on the rates at which these BI shorelines change restricts the consideration of these BIs in the overall coastal planning and management process. This situation can potentially influence the overall coastal sustainability in Pakistan and beyond. Hence, a comprehensive assessment of the coastal dynamics of the BIs in Pakistan becomes essential. In the context of current and future coastal planning processes and management strategies, this study aims to provide comprehensive preliminary baseline information on the spatial-temporal shoreline changes along with the BIs of Pakistan.
While the integrated utilization of remotely sensed earth observations and geographic information systems (GIS) [9,10,12,34] is becoming an essential tool for the management (monitoring and planning) of coastal environments [28,35], developing countries still lack in adapting these tools for several reasons, including data unavailability and higher monitoring coast associated with these labor-intensive and expensive methods. These tools, well known for spatial-temporal change assessment, serve as essential decision support systems vital for corrective local decision making and smart resource allocation to assure coastal sustainability [36]. The factors mentioned above make a large scale and long-term coastal assessment of shoreline changes in developing regions, particularly along with the BIs, unavailable or very limited. The lack of information regarding the rate at which BIs are translocating or their perimeters are being eroded further makes it difficult for developing coastal communities to take appropriate measures for their safety. Therefore, this study aims to use the multitemporal freely available medium resolution satellite observations from the past 30 years (1989-2018) in a GIS environment to systematically evaluate and map the shoreline changes in the BIs of Pakistan. Among two widely used shoreline change assessment methods (i.e., areal-and transect-based) [33,37], this study employs the transect-based assessment approach, as it provides comparatively detailed localized information on the shoreline changes. The results from this study would have significant policy and decision-making implications to aid effective coastal planning and management, which is integral to sustainability in the long-term. Additionally, these preliminary insights into the shoreline changes along the BIs would further help us to understand the impacts of human interventions (i.e., urbanization and intensive coastal development projects) on coasts, providing broader opportunities for sustainable coastal development. These temporal assessments help us to analyze the effectiveness of the previously implemented coastal policies, providing opportunities to update/revise them accordingly.

Study Area
In Pakistan, BIs are located in the form of a long chain along the south-eastern to the southern shoreline adjacent to the Arabian Sea along which the sea level is rising with a rate of > 1.1 mm/y [27]. About~1000 km shoreline of Pakistan is divided into Sindh and Makran coasts. This study is focused along the~320 km Sindh coast. The Sindh coast is further divided into two distinct parts encompassing 5 BIs along the Karachi coast and 14 along with the 6th largest Indus Delta Region (IDR, Figure 1 and Table 1). The Karachi shoreline is the most developed segment along the entire shoreline of Pakistan-presenting the highest anthropogenic/urbanization pressure on this part of the shoreline. However, the IDR shoreline is the most undeveloped segment where these BIs are separated by a unique dendritic network of 17 major creeks, several minor creeks, estuaries, and tidal inlets. Using the tidal classification scheme from Davies [38], the shoreline in the study area is found to be under the influence of mixed tidal conditions (mesotidal 2-4 m-Karachi shoreline, and macrotidal > 4 m-up to Sir creek) [39,40]. Furthermore, the BIs are subjected to a mixed energy regime (western side of IDR as wave-dominated and the eastern side as the tide-dominated) with an annual mean significant wave height (SWH) in the study area of 1.8 m during SW-monsoon season and 1.2 m in NE-monsoon season [39]. Any changes in the energy regime and the wave/tidal height can reshape these unprotected BIs [41]. The wetlands in this deltaic region are protected under the Ramsar Convention's international regulations. The abrupt reduction in freshwater discharge from the Indus River to the Arabian Sea and the sediment loads reaching these natural coastal barriers are contributing to their degradation.

Data Used
For the long-term assessment of the BIs evolution in the study area, the following data sets were used for processing, acquisition of required information for BI perimeter estimation, and interpretation of the results (Table 2).

Satellite Data
This study used freely available and atmospherically corrected 30 m archived Level-2 (surface reflectance) imagery of Landsat satellite (TM: Thematic Mapper and OLI: Operational Land Imager) available from the United States Geological Survey (USGS)'s Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On-Demand Interface. The delineation of multitemporal perimeters (shorelines) of the BIs from remote sensing images is influenced by tides [27] and medium to high energy waves [42]. Therefore, we took great care to acquire the cloud-free Landsat images in the non-flooding pre-monsoon period (December-April).

Data Used
For the long-term assessment of the BIs evolution in the study area, the following data sets were used for processing, acquisition of required information for BI perimeter estimation, and interpretation of the results ( Table 2).

Satellite Data
This study used freely available and atmospherically corrected 30 m archived Level-2 (surface reflectance) imagery of Landsat satellite (TM: Thematic Mapper and OLI: Operational Land Imager) available from the United States Geological Survey (USGS)'s Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On-Demand Interface. The delineation of multitemporal perimeters (shorelines) of the BIs from remote sensing images is influenced by tides [27] and medium to high energy waves [42]. Therefore, we took great care to acquire the cloud-free Landsat images in the non-flooding pre-monsoon period (December-April).

Tide Height Data
Examining the tidal conditions for multitemporal images over the study region from ã 30 year (1989-2018) period showed that the tidal height ranges from~0.20 to 3 m. Based on the data availability in the selected season/period and mean tidal height value constraint (<2.8 m), four out of 17 images acquired over 152/043 path/rows (i.e., from the years 1989, 1999, 2009 and 2018) were selected for the final processing and further interdecadal analysis of the morphological assessment and evolution of the BIs.

Hydrodynamic Data
The seasonal effects of the wave and wind is an essential and strong oceanic factor affecting the nourishment and erosion of the BIs in the study area. Wind creates waves in the ocean, which in turn increases the water level and causes it to penetrate in creeks and BIs in the coastal region. Waves sort the sediments, erode the coast, and transport sediments in and out from one location to other. They also move sediment from near the shore to offshore. Wave set up during summer monsoon along the Pakistan coast makes the water level up to half a meter high (reaching up to 3.5 m in the creeks). The monthly mean (January to December) wind field was adopted from the National Centers for Environmental Prediction (NCEP) to represent the yearly influence and direction of wind-driven waves and currents. To analyze the meteorological forcing influencing the Pakistan coast seasonally, the effects of wind on waves at different tidal stages (spring, neap, high, and low) were observed through the Simulation WAves Nearshore (SWAN) model [43] simulation using regional data from the National Institute of Oceanography (NIO) and analyzing the results for different wind conditions in Near Shore Wave (NSW) modules in MIKE software, powered by the Danish Hydraulic Institute (Denmark) [44]. For the estimation of SWH, we chose December as the representative month for the northerly wind case (January, February, November, and December), April for the westerly wind case (March and April), and June for the southwesterly wind case (May, June, July, August, September, and October).

Sediment Discharge Data
The sand supply from upstream also determines the stability of BIs. A deficit in the sediment discharge/flow reaching the shoreline adversely affects the sediment supply and redistribution to the BIs driven by the combined action of wind and waves, hence their stability. Long term sediment discharge from the Indus River below the Kotri barrage from Kidwai et al. [45] was used to evaluate the influence of the hydrodynamics of the Indus River (IR) on the erosion of the BIs in the study area.

Sea Levels Data
The menace of the rise in sea levels on the Pakistan coast was expected to be more pronounced in the IDR. To analyze its adverse effects on the BIs stability, the long-term historical records of sea levels at Karachi and adjoining IDR were collected during the study period from the Permanent Service for Mean Sea Level (PSMSL) data portal.

Methods
The research work was performed in a series of steps, from data acquisition to data selection, pre-processing, shoreline delineation, change rate analysis, and mapping. For simplification, we divided the overall workflow into two steps. The first step consisted of data acquisition, selection, pre-processing, and perimeter (shoreline) delineation for each BI. The second step included an assessment of changes in the shorelines of each BI during the period 1989-2018 using the transect-based method in the Digital Shoreline Analysis System (DSAS). The DSAS is an ArcGIS tool developed by the USGS and is widely used in coastal assessments to aid policy and decision making in the context of coastal management [15,46]. The spatial assessments were performed to provide a geographical reference and information for the prioritization of different BIs for immediate or gradual measures. The details on the two steps are provided in the following sub-sections.

Step 1: BI Perimeter (Shoreline) Delineation
Most of the shoreline along the BIs in Pakistan can broadly be classified as the muddy shoreline. Details of the classification scheme can be found in Zhang et al. [24]. For muddy coasts, the seaward edge of vegetation is identified as the BI perimeter (shoreline) [20]. The following steps were attempted to precisely delineate the perimeter of the BIs: (a) suitable band selection using a spectral profile curve of the images; (b) histogram threshold shoreline delineation; (c) edge detection; (d) conversion of binary data to shoreline vectors (see Step 1 in Figure 2).

BI Perimeter Delineation
The BIs, along with the Sindh coast, have experienced migratory shorelines over the past three decades ( Figure 3). The image-based BI perimeter uncertainty is subjective and  To check the accuracy of the results, the vector shoreline representing the boundary of each BI was examined with the corresponding false-color combination (bands 2, 3, 4 for TM, and 3, 4, 5 for OLI 8) of Landsat images. This overlay analysis helped us visually verify the alignment of the BI perimeters with the seaward edge of vegetation representing BI boundaries on Landsat images and to ensure the quality/reliability of the results, as practiced by Sajjad et al. [7]. BI-wise visual assessment also enabled us to identify which islands are experiencing disintegration and where new sandbars/BI-segments are appearing over time.

Step 2: BI Perimeter Change Rate Analysis and Implications to Coastal Planning and Development
For change rate analysis, all the delineated shorelines (one for each time-step) were appended to a single shapefile, separately, for each BI, which was later utilized in the DSAS. To ensure that seaward shorelines were only compared with the seaward shorelines and not the landward shorelines, two versions of shoreline coverage representing the BI perimeter were created (i.e., open-ocean side of BI and bay/landside of BI) [47]. We used several editing tools in ArcGIS, including Delete and Merge, to prepare the BI perimeters by removing the unwanted extensions for this study. Using the polygon coverage of the extent of each BI, we clipped larger shorelines for each BI and saved them along with the corresponding baseline segment in its workspace.
After BI perimeter (shoreline) delineation, the following sequential steps were carried out to analyze the evolution and change rate of the BIs: (a) baseline generation; (b) transect generation; (c) computation of distances between BI perimeters and baseline at each transect; and (d) estimation of shoreline change rate statistics. To establish a reference baseline for change rate analysis, a 50 m buffer was created around the oldest shoreline (the year 1989). The baseline for each BI (both main island and neighboring small islands) was segmented into seaside and bay/landside following the overall sinuous shape [28] of the BI perimeter as best as possible [47]. The same process was repeated for all the shorelines of the BIs to create seaside and bay/landside shoreline segments. The division into seaside and landside was made for those parts of BIs which had shifted entirely. This method enabled us to study the BI migration on both sides, in addition to the shoreline change rate. All these data were then imported and managed in the Personal Geodatabase in ArcGIS 10.3. Each BI was analyzed separately, and transects of variable length were cast every 100 m orthogonally along the baseline considering their perimeter, geometry, the direction of translocation, and the extent of migration. Two different methods, endpoint rate (EPR) and linear regression rate (LRR), were used to compute and analyze the longterm and short-term changes in the BIs during the period 1989-2018. With reference to the baseline, the accretion of the shoreline was considered as a positive value and the erosion as a negative value. The results were spatially quantified and classified based on the shoreline migration direction and the rate during the study period. To determine whether erosion rates have accelerated or decelerated temporally, BI perimeters were analyzed using the EPR method for each distinct interval (1989-1999, 1999-2009, and 2009-2018). The confidence of EPR (ECI) was calculated as: where U(A) and U(B) are the uncertainty values for shoreline A and B, respectively; Date A and Date B are the dates of shoreline A and shoreline B, respectively. The uncertainty of the calculated LRR rate, also known as the standard error of the slope of the regression line, is reported as a Confidence Interval of Line of the Regression (LCI). The LCI was considered at a confidence interval of 95% (i.e., LCI95).
To quantitatively assess the vulnerability of the shorelines to erosion in the study area, we adopted a relatively simple but robust classification scheme from Gornitz et al. [48] and divided the BI change rate statistics into 5 classes (i.e., >+2.0 m/y, +1.1 − +2.0 m/y, −1.0 − +1.0 m/y, −1.1 − −2.0 m/y, and <−2.0 m/y). This categorization helped us to prioritize the BIs for immediate or gradual actions and measures, where needed. A quantified histogram map of erosion/accretion intensity was produced for each BI. The interdecadal migration of the Bis was further analyzed using the historical trend of the changes in the BI perimeter distance from the reference baseline (year 1989).

BI Perimeter Delineation
The BIs, along with the Sindh coast, have experienced migratory shorelines over the past three decades (Figure 3). The image-based BI perimeter uncertainty is subjective and probably more difficult, especially when no reference dataset is available for such a large area in undeveloped deltaic regions. In order to evaluate the BI perimeter accuracy, initially, we searched for high-resolution data, such as the reference shoreline or freely available high-resolution satellite images (e.g., Sentinel-2), that coincide with the dates of the Landsat images used. No information was found for the same date. Therefore, the results of this study were verified by authors through personal communications with the local authorities and residents wherever possible [33].
A sequential comparison of the shapes and positions of the BIs revealed that the systematic pattern of the land-loss due to erosion is barrier narrowing, also known as depletion of the BIs (e.g., BI-8 and BI-9). The disintegration of the main island into several small sandbars surrounding the main BI (island segmentation) was observed in some cases (e.g., BI-8 and BI-12 in Figure 3). The corridor created between different island segments due to their disintegration could further give way to increased erosion during storm surges and strong wave-action. Sea-sides of many of these BIs were substantially eroded during the study period. It was observed that most of the BIs have evolved differently (i.e., majority of the BIs are shrinking predominantly on the seaside, some are transgressing landward, and very few are expanding towards the seaside). For some islands, it was not possible to calculate the rate of change throughout the entire perimeter of the BIs because some parts of the BI vanished (e.g., BI-18), and some new areas appeared (e.g., . BI-5 in the Karachi sub-region and BIs 8, 10, 12, and 14-19 in the IDR were disintegrated into smaller units, and some segments appeared later due to new sand deposits/accumulation ( Figure 3). Division of the BIs happened in all the islands located on the east side of the IDR, which formed the frontline of the sparse mangrove cover.

Rate of Shoreline Change in Different BIs
Out of the 19 BIs along the Sindh coast (including the Karachi and Deltaic coasts), 16 were found to be depleting, which calls for intervention and appropriate measures from

Rate of Shoreline Change in Different BIs
Out of the 19 BIs along the Sindh coast (including the Karachi and Deltaic coasts), 16 were found to be depleting, which calls for intervention and appropriate measures from the concerned authorities to be taken on a priority basis (Figure 4). A summary of the overall BI change rate and morphological evolution from 1989 to 2018 is presented in Figure 4. Only one sand body in BI-18 was found to be translocating seaward (seaward accretion and landside erosion), which lies in the eastern IDR containing sparse mangrove cover. It should be noted that for some BIs, the rate of change at a few transects was not computed because the direction of transect casting is perpendicular to the baseline, which does not intersect the shorelines at some places.
the concerned authorities to be taken on a priority basis (Figure 4). A summary of the overall BI change rate and morphological evolution from 1989 to 2018 is presented in Figure 4. Only one sand body in BI-18 was found to be translocating seaward (seaward accretion and landside erosion), which lies in the eastern IDR containing sparse mangrove cover. It should be noted that for some BIs, the rate of change at a few transects was not computed because the direction of transect casting is perpendicular to the baseline, which does not intersect the shorelines at some places.  While the overall mean erosion of the BIs during the period 1989-2018 was estimated to be > 6 m/y (6.9 and 6.4 m per year as per the EPR and LRR methods, respectively), significant spatial heterogeneity among different BIs exists ( Table 3). The Karachi shoreline is the most developed segment along the Sindh shoreline, which is invaginated by four major inlets (i.e., Manora Channel, Korangi creek, Phitti creek, and Khuddi creek). BIs 1-4 are not significantly affected by erosional effects of ocean factors due to being located behind the Manora channel. However, significant translocation was observed for Karachi's outlying twin drumstick-shaped BI-5 (i.e., Buddo and Bundal Islands), located near Korangi-Phitti Creeks (Figure 4). BI-5 migrated landward between 1989 and 2018 at 16.28 m/y, most likely under the erosional effects of wave action and sea level rise. We found very high erosion up to 24.4 m/y ongoing at the seaward side of the island, threatening the eastern parts of Karachi city, where one of the major deep seaports of the country, Port Qasim, is located (Figure 4). The presence of thick vegetation cover on the northern side of the BI-5 (Bundal Island) has protected it against the erosional effects of the port jetties for the passage of the cargo ships through the Phitti Creek to Port Bin Qasim (towards Korangi creek) and Port Karachi (towards Khizri creek) [49]. However, the absence of mangroves on the western and the southern parts have made them vulnerable to the erosional effects of intense wave action. It was found that Buddo (located adjacent to the Bundal Island's southern side-BI-5) is showing signs of accretion at an average rate of 0.49 m/y, particularly evident along its south-eastern lobe. This is most likely due to the accumulation of the eroded material coming from the southern side of Bundal Island.

Transect Level Assessment of Erosion for Local Changes
BIs 1, 2, and 6 experience low (<1 m/y), moderate (<2 m/y), and high erosion (>2 m/y) respectively, at all the transects, hence found to be depleting both the ocean-facing side and land/bay sides. BIs 5, 7, 9, and 11 experience high erosion at more than 80% of the transects on both the ocean side and land side (Table 3). BIs 8, 10, 17, 18, and 19 experience high erosion at more than 60% of the transects. Transects of BIs 3, 13, 14, 15, on the other hand, experience low to moderate erosion (up to 2 m/y). BIs 11-13 and 16-19 are experiencing erosion at a rate of > 2 m/y, except BI-14, which is located on the west side of the IR outfall and BI-15 situated on the east side of the IR outfall ( Figure 4). This could be due to their presence near the IR outfall, which makes it likely to receive some amount of sediments from upstream. BIs 4, 12, and 16 were found to be accreting at more than 50% transects with the highest accretion (> 2 m/y) found for BI-16. The lowest mean erosion rate based on LRR was found in BI-4 (0.25 m/y) located in the Karachi sub-region, whereas the highest mean erosion rate was found in BI-17 (29.12 m/y) located in the administrative unit of Kharo-Chan on the east side of the IR, during the period 1989-2018. For a better comparison among different BIs, a spatial histogram plot of the LRR-based erosion/accretion statistics for each BI was plotted ( Figure 5). These plots help us assess the severity of erosion/accretion in more detail for each BI, providing a broader range of information to coastal planners and decision makers [50,51]. Further analysis of erosion and accretion at different transects facing the Arabian sea and the mainland showed signs of translocation of large sand bodies in different BIs (e.g., BIs 15 and 16). It was found that eroding transects locations predominately exist on the seaside of the majority of the islands, such as BIs 13-17 ( Figure 6). Several BIs have nearly the same number of accreting areas as eroding areas, such as BI-3 (47% eroding, 50% accreting, and 3% no change), BI-4 (56% eroding and 44% accreting), BI-12 (44% eroding and 56% accreting), and BIs 14-15 (51% eroding and 49% accreting) ( Table 3). Accretion progrades the coasts and hence is beneficial, yet it takes a long time to use the accreted zones. This could mean the displacements of the BIs. However, the number of the accreting transects exceeds the number of eroding transects in only a few BIs (i.e., BI-2 (47% eroding and 50% accreting), BI-12 (44% eroding and 56% accreting), and BI-16 (34% eroding and 66% accreting)). However, in the other BIs where erosion is more pronounced, it leads to the destruction of the islands' banks, which lets the seawater intrude, destroy the mangroves, and deplete the island, particularly near the open shoreline of undeveloped BIs, where tidal action is strong (Figure 7).
Among the 19 BIs, the majority of the sand bodies located in BIs 8-9 and 11-19 were found to be translocating towards the mainland (i.e., seaward erosion and landward accretion). Few BIs (i.e., BIs 3, 4, 12, 13, 15, and 18) exhibited mixed behaviors of landward and seaward movement, and only 1 out of the 19 BIs was found to be moving towards the ocean (i.e., seaward accretion and landward erosion; parts of BI-18) ( Table 3). Only BI-16 was found to be accreting on both sea and landside, hence expanding both ways. This presents an interesting research question for further evaluations of these translocations of the BIs and the processes behind this phenomenon. the seawater intrude, destroy the mangroves, and deplete the island, particularly near the open shoreline of undeveloped BIs, where tidal action is strong (Figure 7).    Among the 19 BIs, the majority of the sand bodies located in BIs 8-9 and 11-19 were found to be translocating towards the mainland (i.e., seaward erosion and landward accretion). Few BIs (i.e., BIs 3, 4, 12, 13, 15, and 18) exhibited mixed behaviors of landward and seaward movement, and only 1 out of the 19 BIs was found to be moving towards the ocean (i.e., seaward accretion and landward erosion; parts of BI-18) ( Table 3). Only BI-16 was found to be accreting on both sea and landside, hence expanding both ways. This presents an interesting research question for further evaluations of these translocations of the BIs and the processes behind this phenomenon.   Among the 19 BIs, the majority of the sand bodies located in BIs 8-9 and 11-19 were found to be translocating towards the mainland (i.e., seaward erosion and landward ac cretion). Few BIs (i.e., BIs 3,4,12,13,15,and 18) exhibited mixed behaviors of landward and seaward movement, and only 1 out of the 19 BIs was found to be moving towards the ocean (i.e., seaward accretion and landward erosion; parts of BI-18) ( Table 3). Only BI-16 was found to be accreting on both sea and landside, hence expanding both ways. This presents an interesting research question for further evaluations of these translocations of the BIs and the processes behind this phenomenon.

Interdecadal Translocation of BIs
The shoreline retreat rate varied notably during the years 1989-2018 (Table 4). Each of the BIs had a unique evolution between 1989 and 2018 that altered its shape and position.  Table 4. The mean BI erosion was found to be the highest during the first epoch (1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)13.03 m/y), which gradually reduced during the following period, 1999-2009 (7.76 m/y). However, the average erosion for all BIs of the study area was still very high during the latest period, 2009-2018 (3.81 m/y). In the study area, the uncertainties of the EPR for the periods of 1989-1999, 1999-2009, and 2009-2018 were ±0.62, ±0.62, and ±0.7 m/y, respectively.  (Figure 8), which is likely to be caused by the reduced sand supply from upstream and/or increasing sea level rise (SLR) rate (Figures 8 and 9) [52]. Erosion in BI-13 and BI-15 (located near IR outfall) gradually reduced during the study period, whereas BI-14 located on the west of IR outfall showed signs of expansion seaward.

Hydrodynamics of the Study Area
BIs are wave-and wind-built landforms created in wave-dominated and/or mixed energy environments, typically in micro/meso-tidal environments [57,58]. Results from the hydrodynamical analysis of tidal waves confirm that the tide in the area is mixed, predominantly semidiurnal, with approximately two unequal high and two low water levels occurring in a day. The mean offshore SWH in the west of the east coastline of Pakistan is > 0.6 m during the pre-summer non-flooding monsoon month of April, whereas SWH during storms can reach up to 5.5 m. The tidal wave propagates eastward produced by the westerly wind (Supplementary Figure S2). Storm surges can significantly increase water levels along with the BIs and the coast, potentially leading to severe flooding of unprotected areas. In the south-westerly monsoon season, the wave breaking wind direction causes west to east littoral drift (alongshore transport of sediments), which could be depositing the eroded materials on the south-eastern side of the BIs on the west of the IR (Supplementary Figures S1 and S2) [28], thus not only changing the shape and position of these BIs located in the west of the IR but also creating new sand bars to appear along with a reduced erosion rate compared to those islands located on the west of the IR ( Table  3). The hydrodynamics of an area also actively control the sediment budget of the coast, and hence, its stability. There has been a significant reduction in the sediments in the Indus Delta coastal plains ( Figure 10) due to the construction of dams, barrages, and channeling in the up-stream areas [45]. Limited/reduced sediment supply from the upstream distributaries coupled with rising relative sea levels and shoreline erosion will eventually lead to the separation of such natural coastal barrier islands from the mainland of Mississippi delta [19,59] and further narrowing of the deltaic islands in the case of non-episodic flooding and storm events. In the same way, the drastic reduction in the sediment discharge from the IR below the Kotri barrage has increased the erosional effects of waves and tides, resulting in seawater intrusion into coastal lands [28], and is expected to cause the erosion of the BIs in IDR-particularly on the landward sides.

Discussion
This study demonstrates the application of earth observation data by combining the well-established shoreline evaluation methods to evaluate the shoreline changes in the BIs of Pakistan. The study provides important information to coastal planners and developers to direct the efforts for the stability of coastal areas. Communities who maintain their ecological features, such as BIs and coastal marshes, are resilient against storm surgeinduced flooding and inundation [52][53][54][55]. In this context, the results from this study provide relevant references for risk-informed decision making in highly eroding BIs.
The effectiveness of coastal protection measures is closely linked to the decisions made by the concerned authorities and stakeholders (public/private), which could lead to coastal adaptation through the utilization of proper information on coastal dynamics (i.e., erosion results from this study) and policy implications [56]. The Government of Pakistan has planned two coastal cities (i.e., Island City and Zulfiqarabad) along the IDR to be built on BI-5 and land covering four administrative units of Keti Bandar (BIs 13,14), Kharo-Chhann (BIs 15-16), Shah Bandar (BIs 17-19), and Jati. Given the high erosion status of BIs 13 and 16-19 ( Figure 6), future depletion of these BIs can make the planned cities vulnerable to inundation, erosion hazards, and coastal storms. The ongoing high erosion should be considered in the current and proposed development project design in the pre-construction phase to ensure the reduction in the current erosion rate as well as in the post-construction phase to assess the project's behaviors and impacts. Hence, in view of regional coastal sustainability, we recommend the consideration of these depleting natural resources in current and proposed coastal development project designs in the pre-construction phase on a priority basis. The study is a practical and applied demonstration of earth observations in the context of coastal sustainability. The historical sea level at the Karachi tide-gauge station showed an accelerated trend during the period 1989-2018 (PSMSL) (Figure 9). This rise in sea levels is likely to increase the inundation, accelerating erosion and even BI loss in future in the low-lying deltaic areas of Sindh Province if no protection measures are taken. Therefore, the results from this study can be utilized in planning and design while taking the erosion into account to tackle these. In addition to climate change, human activities (such as land reclamation for coastal development) are the triggers of climatic hazards, such as coastal flooding, rising sea levels, and cyclones, in the low-lying deltaic regions (e.g., IDR) increasing the vulnerability of coastal ecosystems and landmarks. Furthermore, the depletion of the BIs in the study area might not only affect Pakistan but also have regional consequences due to their various services. During pre-monsoon (March-May) and post-monsoon (October-November) periods, tropical cyclones and storms usually develop in the Arabian sea. The prevailing monsoon winds push them along. The Arabian Sea is northwest of the Indian ocean and shares the coast with India, Pakistan, Oman, Iran, Siri-Lanka, Maldives, and Somalia.

Hydrodynamics of the Study Area
BIs are wave-and wind-built landforms created in wave-dominated and/or mixed energy environments, typically in micro/meso-tidal environments [57,58]. Results from the hydrodynamical analysis of tidal waves confirm that the tide in the area is mixed, predominantly semidiurnal, with approximately two unequal high and two low water levels occurring in a day. The mean offshore SWH in the west of the east coastline of Pakistan is > 0.6 m during the pre-summer non-flooding monsoon month of April, whereas SWH during storms can reach up to 5.5 m. The tidal wave propagates eastward produced by the westerly wind (Supplementary Figure S2). Storm surges can significantly increase water levels along with the BIs and the coast, potentially leading to severe flooding of unprotected areas. In the south-westerly monsoon season, the wave breaking wind direction causes west to east littoral drift (alongshore transport of sediments), which could be depositing the eroded materials on the south-eastern side of the BIs on the west of the IR (Supplementary Figures S1 and S2) [28], thus not only changing the shape and position of these BIs located in the west of the IR but also creating new sand bars to appear along with a reduced erosion rate compared to those islands located on the west of the IR ( Table 3). The hydrodynamics of an area also actively control the sediment budget of the coast, and hence, its stability. There has been a significant reduction in the sediments in the Indus Delta coastal plains ( Figure 10) due to the construction of dams, barrages, and channeling in the up-stream areas [45]. Limited/reduced sediment supply from the upstream distributaries coupled with rising relative sea levels and shoreline erosion will eventually lead to the separation of such natural coastal barrier islands from the mainland of Mississippi delta [19,59] and further narrowing of the deltaic islands in the case of non-episodic flooding and storm events. In the same way, the drastic reduction in the sediment discharge from the IR below the Kotri barrage has increased the erosional effects of waves and tides, resulting in seawater intrusion into coastal lands [28], and is expected to cause the erosion of the BIs in IDR-particularly on the landward sides.
Coastal adaptation measures through building-with-nature are gaining increasing attention across the globe recently [56]. The ecosystem-based measures, such as the active plantation of coastal forests on the back-barrier basin/backshore side of islands, can support the stability of these fragile features through capturing wash over sand within the island system and re-building a platform onto which the island can translocate [60]. For example, mangrove might have protected the northern side of the outlying BI-5 and slowed down erosion as compared to the southern side ( Figure 6). Sand nourishment through placing dredged material and sediments on the backshore of the island fronting mainland shore can naturally develop the island, thereby slowing down the land loss from these deltaic islands. These measures will reduce not only coastal erosion but also provide multiple benefits, such as biodiversity conservation and provisioning of economic, biophysical, and social services [7,54,55,61]. Coastal adaptation measures through building-with-nature are gaining increasing attention across the globe recently [56]. The ecosystem-based measures, such as the active plantation of coastal forests on the back-barrier basin/backshore side of islands, can support the stability of these fragile features through capturing wash over sand within the island system and re-building a platform onto which the island can translocate [60]. For example, mangrove might have protected the northern side of the outlying BI-5 and slowed down erosion as compared to the southern side ( Figure 6). Sand nourishment through placing dredged material and sediments on the backshore of the island fronting mainland shore can naturally develop the island, thereby slowing down the land loss from these deltaic islands. These measures will reduce not only coastal erosion but also provide multiple benefits, such as biodiversity conservation and provisioning of economic, biophysical, and social services [7,54,55,61].

Policy Implications for Coastal Management
The results from this study have important implications relevant to decision making and policy formulation for better coastal planning and designing conservation strategies to reduce erosion and protect the BIs in Pakistan. For example, the BIs experiencing high erosion (categorized in the highest erosion group, erosion > 2 m/y) must be prioritized for appropriate measures. They should be considered in the process of coastal planning and development projects [28]. Otherwise, defense measures against climate change remain unproductive. The localized transect-based information of erosion/accretion rate provided in this study, which is not presented by Waqas et al. [33], is vital for local actions and decision making in the context of the long-term sustainability of coastal areas in Pakistan. These localized findings on coastal erosion can be used to scale down the study area and combined with other datasets (e.g., population, infrastructure, and natural systems) for further detailed analysis, such as the risk/vulnerability assessment of coastal communities to coastal erosion, which could progressively contribute towards coastal sustainability in the wake of global environmental change [62,63].

Policy Implications for Coastal Management
The results from this study have important implications relevant to decision making and policy formulation for better coastal planning and designing conservation strategies to reduce erosion and protect the BIs in Pakistan. For example, the BIs experiencing high erosion (categorized in the highest erosion group, erosion > 2 m/y) must be prioritized for appropriate measures. They should be considered in the process of coastal planning and development projects [28]. Otherwise, defense measures against climate change remain unproductive. The localized transect-based information of erosion/accretion rate provided in this study, which is not presented by Waqas et al. [33], is vital for local actions and decision making in the context of the long-term sustainability of coastal areas in Pakistan. These localized findings on coastal erosion can be used to scale down the study area and combined with other datasets (e.g., population, infrastructure, and natural systems) for further detailed analysis, such as the risk/vulnerability assessment of coastal communities to coastal erosion, which could progressively contribute towards coastal sustainability in the wake of global environmental change [62,63].
BIs depletion, particularly on the low-lying deltaic region, can make the coastal natural systems vulnerable to the frequent flooding and erosive effects of the rising sea levels. The maps and spatial database comprise a current and valuable source of information on the BIs evolution in Sindh, including historical perimeters and their change rate statistics. This spatial database can serve as a solid foundation for BIs' future assessment and applications. Similarly, the approach presented here can be adopted by academics, urban planners, developers, conservationists, and government agencies for a variety of purposes, such as long-term monitoring of the coastal natural systems in the context of coastal sustainability. Furthermore, the results of the BIs morphological assessment presented in this study can be integrated with the mainland shoreline analysis to provide a more comprehensive evaluation of the dynamics of Pakistan's shoreline-left for future studies.

Limitations and Future Directions
The authors acknowledge the current limitations of this study. For instance, while the transect-based method works well in a small subset of the areas of interest, the approach is not appropriate for assessments along the long and highly sinuous shoreline. For this reason, the whole study area was subdivided into smaller subsets, one for each island in this study. Though this process increases the computation time, the results provide detailed local information, which is overlooked by other in-practice shoreline change assessment methods, such as an area-based approach [33]. Some other parameters, not addressed within this study, may affect BI shoreline changes, such as creek and tidal inlet modification as a result of low sediment supply and increased SLR. Previous research has concluded the modification (narrowing/widening) in the creeks of the IDR [27], which could be a possible cause of the BI disintegration and transgression [52] in the study area. The influence of creek modification on the evolution and sustainability of the BIs should also be investigated.
Similarly, if possible, higher resolution satellite observations could be used for similar assessments, which could enhance the reliability of the results. However, it could significantly increase the computation costs and may require more time and higher capacities/capabilities, which could make its adoption difficult in developing countries. Since changes in the BIs' morphology impact the flood hazard and long-term coastal management depends on the changing shoreline position, depletion of the barrier island and flood hazard can be jointly studied for comprehensive coastal hazard risk/vulnerability assessment. The current study provides a basis for such further evaluations.
The sea level is rising along the IDR region, which could be the reason behind the estimated accelerated erosion of the BIs fronting the IDR, particularly those adjacent to the eastern side of the IR outfall at Khobar Creek, reaching up to 29 m/y ( Figure 3). However, further in-depth evaluations in this regard are a pre-requisite to improve our understanding regarding the correlations among sediment supply, sea-level rise, and the accelerated erosion of the BIs. This represents an interesting research question for future studies, such as the storm-induced risk assessment of these coastal barriers.

Conclusions
The BIs are of both regional as well as local significance due to the provisioning of multiple services, which makes them a useful natural asset. This study used satellite observations to evaluate the morphological changes in the BI chain of Pakistan from 1989 to 2018. The approach used here integrated GIS and remote sensing to assess the spatial-temporal erosion of 19 BIs in order to inform the coastal management (assessment, monitoring, and planning) in Pakistan. The rates of erosion were quantified at the transect level (n = 2530) and reported as the average statistics for each island. This spatial information on erosion is important for coastal planners and decision makers for taking appropriate actions and devising effective adaptation measures to mitigate the erosion process for coastal sustainability, particularly in the regions which are highly vulnerable to erosion. The results show that while the erosion increased in the majority of the islands on both sides (i.e., the side facing the Arabian sea and the mainland) during the period 1989-2018, the shorelines of some of the islands facing the mainland experienced some accretion.
It was found that~65% of BIs in the study area experienced high erosion,~15% experienced moderate, and~20% experienced low erosion during the period 1989-2018. This situation calls for consideration of the high erosion of BIs in the coastal planning and development plans in Pakistan. While the high erosion rates make the coastal communities vulnerable to erosion hazard, the depletion of the BIs in the study area can make the coastal population and infrastructure more susceptible to climate change-driven SLRinduced flooding and affect biodiversity due to the loss of natural systems. Therefore, it is imperative to assure the stability of these BIs in the context of coastal resilience. While the study acts as a roadmap to establish a coastal-change database for the monitoring of coastal environments in Pakistan, the approach presented here is equally applicable to other coastal regions-particularly data-scarce developing countries.
Supplementary Materials: The following are available online at https://www.mdpi.com/2077-131 2/9/3/295/s1, Figure S1: Wind climatology of the study area from January-December. Figure S2: Significant wave height and wave direction along the coast of Pakistan, driven by the monthly mean wind in representative months for northerly (January, February, November, and December), westerly (March and April), and south-westerly (May, June, July, August, September, and October) wind cases.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are openly available in the repositories, acknowledged and mentioned in the data source table within this article.