Spatial and Temporal Variability of Open-Ocean Barrier Islands along the Indus Delta Region

Barrier islands (BIs) have been designated as the first line of defense for coastal human assets against rising sea level. Global mean sea level may rise from 0.21 to 0.83 m by the end of 21st century as predicted by the Intergovernmental Panel on Climate Change (IPCC). Although the Indus Delta covers an area of 41,440 km2 surrounded by a chain of BIs, this may result in an encroachment area of 3750 km2 in Indus Delta with each 1 m rise of sea level. This study has used a long-term (1976 to 2017) satellite data record to study the development, movement and dynamics of BIs located along the Indus Delta. For this purpose, imagery from Landsat Multispectral Scanner (MSS), Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) sensors was used. From all these sensors, the Near Infrared (NIR) band (0.7–0.9 μm) was used for the delineation and extraction of the boundaries of 18 BIs. It was found that the area and magnitude of these BIs is so dynamic, and their movement is so great that changes in their positions and land areas have continuously been changing. Among these BIs, 38% were found to be vulnerable to oceanic factors, 37% were found to be partially vulnerable, 17% remained partially sustainable, and only 8% of these BIs sustained against the ocean controlling factors. The dramatic gain and loss in area of BIs is due to variant sediment budget transportation through number of floods in the Indus Delta and sea-level rise. Coastal protection and management along the Indus Delta should be adopted to defend against the erosive action of the ocean.


Introduction
Global warming, being caused by the increase in temperature and atmospheric CO 2 , acts as a catalyst in the melting of glaciers, expansion of ocean water, increase in sea surface temperature, rise in sea-level, and increase in tropical storms intensity [1][2][3][4][5][6].Accelerated sea level rise (SLR) threatens human settlements, and environmental and economic assets which have been tremendously developed in the coastal zones over the last five decades.SLR also shows alarming situation to the low lying sandy beaches and barrier islands (BIs) and intensifies erosion along the coastal areas.It has been predicted that, globally, up-to one meter rise in sea-level for the next hundred years can severely increase salinity of estuarine water, disturb coastal sediment budget supply, disrupt marine life, damage sub-surface and surface fresh water supplies, and damage agricultural and industrial areas [7][8][9].
Coastal topography may constantly change with numerous oceanic and physical processes including tidal flooding, SLR, and tsunami effect.This may cause sinking of low-land areas, and erosion of sediments which is a great concern for coastal researchers and investigators [10][11][12][13][14].Further sediments excavation, river modification and construction along coastal areas also play an important role in coastal areas including the BIs [15][16][17][18][19].
BIs are generally defined as the long offshore sand deposits which are 3 to 160 km in length and 1 to 3 km in width, and are separated but parallel to the main coastal area by the water inlets, bays and inland sea [20].BIs form a first line of defense for the main-land against high-energy tropical cyclones and rise in sea-level, which would otherwise directly transfer their full energy against the seashore.Provided coastal sand supplies are plentiful, and the gradient of sea-bed is gentle, the main-land is usually bordered by a chain of BIs.The water in the estuaries and bays is continuously exchanged with the water of the open shelf by flow through tidal inlets, water gaps separating BIs from one another [21].High wave energy is likely to close the tidal inlets while strong tidal flux keeps them open [22].A recent survey reveals that 20,783 km of coastline is occupied by 2149 BIs, worldwide [23].Therefore, rise in sea-level poses a considerable risk of 10% to the world's total population (7.6 billion) living in lowland coastal areas [24].Furthermore, tropical cyclones hitting BIs and the coastal systems will cause more devastation to the coastal-urban environment and infrastructure, and also disturb the geometries of the BIs [25][26][27].In order to track the changes in BIs shape, area, and their position, several methods have been adopted.For example, the topographic surveys, aerial photogrammetry, and GPS (Global Positioning System) surveys.Recently, the unmanned airial systems (UAS) have revolutionized the science and have proved their practice in remote sensing for short term change detection of coastal land masses, which include both autonomous and remotely piloted aircraft [28].Although these methods have high spatial and temporal resolution for coastal assessment of short-term changes as well as long-term trends, but they are labor intensive and expensive for non-funded research projects.Thus, the free availability of the medium spatial resolution satellite data (Landsat and Sentinel) in the last decade has enabled the use of multi-sensor, multi-spectral, and multi-temporal satellite imagery of extensive coastal areas to detect and monitor long-term trends of coastal land masses and coastline changes.

Background of the Study
Pakistan's shoreline extends more than 1000 km along the Arabian Sea which is divided into two coastal areas (i.e., coastal areas of Sindh and Makran).More than 10% of the population of Pakistan is living in the locality of coastal zone, about 20% of the coastal area of Pakistan is comparatively developed, and 40% of the country's industrial areas are located somewhere on or near the shore.The coast of Karachi city is about 70 km long surrounded by a chain of BIs, oriented NW-SE.Currently, 6.8% of the total population of Pakistan (15 million people) is living in the vicinity of the coastal city Karachi of the Indus Delta Region (IDR) which makes it the fifth largest coastal city [29] where a rate of SLR has increased from 1.1 mm/year [30,31] to 1.8 mm/year.IDR is considered as economic hub for Pakistan because Port Qasim and Karachi Port handle about 90% of all external trades.The IDR is stretched more than 200 km and always affected by high southwestern summer monsoon (May-September) winds having average speed of 15 m/s but during the northeastern winter monsoon winds below with an average speed of 5 m/s [32].Wave measuring at 20 m water depth offshore of Karachi city show a mean significant wave height of 1.8 m during SW monsoon with a mean wave period of 9 s and during the NE monsoon winds, significant wave height is about 1.2 m with a wave period of 6.5 s [33].The geomorphological changes have occurred at the center of Karachi coast along Bundal and Buddo Islands [34].IDR is experiencing geomorphological changes in all areas along its major and minor creeks in the form of erosion processes and low deposition rate of sediment budget from Indus River during post-damming era [35].IDR has been considered as one of the most dynamic deltaic systems which consists of 16 major tidal inlets having 17,000 km 2 area with an active tidal flat area of 10,000 km 2 , significant number of BIs [36] and hosts world's largest arid mangrove forest [37].
The IDR is ranked as the fourth highest delta in receiving average high wave energy as compared with other famous deltas around the world (e.g., Nile, Mississippi Niger, Ganges, and Ebro deltas) [38].The soil of IDR is less resistant to ocean hydro dynamics so the deltaic shoreline recede at average rate of 50 m/y along the central delta coast [35].SLR followed by sea-water intrusion is claiming more land and changing the geomorphology of the low-lying land masses in the Indus Delta [32,39,40].A SLR of about 1 meter is expected to sink or sea-encroach an area of 3750 km 2 in the Indus Delta [30].The global mean sea-level rose at a mean rate of 1.7 mm/year from 1900 to 2010 and at a rate 3.2 mm/year between 1993 and 2010 [41].Moreover, Intergovernmental Panel on Climate Change (IPCC) f claimed that globally, sea-level will rise from 0.18 m to 0.59 m by the end of 21st century [42] but in the current scenario it will rise to 0.21 m to 0.8 m [43].Therefore, if the sea-level rose at considerable rate, these narrow and low-lying landforms (i.e., BIs) would be extremely vulnerable so, they either migrate landward; experience erosion or accretion in its land mass in addition to other oceanic factors, but its overall response is unpredictable.
Several studies have investigated the behavior of BIs under different climatic and geographic setting, i.e., at global scale, Stutz and Pilkey (2011) [23] have studied the influence of sediments depositional settings, ocean climatic system and wave-tide regime on open-ocean BIs distribution and morphology.They revealed that 20,783 km of total worldwide coastal shoreline is occupied by 2149 open-ocean BIs.While at local scales, Madurapperuma et al. (2017) [44] used kite aerial photographs for mapping shoreline vulnerabilities at the Oluvil Harbor in Ampara, Sri Lanka and found that a 2 m rise in sea level inundated 90% area of the harbor.Kundu et al. (2014) [45] carried out shoreline mapping of the Sagar Island in West Bengal, India for the period of 1951 to 2011 using geospatial techniques to estimate morphological change and its future prediction.In a study of nearshore and foreshore influence on over-wash of a BI proposed by Matias et al. (2014) [46], identified and described longshore differences in storm impacts along a BI and evaluated the role of sub-aerial and submerged morphological variations in over-wash events on the western segment of Barreta Island, which is part of the Ria Formosa BI chain, in southern Portugal.Morton (2008) [47] has investigated that the Mississippi-Alabama BIs in the north-central Gulf of Mexico are undergoing rapid systematic land loss and translocation because of disruption of the sand-transport system.Moore et al. (2007) [48] has carried out sensitivity analyses on the outer banks of North Carolina and recommended that if sea-level rose by 0.9 m by the end of 21st century, the outer banks may translocate up to 2.5 times more quickly than at the present rate.

Purpose of the Study
Land reclamation has been a major anthropogenic factor affecting the Pakistan's coast and BIs and there has not been any systematic study for the evaluation.Therefore, this study aims at using the satellite remote sensing as a mean for (i) assessing the development, movement and change (morphology) in position of the BI chain along the Indus Delta from 1976 to 2017 and (ii) assessment of the rate of systematic land-loss due to physical processes.

Study Area
IDR is located between the Indian border along 'Sir Creek' on the east to the 'Hub River' on the west having length of 320 km and forming IDR (Figure 1).IDR is located at the head of the Arabian Sea, between Korangi Creek and the Rann of Kutch [30,40].IDR has 22 tidal creeks which supply the sediment budgets to the BI chain.Main creeks of IDR include Phitti, Waddi Khuddi, Dabbo, Hajambro, Wari, Khar, Keti Bunder, Chann, and Khobar connecting Indus River to Indian Ocean [8].
IDR has 2.6% of the world BI with total length of 567 km [23] (Figure 1).Indus Delta is mixed wave dominated delta [23] where continuous accretion and erosion due to waves and tides up to 3 m in height modify the BIs.Bundal Island is a triangular shaped island located at the eastern side of Karachi harbor and at the intersection of the three major creeks, Gizri creek, Korangi creek, and Phitti creek.Its northern side is wide, and the southern side is narrow.Bundal Island has great social-economic importance in the study area because it protects Phitti creek, a navigational channel of Port Bin Qasim, and the human settlements from the high energy wave.

Satellite Data
This study has used the archived imagery of Landsat-2 (L2) Multispectral Scanner (MSS), Landsat-5 (L5) Thematic Mapper (TM), Landsat-7 (L7) Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 (L8) Operational Land Imager (OLI).L2, L5, L7, and L8 satellites were launched during January 1975, March 1984, April 1999, and February 2013, respectively.L2 and L5 were decommissioned during July 1983 and January 2013, respectively.The MSS sensor had spatial resolution of 80 m with a revisit frequency of 18 days and acquired the imagery on the Worldwide Reference System-1 (WRS-1), while the TM, ETM+ and OLI sensors have spatial and temporal resolutions of 30 m and 16 days, respectively on Worldwide Reference System-2 (WRS-2).All the images acquired from ETM+ after 31 May 2003 suffer from the Scan Line Corrector (SLC) error which causes 22% loss of the data [49].Therefore, the SLC error was corrected using the "Fill nodata tool" under the raster tools available in QGIS 2.8.8 software.
Overall, this study has used 10 Landsat Collection 1 (C1) images from 1976 to 2017 acquired through L2 MSS, L5 TM, L7 ETM+, and L8 OLI sensors (Table A1).The Level-1 (radiometrically calibrated and orthorectified) MSS image were obtained from the United States Geological Survey (USGS) Earth Explorer website (http://earthexplorer.usgs.gov/)while the atmospherically corrected Level-2 (surface reflectance) images of TM, ETM+, and OLI were obtained from the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On Demand Interface (https://espa.cr.usgs.gov/).Images from different months of the years enabled the investigation of the morphological changes of the BIs over the study area according to the same oceanic and cloud free condition.

Satellite Data
This study has used the archived imagery of Landsat-2 (L2) Multispectral Scanner (MSS), Landsat-5 (L5) Thematic Mapper (TM), Landsat-7 (L7) Enhanced Thematic Mapper Plus (ETM+), and Landsat-8 (L8) Operational Land Imager (OLI).L2, L5, L7, and L8 satellites were launched during January 1975, March 1984, April 1999, and February 2013, respectively.L2 and L5 were decommissioned during July 1983 and January 2013, respectively.The MSS sensor had spatial resolution of 80 m with a revisit frequency of 18 days and acquired the imagery on the Worldwide Reference System-1 (WRS-1), while the TM, ETM+ and OLI sensors have spatial and temporal resolutions of 30 m and 16 days, respectively on Worldwide Reference System-2 (WRS-2).All the images acquired from ETM+ after 31 May 2003 suffer from the Scan Line Corrector (SLC) error which causes 22% loss of the data [49].Therefore, the SLC error was corrected using the "Fill nodata tool" under the raster tools available in QGIS 2.8.8 software.
Overall, this study has used 10 Landsat Collection 1 (C1) images from 1976 to 2017 acquired through L2 MSS, L5 TM, L7 ETM+, and L8 OLI sensors (Table A1).The Level-1 (radiometrically calibrated and orthorectified) MSS image were obtained from the United States Geological Survey (USGS) Earth Explorer website (http://earthexplorer.usgs.gov/)while the atmospherically corrected Level-2 (surface reflectance) images of TM, ETM+, and OLI were obtained from the USGS Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) On Demand Interface (https://espa.cr.usgs.gov/).Images from different months of the years enabled the investigation of the morphological changes of the BIs over the study area according to the same oceanic and cloud free condition.

Tidal Data
Tide height was used in the satellite image selection process for the identification and extraction of geometry of the BIs.Ocean tide height along the coast of Sindh, Pakistan is monitored using a fixed tide gauge located at Port Muhammad Bin Qasim (24.7833 • N, 67.3500 • E) by National Institute of Oceanography (NIO) Karachi, Pakistan (Figure 1b).The maximum and minimum tide heights recorded at this station are 3.8 m and little less than 0.1 m between 1976 to 2017.

Satellite Data Selection Criteria
Due to seasonal variations in the ocean controlling factors (wave, tides, sea level rise, storms and sediment budget) and suitable atmospheric conditions, selection of satellite images at the same time of the year is important for studying the morphology of BIs [50].Therefore, only those Landsat (MSS, TM, ETM+, and OLI) images were selected which had 0% cloud cover over the study area duly verified from the metadata file (.MTL) which is provided along each image.The objective for taking only those scenes with 0% cloud cover is due to the limitation of the optical imagery.This criterion can be relaxed to 10% to 20% provided the studied area (or the specific barrier island) is not affected.This limitation can be overcome by using the Synthetic Aperture Radar (SAR) imagery and authors do not exclude this possibility.
Furthermore, the images from non-flooding months (December to April) were selected because tropical storms affect the coastal areas of Pakistan during May to June and September to November.Only those images were selected which had a tide height of <1 m within a time window of ±1 h before or after the image acquisition time.This image selection criteria (Figure 2) resulted in ten images including one image from L2 MSS, two images from L5 TM, six images from L7 ETM+, and one image from L8 OLI sensor from 1976 to 2017.The selected images were obtained as Level-2 surface reflectance products for the Landsat TM, ETM+ and OLI sensors, due to the unavailability of the surface reflectance product Landsat MSS, surface reflectance was estimated in house using the 6S atmospheric correction method [51] (supplementary materials).

Tidal Data
Tide height was used in the satellite image selection process for the identification and extraction of geometry of the BIs.Ocean tide height along the coast of Sindh, Pakistan is monitored using a fixed tide gauge located at Port Muhammad Bin Qasim (24.7833°N, 67.3500°E) by National Institute of Oceanography (NIO) Karachi, Pakistan (Figure 1b).The maximum and minimum tide heights recorded at this station are 3.8 m and little less than 0.1 m between 1976 to 2017.

Satellite Data Selection Criteria
Due to seasonal variations in the ocean controlling factors (wave, tides, sea level rise, storms and sediment budget) and suitable atmospheric conditions, selection of satellite images at the same time of the year is important for studying the morphology of BIs [50].Therefore, only those Landsat (MSS, TM, ETM+, and OLI) images were selected which had 0% cloud cover over the study area duly verified from the metadata file (.MTL) which is provided along each image.The objective for taking only those scenes with 0% cloud cover is due to the limitation of the optical imagery.This criterion can be relaxed to 10% to 20% provided the studied area (or the specific barrier island) is not affected.This limitation can be overcome by using the Synthetic Aperture Radar (SAR) imagery and authors do not exclude this possibility.
Furthermore, the images from non-flooding months (December to April) were selected because tropical storms affect the coastal areas of Pakistan during May to June and September to November.Only those images were selected which had a tide height of <1 m within a time window of ±1 h before or after the image acquisition time.This image selection criteria (Figure 2) resulted in ten images including one image from L2 MSS, two images from L5 TM, six images from L7 ETM+, and one image from L8 OLI sensor from 1976 to 2017.The selected images were obtained as Level-2 surface reflectance products for the Landsat TM, ETM+ and OLI sensors, due to the unavailability of the surface reflectance product Landsat MSS, surface reflectance was estimated in house using the 6S atmospheric correction method [51] (supplementary materials).

Barrier Island Identification and Extraction
A BI can be identified based on its association with different features i.e., a BI always has a back-barrier lagoon, a tidal inlet, it faces the shore and it is distant from the land [52].Taking advantage of the higher reflectivity of land features in the Near Infrared (NIR) wavelength, the NIR band of L2 MSS (0.7-0.8 µm), L5 TM (0.76-0.90 µm), L7 ETM+ (0.77-0.90 µm), and L8 OLI (0.85-0.88 µm) sensors was used for developing water/land mask through Iso-cluster image classification method [53].
The classified images (water/land mask) were then visually examined with their corresponding true color images by displaying at a scale of 1:75,000 for visually identifying the boundary of BIs.Shorelines of the wet marshy sand beaches around the perimeter of BIs are less reliable and vary more with tide level than sandy shorelines of island [54].So, only elongated sand deposits which fulfill the criteria for BI identification [52] were extracted.The boundaries of the classified BIs were then converted to polygon features by using the raster to vector conversion tool available in QGIS 2.8.8 software.Where necessary, minor edits were made to the boundaries of BIs to extract the more accurate boundary.Small polygons found near the main perimeter of an island were merged to that specific island polygon.Furthermore, polygon areas were calculated and added to the polygon attribute tables.Through this process 18 BIs (Figure 1) were delineated for the whole IDR.

Pixel by Pixel Frequency of Barrier Island
After finalizing the boundary of each BI from each of the ten selected images (Table A1), raster overlay operation was applied to determine the pixel by pixel frequency of each existing valid pixel of island body.For this purpose, polygon feature of each BI, for each year was converted to a binary raster (with a native grid size of 30 × 30 m) using the vector to raster conversion tool available in QGIS 2.8.8 software where a filled value of 1 was assigned to each pixel of BI, and 0 for each pixel representing any feature(s) other than BI.
The boundary of each BI was different in each of the 10 images but for raster overlay operation, one need to have the same spatial extent for each BI.Therefore, a fixed spatial extent (suitable for all years for each BI) was created by drawing a rectangular polygon to the maximum extent on the east, west, north and south directions of that specific BI (Figure 1-red rectangles).At this point when spatial extent of each BI was same, raster overlay by addition operation was applied, which calculated the frequency value of 1 to 10 for each pixel of the BI.The pixel which shows frequency value of 10, indicates that the island body exists at that pixel (location) in each of the ten rasters (images) while descending number of frequency values (9 to 1) shows that existence of valid island pixels for each image is decreasing.

Boundary Delineation of Barrier Islands
The variant changes observed in the boundary of BIs during study period indicate that generally, there is a retreat of the boundary, and shapes of the BIs have evolved (Figure 3).In Figure 3, the color variations represent the shape of a BI in a specific year.For BIs 1 to 13 the 1976 boundary (filled grey colored polygons) was used as a reference to gauge the movement of the BIs, while it should be noted that for BIs 14 to 18 there is no boundary for the year 1976 which was due to the unavailability of the L2 MSS image satisfying the image selection criteria (Figure 2), hence 1990 boundary (filled grey colored polygons with red outline) was used as a reference.It is evident that, except BIs 1 to 4, 6, and 13, the shape of the BIs has been evolving significantly which can be due to the climatic and oceanic factors.It has been observed that some of the BIs are shrinking while others are expanding in both directions-i.e., oceanside and landside-against the accelerated SLR.
Among all the 18 BIs, BIs 5 and 9 (Figure 3) have lost significant areas and moved towards the land.BI5 also known as Bundal Island (local name) is the largest island located near the most western part of the Indus Delta.In personal communication with the NIO officials we came to know that the island has faced erosion from tidal currents passing through small prolonged tidal inlets on the eastern side and due to the passage of the cargo ships from the southern part of the island at the opening of the Phitti Creek, which produces constant backward and forward ocean currents.The island is abundant in mangroves on the northern side which helps it to retain its northern side.The boundary of the island is shrinking slowly from the eastern and southern sides while expanding from the western side as it is protected from any direct erosive factor (Figure 3).An overall summary of BI wise morphological changes from 1976 to 2017 has been presented in Table 1.

Sustainability of Barrier Islands
The vulnerability or sustainability of BI over the study period against the climatic and oceanic factors such as storm erosion, reductions in sediment, longshore drift and SLR shows that the BI are mostly vulnerable to greater extent (Figure 4).The pixel frequency values from 1 to 10 for each BI revels about the existence of the island body covering 10 years of study period at that pixel which is important to assess the sustainability of the BI chain of the Indus Delta.Increasing pixel value from 1 to 10 shows more stable part of the island, such that 10 indicates an area present throughout the study period whereas 1 indicates a part of the island that appears only once in study period (Figure 4).

Sustainability of Barrier Islands
The vulnerability or sustainability of BI over the study period against the climatic and oceanic factors such as storm erosion, reductions in sediment, longshore drift and SLR shows that the BI are mostly vulnerable to greater extent (Figure 4).The pixel frequency values from 1 to 10 for each BI revels about the existence of the island body covering 10 years of study period at that pixel which is important to assess the sustainability of the BI chain of the Indus Delta.Increasing pixel value from 1 to 10 shows more stable part of the island, such that 10 indicates an area present throughout the study period whereas 1 indicates a part of the island that appears only once in study period (Figure 4).

Change in Area of Barrier Islands
The box and whisker plot (Figure 5) show the percent change in area from 1976 to 2017 (for each BI) on the vertical axis, while the horizontal axis shows the BIs from 1 to 18.The red line shows the zero line of the data distribution.The value above the zero line shows the gain in area of a specific BI while the values below show the loss in area.It can be observed that Island 4 sustained, so we can say that the size of the box varies as there is variant change in area of the BI from 1976 to 2017.The Island 11 (Figure 5) shows the most vigorous change above the median value means the variance in data is higher in the 50th to 75th percentile of the data.

Change in Area of Barrier Islands
The box and whisker plot (Figure 5) show the percent change in area from 1976 to 2017 (for each BI) on the vertical axis, while the horizontal axis shows the BIs from 1 to 18.The red line shows the zero line of the data distribution.The value above the zero line shows the gain in area of a specific BI while the values below show the loss in area.It can be observed that Island 4 sustained, so we can say that the size of the box varies as there is variant change in area of the BI from 1976 to 2017.The Island 11 (Figure 5) shows the most vigorous change above the median value means the variance in data is higher in the 50th to 75th percentile of the data.

Translocation of Barrier Islands
The coastal erosion over centuries is a result of natural processes and sea level change.The rate of erosion seems to have increased at some points along the BIs of the Indus Delta.Translocation of BIs was estimated based on the movement of that specific BI relative to the reference shape.The reference shape (Figure 6, black colored polygon) for each BI was the shape delineated from the L2 MSS (from 1976) and L5 TM (from 1990) images for BIs 1 to 13 and 14 to 18, respectively.In Figure 6, the shades of blue and red colors represent less sustainable and highly sustainable pixel of the BIs.Due to the existence of BIs 1 to 4 on the back side of the Manora Beach, they were not significantly migrated or translocated (Figure 6).A significant translocation was observed for Bundal Island (BI 5) which migrated landward approximately 1.5 km (Table 2) from the ocean exposing side of the island between 1976 and 2017 with a rate of 38 m per year due to the action of ocean dynamical factors.BIs 6, 7, 8, 15, 16, and 17 have gained a significant area since 1976 and overall, they have started their

Translocation of Barrier Islands
The coastal erosion over centuries is a result of natural processes and sea level change.The rate of erosion seems to have increased at some points along the BIs of the Indus Delta.Translocation of BIs was estimated based on the movement of that specific BI relative to the reference shape.The reference shape (Figure 6, black colored polygon) for each BI was the shape delineated from the L2 MSS (from 1976) and L5 TM (from 1990) images for BIs 1 to 13 and 14 to 18, respectively.In Figure 6, the shades of blue and red colors represent less sustainable and highly sustainable pixel of the BIs.Due to the existence of BIs 1 to 4 on the back side of the Manora Beach, they were not significantly migrated or translocated (Figure 6).A significant translocation was observed for Bundal Island (BI 5) which migrated landward approximately 1.5 km (Table 2) from the ocean exposing side of the island between 1976 and 2017 with a rate of 38 m per year due to the action of ocean dynamical factors.BIs 6, 7, 8, 15, 16, and 17 have gained a significant area since 1976 and overall, they have started their development towards the land.While BIs 5, 9, 11, and 12 were found to be highly vulnerable to the oceanic conditions and have lost a significant area.
development towards the land.While BIs 5, 9, 11, and 12 were found to be highly vulnerable to the oceanic conditions and have lost a significant area.A BI wise status of its vulnerability or sustainability has been devised based on the frequency of the pixels-i.e., a BI whose pixel frequency count of 1 is high indicates that specific BI is 'vulnerable' (V) to oceanic factors, frequency pixel counts of 2 to 4, 5 to 8, and 9 to 10 indicate that specific BI is 'partially vulnerable' (PV), 'partially sustainable'(PS), 'and sustainable' (S) during the study period respectively (Table 2).The movement either it is landward, oceanward or in both directions, has also been reported in Table 2 along with the fate of BI whether it is shrinking or expanding.It can be seen that only islands 1 to 4 have more than 40% of their areas which are sustainable, islands 5 to 14 have areas majority of the areas between partially vulnerable to partially sustainable while islands 15 to 18 have more than 50% of their areas which are vulnerable to oceanic factors.Among the 18 BI of the Indus Delta BI chain, majority (11 out of 18) have been found to be moving towards land, a smaller amount (6 out of 18) of BIs has been found which exhibited the mixed behavior of movement (i.e., some parts are moving towards land and some towards ocean) and only 1 out of 18 are found to be moving towards ocean (Table 2).This is an alarming situation and indicates that the sea-level is rising and encroaching the land areas of Indus Delta BI chain.Furthermore, a notable remark is that 10 out of 18 BIs have been found which are expanding while 8 out of 18 BIs are shrinking (Table 2).' represent shrinkage, expansion, landward movement, oceanward movement, and movement in both directions, respectively.

Discussion
This study attempts to study the historical variability of the barrier islands (BIs) located along the Indus Delta Region by employing the Landsat satellite imagery from 1976 to 2017.Overall, 18 BIs were delineated, and their morphological behavior was studied.The spatial resolution is very important in detection and monitoring of the trends of costal land masses and coastline changes.Therefore, taking the advantage of the medium spatial resolution and free availability of the Landsat imagery, this study has used the Landsat TM, ETM+, and OLI imagery at 30 m spatial resolution (except for MSS sensor which has a spatial resolution of 60 m).In recent years, image acquisition systems with high spatial resolution such as spaceborne sensors and unmanned aircraft systems (UAS) or drones have emerged as an important platform for high spatial resolution data collection [28].The evaluation of detection accuracy is always the focus of discussion.Therefore, results of this study have been verified by authors themselves during frequent field visits with NIO officials and personal communications with the residents, where applicable.
It was found that BIs 5 and 9 faced significant erosion and moved towards the land.BI 5 also known as Bundal Island (local name) is the largest island located near the most western part of the Indus Delta.This island is of high socio-economic importance for the region as it has human settlements from past three decades.Significant net loss of 61% in the area of the Bundal Island and translocation of about 1.56 km of its boundary landward with a rate of 38 m/year shows vulnerability of the socioeconomic coastal environment associated to the island.Historical morphological changes show maximum of 36% of the Bundle Island was vulnerable to oceanic controlling factors with only 6% of the island sustained throughout the study period.Tidal currents contributed to the erosive and accretion action along the Bundal Island.The tidal inlets with narrow opening at the mouth and a large inside the Island let the high tide water enters makes island marshy and erosive.Among the other BIS, BIs 6, 7, 8, 15, 16, and 17 had gained a significant area since 1976 and overall, they had started their development towards the land.While BIs 5, 9, 11, and 12 were found to be highly vulnerable to the oceanic conditions and have lost a significant area.
Results suggested that the BIs 1 to 4, 5 to 14, and 15 to 18, are sustainable, partially sustainable and vulnerable to the oceanic factors, respectively.Majority of the Indus Delta BIs have been found to be moving towards land which alerts the local authorities and indicates that the sea-level is rising and encroaching the land areas of the BIs.Climate change and anthropogenic activities (land reclamation and modification) can trigger extreme climatic events such as flooding, tropical cyclone and SLR along the deltaic regions (e.g., Indus Delta Region) leading to an increased vulnerability of coastal landmarks.Cyclones and storms usually develop in the Arabian sea during pre-monsoon (March-May) and post-monsoon (October-November) seasons with favorable month as October and pushed by the prevailing monsoon winds along their direction.Arabian Sea is in northwest of the Indian ocean and share coast among India, Pakistan, Oman, Iran, Siri-Lanka, Maldives, and Somalia.

Conclusions
An archive of Landsat imagery from 1976 to 2017, has been used in this study to explore the morphological changes in the barrier island of the Indus Delta Region (IDR).Unlike global scale assessment of the barrier islands, a regional study allows in depth spatial assessment of the island.In general, it has been found that, accretion is more significant in the most parts of the islands, but erosion has also remarkably increased in most of the islands.The relative rates of change of islands morphology and migration are recorded in the IDR barrier islands chain.The natural future trends of change for the IDR barrier islands will be continued in rapid land mass change, change in perimeter, and migration as a result of accelerated sea level rise, frequent intense storms, and sand budget supply.As a result of global warming which causes the sea level rise increases beyond the expected rate will also likely to increase the land loss of the barrier islands regionally and worldwide.The total percent area of barrier islands relatively changed from 1976 to 2017 was found −16.73%.Overall, 44.4% of the total barrier island are shrunk while the rest of them are expanded due to heavy sediments budget supply.The translocation of the barrier island found about of 66.1% of the total islands moved landward, 5.6% along seaward and the rest of the islands moved along both sides.A major part of these barrier islands chain is vulnerable which is about 38.3% and 36.7% is partially vulnerable while 17.4% is partially sustainable and only 7.6% sustained in the study area against the ocean controlling factors.
Most of the barrier islands along the IDR are still undisturbed by the human intervention.The existence of the barrier islands at a place is mostly determined by the sand supply and sea level change.Protection against the sea level rise is an especially important issue which can damage barrier beaches and coastal system.Coastal protection and coastal management along the IDR should be adopted to defend against the erosive action of the ocean.In future, increase in land loss can be mitigated by placing dredged material on the backshore and shore facing side of the island so the island experiences nourishment and rebuilding.

18 Figure 1 .
Figure 1.(a) Location of study area (filled red rectangle) and worldwide Barrier Island (BI) chains (black dots) derived from [23] and (b) barrier islands located within Indus Delta Region (IDR).

Figure 1 .
Figure 1.(a) Location of study area (filled red rectangle) and worldwide Barrier Island (BI) chains (black dots) derived from [23] and (b) barrier islands located within Indus Delta Region (IDR).

Figure 2 .
Figure 2. Flowchart of the methodology for the morphological change detection in the IDR barrier islands.

18 Figure 3 .
Figure 3. Barrier islands shapes as delineated from 1976 to 2017 from Landsat images.

Figure 3 .
Figure 3. Barrier islands shapes as delineated from 1976 to 2017 from Landsat images.

Figure 4 .
Figure 4. Frequency pixel count of Barrier Island.The inset views for BI 1&2, BI 3&4, and BI 6 show the frequency of the data between 0 and 1 × 10 3 .

Figure 4 .
Figure 4. Frequency pixel count of Barrier Island.The inset views for BI 1&2, BI 3&4, and BI 6 show the frequency of the data between 0 and 1 × 10 3 .

Figure 5 .
Figure 5. Barrier island percent area change from 1976 to 2017.

Figure 5 .
Figure 5. Barrier island percent area change from 1976 to 2017.

Figure 6 .
Figure 6.Pixel by pixel frequency.The black colored polygons show the extent of a barrier island in 1976 (for barrier islands 1 to 14) and 1990 (for barrier islands 14 to 18).

1 )
Net translocation and translocation rate for barrier island 1 to 13 and barrier islands 14 to 18 are with reference to the periods 1976 to 2017 and 1990 to 2017, respectively.(2) V, PV, PS, and S represent vulnerable, partially vulnerability, partially sustainability and sustainable status of barrier islands, respectively and the signs '−', '+', '↑', '↓' and ' Remote Sens. 2018, 10, x FOR PEER REVIEW 1 of 1 ↕

Table 1 .
Morphological changes in barrier islands with respect to area and perimeter from 1976 to 2017 in different temporal periods.Net area change, and net perimeter change for barrier island 1 to 13 and barrier islands 14 to 18 are with reference to the year 1976 and 1990, respectively. Note:

Table 2 .
Net translocation, translocation rate, vulnerability/sustainability status, movement and the fate of barrier islands.