Three Decades of Coastal Changes in Sindh, Pakistan (1989–2018): A Geospatial Assessment

: Coastal erosion endangers millions living near-shore and puts coastal infrastructure at risk, particularly in low-lying deltaic coasts of developing nations. This study focuses on morphological changes along the ~320-km-long Sindh coastline of Pakistan over past three decades. In this study, the Landsat images from 1989 to 2018 at an interval of 10 years are used to analyze the state of coastline erosion. For this purpose, well-known statistical approaches such as end point rate (EPR), least median of squares (LMS), and linear regression rate (LRR) are used to calculate the rates of coastline change. We analyze the erosion trend along with the underlying controlling variables of coastal change. Results show that most areas along the coastline have experienced noteworthy erosion during the study period. It is found that Karachi coastline experienced 2.43 ± 0.45 m / yr of erosion and 8.34 ± 0.45 m / yr of accretion, while erosion on the western and eastern sides of Indus River reached 12.5 ± 0.55 and 19.96 ± 0.65 m / yr on average, respectively. Coastal erosion is widespread along the entire coastline. However, the rate of erosion varies across the study area with a general trend of erosion increasing from west to east in the Indus Delta region (IDR), and the highest average erosion rate is 27.46 m / yr. The interdecadal change during 1989–1999, 1999–2009 and 2009–2018 periods depicted an increasing linear trend ( R 2 = 0.78) from Karachi to Indus River (IR) East zone. The spatial trend from west to east is positively correlated with mean sea level rise, which has increased from 1.1 to 1.9 mm / year, and negatively correlated with topographic slope, which is found to be decreasing eastward along the coastline. The ﬁndings necessitate appropriate actions and have important implications to better manage coastal areas in Pakistan in the wake of global climate change.


Introduction
Coastal degradation endangers millions living near-shore and puts coastal infrastructure at risk. Impacts of coastal erosion are often sudden, making it a major coastal hazard [1,2]. The situation is more appalling in low-lying deltaic regions in developing nations (such as Pakistan, Philippines, and Bangladesh), which are less prepared to cope with the risks. Some long-term environmental processes such as sea level rise (SLR) [3] and increasing intensity of short-term events, such as storm surge and coastal flooding [4], manifest themselves in the form of coastal morphological change and intensify coastal erosion [5]. The trends and rates of morphological change of the highly dynamic land-water transition interface, known as coastline, represent the summaries of various coast-influencing processes such as SLR and storm surges [6]. Anthropogenic activities are another major factor accelerating coastal river's lower reaches, causing degradation of the sixth largest delta of the world. The abrupt cut-down in the sediment load and water discharge to the Arabian Sea has increased the impact of waves and tides and impeded the growth of the mangroves, resulting in seawater intrusion and accelerated coastal erosion. The IR drains into the Arabian Sea at Khobar Creek [25], bisecting the deltaic region into two zones, IR West and IR East. Following Ijaz et al. [26], we also divide the Sindh coastline into three zones: Karachi, IR West, and IR East (Figure 1). Each of the coastal zones contains a varying number of administrative units (Table 1). We will not include the barrier islands (BIs) in the analysis as they are very localized features along the coastline.
Remote Sens. 2020, 12, x FOR PEER REVIEW 3 of 24 down in the sediment load and water discharge to the Arabian Sea has increased the impact of waves and tides and impeded the growth of the mangroves, resulting in seawater intrusion and accelerated coastal erosion. The IR drains into the Arabian Sea at Khobar Creek [25], bisecting the deltaic region into two zones, IR West and IR East. Following Ijaz et al. [26], we also divide the Sindh coastline into three zones: Karachi, IR West, and IR East (Figure 1). Each of the coastal zones contains a varying number of administrative units (Table 1). We will not include the barrier islands (BIs) in the analysis as they are very localized features along the coastline.

Zone Administrative Units Karachi
Karachi West, Karachi South, and Karachi East

IR East
Kharo-Chhann and Shah Bandar, located between Ghara-Chan Creek and Sir Creek

Data Set for Coastline Delineation
The IDR creek system is under the influence of high tides [26], medium to high energy waves [32], and sea level rise of over 1.1 mm/yr [27,33]. The coastline delineated from satellite images varies with the tidal height when the satellite passes [34] (Figure 2). The National Institute of Oceanography (NIO), Pakistan, monitors the ocean tide height using a fixed tide gauge located at Port Muhammad Bin Qasim (24.7833°N, 67.3500°E). The tidal status at the time of image acquisition was based on

Zone
Administrative Units the tidal height when the satellite passes [34] ( Figure 2). The National Institute of Oceanography (NIO), Pakistan, monitors the ocean tide height using a fixed tide gauge located at Port Muhammad Bin Qasim (24.7833 • N, 67.3500 • E). The tidal status at the time of image acquisition was based on astronomical tidal height predictions from Tides 4 Fishing (https://tides4fishing.com/), which provides expected tidal height values and the times they are likely to occur. The data from the tide gauge show that the tidal height varied approximately between 0.20 to 2.80 m during the passes of Landsat TM (Thematic Mapper) and OLI (Operational Imager) over the past three decades (Table 2).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 24 astronomical tidal height predictions from Tides 4 Fishing (https://tides4fishing.com/), which provides expected tidal height values and the times they are likely to occur. The data from the tide gauge show that the tidal height varied approximately between 0.20 to 2.80 m during the passes of Landsat TM (Thematic Mapper) and OLI (Operational Imager) over the past three decades (Table 2). We selected only cloud-free and atmospherically corrected 30-m Landsat (TM: Thematic Mapper and OLI: Operational Image Processor) images acquired in non-flooding pre-monsoon season (Dec-April) to evaluate the geomorphological change of the coastline [27,35]. All the Landsat CDR (climate data records) images were pre-processed to L1T level (i.e., after standard terrain correction) with consistent geo-registration and within the prescribed tolerances (i.e., 12 m root-mean-square error) [11]. Due to the data availability and tidal level constraint (≤2.8 m), only some of the cloud-free images were suitable. The final data set contained 17 images from path 152, row 043. Four of these images were used for coastline delineation and change analysis. The data used, the software for processing the data, and related information are provided in Table 2.  We selected only cloud-free and atmospherically corrected 30-m Landsat (TM: Thematic Mapper and OLI: Operational Image Processor) images acquired in non-flooding pre-monsoon season (Dec-April) to evaluate the geomorphological change of the coastline [27,35]. All the Landsat CDR (climate data records) images were pre-processed to L1T level (i.e., after standard terrain correction) with consistent geo-registration and within the prescribed tolerances (i.e., 12 m root-mean-square error) [11]. Due to the data availability and tidal level constraint (≤2.8 m), only some of the cloud-free images were suitable. The final data set contained 17 images from path 152, row 043. Four of these images were used for coastline delineation and change analysis. The data used, the software for processing the data, and related information are provided in Table 2.

Coastline Delineation
Toure et al. [37] reviewed various approaches used for coastline delineation from remote sensing imagery using various image processing methods. These methods to determine the coastal changes include pixel frequency count, image histogram thresholding [38], edge filter algorithm and tasseled cap transformation [18], and image classification through iterative self-organizing data analysis (ISODATA) [39]. We adopted an automatic coastline delineation method that consists of the following steps [40]: (a) band selection based on spectral profile curve of the images; (b) selection of the optimal index for coastline delineation; and (c) coastline delineation based on histogram threshold. The comprehensive framework is shown in Figure 3.
Prior to coastline delineation, all the satellite images are georeferenced and projected to a common Universal Transverse Mercator (UTM) map system (42N Zone) using ENVI 5.3 (ENvironment for Visualizing Images). ENVI's X and Y coordinates and arbitrary profiles (transects) are used to generate terrain profiles perpendicular to the coastline and analyze the spectral profile curves of the final composite normalized images (i.e., image composite of all spectral bands) to identify the suitable signal bands for coastline delineation. Shortwave infrared Band 5 (1.55-1.75 µm) is well-suited to extract the land-water interface from Landsat TM [41] and OLI images [18], as it shows the highest gradient between the seawater and the land ( Figure S1). The modified normalized difference water index (MNDWI) was implemented to the images of 1989, 1999, 2009, and 2018, using the equation: where Green is the green band (band 2 for Landsat TM/ETM+ data and band 3 for Landsat OLI data) and MIR is the middle infrared band (band 5 for Landsat TM/ETM+ data and band 6 for Landsat OLI data). MNDWI is a modification of the normalized difference water index (NDWI) and has been widely used in detecting water surfaces. The object-based image analysis (OBIA) was finally carried out for automatic delineation of the land-water interface by applying the histogram threshold to the MNDWI images. Suitable thresholds were determined for each index image in a trial-and-error process. A threshold value of −0.5 was found to be stable and repeatable for coastline mapping based on the MNDWI.
include pixel frequency count, image histogram thresholding [38], edge filter algorithm and tasseled cap transformation [18], and image classification through iterative self-organizing data analysis (ISODATA) [39]. We adopted an automatic coastline delineation method that consists of the following steps [40]: (a) band selection based on spectral profile curve of the images; (b) selection of the optimal index for coastline delineation; and (c) coastline delineation based on histogram threshold. The comprehensive framework is shown in Figure 3. Prior to coastline delineation, all the satellite images are georeferenced and projected to a common Universal Transverse Mercator (UTM) map system (42N Zone) using ENVI 5.3 (ENvironment for Visualizing Images). ENVI's X and Y coordinates and arbitrary profiles (transects) are used to generate terrain profiles perpendicular to the coastline and analyze the spectral profile curves of the final composite normalized images (i.e., image composite of all spectral bands) to identify the suitable signal bands for coastline delineation. Shortwave infrared Band 5 (1.55-1.75 µm) is well-suited to extract the land-water interface from Landsat TM [41] and OLI images [18], as it shows the highest gradient between the seawater and the land ( Figure S1). The modified normalized difference water index (MNDWI) was implemented to the images of 1989, 1999, 2009, and 2018, using the equation: Each land-water classified image was compared with their corresponding Landsat image for visual assessment and to ensure the quality of the classification results. The continuous edge representing the coastline was delineated using an edge detection approach based on a morphological filter [40]. The binary data were then converted into vector format for further assessment. All the inland water bodies, such as lakes, that are not adjacent to the Arabian Sea were removed. It was noted that the study area has very dense dendritic networks of creeks extending from~180 km along the coast and 32 km landward. Using prior spatial knowledge, about 2.5 km inward passage was identified as a coastline.

Analysis of Coastline Change
Decadal coastline change was analyzed systematically using the Digital Shoreline Analysis System (DSAS) [18,42] developed by the USGS. The following computation was carried out to compute the change rate: (a) coastline delineation; (b) baseline generation; (c) transects generation; (d) computation of distances between the coastline and the baseline at each transect; and (e) estimation of coastline change rate at each transect. Each coastline computed was assigned five attributes, i.e., date, length, ID, shape, and uncertainty. We selected the 1989 coastline as a benchmark, as 1989 is the start of the data series and the year that the United Nations Environment Programme (UNEP) included Pakistan in the group of countries vulnerable to the impact of rising sea level [29]. Among the three existing methods for baseline demarcation (i.e., creating a baseline from a specific distance of a coastline, using a pre-existing baseline, and buffering method), we adopted the buffering method, as it is considered to be the most reliable and accurate [43].
Following Nandi et al. [43], a buffer of 200 m was generated around the baseline coastline. Two hundred meters was chosen as the optimal baseline distance after a trial-and-error process. Baseline was segmented into onshore and offshore following the coastline as best as possible, and transects were cast using the delineated coastlines and the baseline. A total of 1232 transects were cast along the coastline at intervals of 500 m. The spacing was determined after trying different options, including 200, 300, 400, 500, 600, and 700 m. The spacing of 500 m was found to be optimal, as it caused less intersection between the transects. Transect length was chosen as the maximum distance between the baseline and the farthest coastline (2000 m in the study). The transects were numbered from 1 to 1232, with 1 on the far west and 1232 at the far east of the coastline. Three different statistical models (i.e., end point rate (EPR), linear regression rate (LRR), and least median of squares (LMS)) in the DSAS program were used to compute the rate of coastline change. The EPR method calculates the coastline change between two coastlines separated by a certain time period (years), using: where d1 and d2 are the distances separating the coastlines and the baseline, and t1 and t2 are the dates when the coastlines are observed. EPR was calculated based on three time periods, i.e., 1989-1999, 1999-2009, 2009-2018, and for 1989-2018. The confidence of EPR (ECI) is calculated as: where Unc (A) and Unc (B) are the uncertainty values for coastlines A and B, respectively, Date A and Date B are the dates of coastlines A and B, respectively. The LRR method computes the coastline change by linear least-squares regression of multiple coastline positions (i.e., 1989,1999,2009, and 2018) for each transect: where d represents the distance from the baseline, X is the time interval, m is the slope of the fitted line, and b is the y-intercept. The slope of the fitted line (m) represents the rate of coastline change. The R-squared (R 2 ) value of the linear regression is calculated for each transect: where d is the known distance of an observed coastline from the baseline, dp is the predicted value of d based on the regression equation, d is the mean position of the coastline, and n is the number of observed coastlines. In this study, R 2 > 0.87 is used as the limit of certainty, following [18]. The uncertainty of the calculated rate, also known as the standard error of the slope of the regression line, is reported as confidence interval of the line of regression (LCI). The LCI is considered at confidence interval of 95% (p = 0.05), i.e., LCI95.
In the LRR method, the mean offset of the sample is used to determine the best-fit equation for the regression line. Whereas, the LMS method uses the median value of the squared residuals instead of the mean offset in order to determine the equation for the slope line.
The observed coastlines were superimposed on the corresponding Landsat images to check their accuracy. The coastline for the year 2018 was also predicted with the LRR model based on coastlines of the previous years (1989, 1999, and 2009). Correlation between the predicted and the actual coastlines demonstrates the accuracy of the observed coastline and the regression model.

Classification of Erosion
We adopted a classification scheme from Gornitz et al. [44] and divided the coastline change rate into five classes (i.e., >2.0, 1.1-2.0, −1.0 to 1.0, −1.1 to −2.0, and <−2.0 m/yr), where the positive values indicate accretion and the negative values indicate erosion. An intensity curve was plotted based on the LRR for each subzone. A linear regression statistics-based decision matrix was also devised to assess the severity of erosion/accretion in each zone and subzone [17,45].

Influencing Factors for Coastline Change
The terrain slope along each transect is measured as a ratio of elevation difference between upslope point A and downslope point B to the distance. The slope, calculated based on SRTM (Shuttle Radar Topography Mission) 30-m DEM in terms of percentage rise, was plotted and analyzed using the spatial analyst tool in ArcGIS. A topographic position tool (TPI) developed by Jenness [46] was further used to analyze the spatial variation in the nearshore topography. Nearshore topography and the location of the coastlines in 1989, 1999, 2009, and 2018 were compared to study the patterns of erosion/accretion in relation to the topography.
Statistical analysis of the mean sea level data for the period of 1916-2015, obtained from the Permanent Service for Mean Sea Level (PSMSL) archive, was carried out to evaluate the trend of SLR. A non-parametric Mann-Kendall test for statistically significant increasing/decreasing trend [47,48] and Sen's slope method to estimate the magnitude of the trend of sea level rise were applied [49,50]. The information was coupled with coastal erosion to discuss the possible impacts of SLR on coastal erosion. It is noted that the data from PSMSL considers both the climatic SLR and geological SLR. Hence, the results should be interpreted as such if any correlation between erosion and PSMSL data exists.

Coastline Delineation
The image-based coastlines of the Sindh over 1989-2018 are shown in Figure 4. The results show that almost the entire coast experienced either accretion or erosion, except for the Karachi-West, which has been converted into the artificial hard coastline. The IR West and IR East coastal regions experienced erosion up to 14.17 and 19.76 m/yr, respectively. This erosion rate is greater than that of the Karachi coastline (0.85 m/yr). The assessment of image-based delineated coastline through geometric evaluation remains the focus of scientific discussion. Image-based coastline uncertainty is subjective and probably more difficult, especially when no reference datasets are available for such a large area in undeveloped deltaic regions. In order to evaluate the coastline accuracy, initially we searched for high-resolution data such as aerial orthophotos or freely available high-resolution satellite images such as Sentinel-2 that coincide with dates of Landsat images used. Little information is found that matches the close date and no information for the same date. Therefore, results of this study have been verified by authors through personal communications with the local authorities and residents wherever possible. large area in undeveloped deltaic regions. In order to evaluate the coastline accuracy, initially we searched for high-resolution data such as aerial orthophotos or freely available high-resolution satellite images such as Sentinel-2 that coincide with dates of Landsat images used. Little information is found that matches the close date and no information for the same date. Therefore, results of this study have been verified by authors through personal communications with the local authorities and residents wherever possible.

Analysis of the Rate of Coastline Change
The obtained R 2 value between EPR and LRR is very high (R 2 = 0.94, Figure S1). Since the affinity between EPR-LRR is high, compared to that of EPR-LMS (R 2 = 0.80) and LMS-LRR (R 2 = 0.91), and because LRR is the most widely used method for coastline change rate estimation, we use the results from the LRR model to further analyze the coastline changes. The LRR method is considered the most reliable estimator of the rate of change in coastline position [51] and forecaster of trends for the future, as it minimizes the potential random error and short-term variability [52,53].
A validation step was carried out to assess the quality of the LRR-based statistical results through correlation analysis between the remotely sensed coastline and the model-predicted coastline for the year 2018. The correlation between the model-predicted and the actual coastlines was found to be 90%. Overlapped areas along the Sindh coastline between predicted and remotely sensed coastlines for the year 2018 can be found in Figure 5. Both the statistical and graphical

Analysis of the Rate of Coastline Change
The obtained R 2 value between EPR and LRR is very high (R 2 = 0.94, Figure S1). Since the affinity between EPR-LRR is high, compared to that of EPR-LMS (R 2 = 0.80) and LMS-LRR (R 2 = 0.91), and because LRR is the most widely used method for coastline change rate estimation, we use the results from the LRR model to further analyze the coastline changes. The LRR method is considered the most reliable estimator of the rate of change in coastline position [51] and forecaster of trends for the future, as it minimizes the potential random error and short-term variability [52,53].
A validation step was carried out to assess the quality of the LRR-based statistical results through correlation analysis between the remotely sensed coastline and the model-predicted coastline for the year 2018. The correlation between the model-predicted and the actual coastlines was found to be 90%. Overlapped areas along the Sindh coastline between predicted and remotely sensed coastlines for the year 2018 can be found in Figure 5. Both the statistical and graphical assessments confirmed that the LRR model-predicted coastline values were in close agreement with the remotely sensed coastline values for all the zones. Therefore, results are acceptable for further analysis.
Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 24 assessments confirmed that the LRR model-predicted coastline values were in close agreement with the remotely sensed coastline values for all the zones. Therefore, results are acceptable for further analysis.

Coastline Change Rate Statistics during the Last Three Decades (1989-2018)
All the zones exhibit both accretion and erosion along the coastline during the study period. The demarcations of these two classes (erosion/accretion) from the LRR model-based statistical output for

Coastline Change Rate Statistics during the Last Three Decades (1989-2018)
All the zones exhibit both accretion and erosion along the coastline during the study period. The demarcations of these two classes (erosion/accretion) from the LRR model-based statistical output for Karachi, IR West, and IR East zones are presented in Figure 6. Detailed information on the minimum, maximum, and average in each class (erosion/accretion) for each zone and subzone in this study is given in Table 1.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 24 Karachi, IR West, and IR East zones are presented in Figure 6. Detailed information on the minimum, maximum, and average in each class (erosion/accretion) for each zone and subzone in this study is given in Table 1. The Karachi zone falls in the low erosion category (−0.85 ± 0.45 m/yr), which is mainly contributed to by the erosion of the outlying islands (Table 3). Results show that the western parts of the Karachi coastline have been relatively less unstable as compared to other subzones during the last three decades. The IR West zone and the IR East zone (covering IDR) were highly eroded during 1989-2018 with mean erosion levels of −14.17 ± 0.45 and −19.96 ± 0.55 m/yr, respectively. The highest average erosion rate among all three zones was experienced by the subzone Kharo-Chhann, located in IR East (−27.46 m/yr), during the past three decades. The highest accretion rate (12.14 m/yr) is The Karachi zone falls in the low erosion category (−0.85 ± 0.45 m/yr), which is mainly contributed to by the erosion of the outlying islands (Table 3). Results show that the western parts of the Karachi coastline have been relatively less unstable as compared to other subzones during the last three decades. The IR West zone and the IR East zone (covering IDR) were highly eroded during 1989-2018 with mean erosion levels of −14.17 ± 0.45 and −19.96 ± 0.55 m/yr, respectively. The highest average erosion rate among all three zones was experienced by the subzone Kharo-Chhann, located in IR East (−27.46 m/yr), during the past three decades. The highest accretion rate (12.14 m/yr) is observed at the transects along the subzone Karachi South coastline, which is a potential result of coastal reclamation. While erosion is prevalent in the IDR, the change rate assessment results show that rates of the landward erosion experienced an increasing trend from the west to the east along the Pakistan coastline (Figure 7). The distance (in meters) values from the baseline at each transect are analyzed for the interdecadal variations during 1989-1999, 1999-2009, and 2009-2018. Large segments of the coastline are found to be receding landward from 1989 to 2018. This statement is supported by a linear increase in average change in distance from baseline (R 2 = 0.78, Figure 8). observed at the transects along the subzone Karachi South coastline, which is a potential result of coastal reclamation. While erosion is prevalent in the IDR, the change rate assessment results show that rates of the landward erosion experienced an increasing trend from the west to the east along the Pakistan coastline (Figure 7). The distance (in meters) values from the baseline at each transect are analyzed for the interdecadal variations during 1989-1999, 1999-2009, and 2009-2018. Large segments of the coastline are found to be receding landward from 1989 to 2018. This statement is supported by a linear increase in average change in distance from baseline (R 2 = 0.78, Figure 8).

Karachi Coastline
Karachi zone is the most developed segment along the Sindh coastline. The prominent features of the Karachi coastline are raised beaches, shallow lagoons, marine terraces, sea cliffs, and sand dunes. A closer inspection of the coastline change statistics for Karachi zone shows that the coastline oscillates between erosion and accretion ( Figure 9). The intensity plot shows that the NW-SE oriented Karachi coastline is largely disparate to the subzone Karachi South, as compared to other two subzones (Karachi West and Karachi East). Very high accretion rate (up to 8.34 m/yr) is observed at the south of Korangi Creek entrance and the southern part of Karachi coastline, which forms the coastal segment of subzone Karachi South (Figure 9). These obvious changes in land accretion at the subzone Karachi South are evident in Figure 10 and appeared prominently in the corresponding Landsat TM and OLI images in different years. Accretion measured at the Karachi South subzone is artificial because a large part of the area has been reclaimed from Arabian Sea for coastal infrastructure development [54]. Luijendijk et al. [55] also found this area near Karachi south accreting at very high rate during the 1984-2016 period. A comparison of the satellite images before and during the study period ( Figure 10) compliments the results for coastline assessment in this sector (Karachi South). Therefore, we emphasize that the coastal development and construction projects should be designed with great care, as studies confirm that the land in this subzone is susceptible to liquefaction in case of any major earthquake events.
Subzone Karachi East (adjacent to IR West zone) is experiencing high erosion (up to −2.43 m/yr). This area hosts the country's deep-water seaports, Port Bin Qasim (located towards Korangi Creek) and Port Karachi (towards Khizri Creek). Port Karachi is located more inward, whereas Port Bin Qasim is situated on the coastline along the Arabian Sea in subzone Karachi South Widening and deepening of the navigational channel to facilitate the port jetties for the passage of the cargo ships through the Phitti Creek to the seaport Muhammad Bin Qasim has resulted in the coastal erosion in this area. In 1970, when this port started functioning, the presence of barrier islands (Bundal and Buddo Islands) next to the port acted as first line of defense against coastal hazards such as storms and coastal flooding. The high erosion rate in Karachi East subzone is an indication that these natural

Karachi Coastline
Karachi zone is the most developed segment along the Sindh coastline. The prominent features of the Karachi coastline are raised beaches, shallow lagoons, marine terraces, sea cliffs, and sand dunes. A closer inspection of the coastline change statistics for Karachi zone shows that the coastline oscillates between erosion and accretion ( Figure 9). The intensity plot shows that the NW-SE oriented Karachi coastline is largely disparate to the subzone Karachi South, as compared to other two subzones (Karachi West and Karachi East). Very high accretion rate (up to 8.34 m/yr) is observed at the south of Korangi Creek entrance and the southern part of Karachi coastline, which forms the coastal segment of subzone Karachi South (Figure 9). These obvious changes in land accretion at the subzone Karachi South are evident in Figure 10 and appeared prominently in the corresponding Landsat TM and OLI images in different years. Accretion measured at the Karachi South subzone is artificial because a large part of the area has been reclaimed from Arabian Sea for coastal infrastructure development [54]. Luijendijk et al. [55] also found this area near Karachi south accreting at very high rate during the 1984-2016 period. A comparison of the satellite images before and during the study period ( Figure 10) compliments the results for coastline assessment in this sector (Karachi South). Therefore, we emphasize that the coastal development and construction projects should be designed with great care, as studies confirm that the land in this subzone is susceptible to liquefaction in case of any major earthquake events.
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 24 barriers against the erosional effects of Arabian Sea water are compromised and adversely affected. However, further evaluation in this regard is needed. Figure 9. Coastline assessment for Karachi coastline based on linear regression rate (LRR) statistical method over the period 1989-2018. The spatial location of the zone is presented in the first column. The second column represents plot of rate of accretion/erosion against each subzone labeled as "sector". The left axis on the plot shows the transect number increasing from west to east, and the top axis represents the rate of the accretion/erosion in mm/yr. Further, each subzone is classified into different levels based on the severity of erosion in that sector. Given the severity level, it is analyzed whether some action is required to tackle the erosion problem or not. The second column represents plot of rate of accretion/erosion against each subzone labeled as "sector". The left axis on the plot shows the transect number increasing from west to east, and the top axis represents the rate of the accretion/erosion in mm/yr. Further, each subzone is classified into different levels based on the severity of erosion in that sector. Given the severity level, it is analyzed whether some action is required to tackle the erosion problem or not.
Subzone Karachi East (adjacent to IR West zone) is experiencing high erosion (up to −2.43 m/yr). This area hosts the country's deep-water seaports, Port Bin Qasim (located towards Korangi Creek) and Port Karachi (towards Khizri Creek). Port Karachi is located more inward, whereas Port Bin Qasim is situated on the coastline along the Arabian Sea in subzone Karachi South Widening and deepening of the navigational channel to facilitate the port jetties for the passage of the cargo ships through the Phitti Creek to the seaport Muhammad Bin Qasim has resulted in the coastal erosion in this area. In 1970, when this port started functioning, the presence of barrier islands (Bundal and Buddo Islands) next to the port acted as first line of defense against coastal hazards such as storms and coastal flooding. The high erosion rate in Karachi East subzone is an indication that these natural barriers against the erosional effects of Arabian Sea water are compromised and adversely affected. However, further evaluation in this regard is needed.

IR West Zone
Based on the intensity of the advance/retreat curve plot (Figure 11), it is observed that the IR West zone exhibits a very high rate of erosion towards inland with average LRR-based erosion (14.174 ± 0.55 m/yr) increasing from west to east along the IDR coastal stretch. This trend is observed to be correlated with the flat topography. The IR West zone is an embayed coast of tidal mudflats, lagoons, barriers islands, mangrove canopy, and various major creeks. The sector between Waddi/Khuddi and Dabbo creeks in the Mirpur-Sakro subzone is found to be highly erosional and receding at a rate reaching up to 70 m/yr during the last three decades. This result also agrees with earlier observations from Giosan et al. [24] and Hashmi and Ahmad [57], who showed that this sector is eroding at very high and similar rates during the post-damming period, mainly due to a reduction in the sediment influx to the Indus River and SWI. Towards the farther east, the coastal sector between Dabbo and Chhan creeks covering the Ghorabari administrative unit also seems to be receding at an average rate of 15 m/yr, significantly less compared to the subzone Mirpur Sakro However, the maximum erosion rate in Ghorabari coastline segment is found to be increased (up to 50 m/yr in most parts).
The coastline segment from Chhan to Hejamaro creeks in subzone Keti-Bandar shows a similar pattern. The results show an interesting and mixed behavior of accretion and erosion at the mouths of Chhan, Hejambaro, and Turshian creeks. The western (left) sides of these creeks are found to be accreting, whereas the eastern and northeastern sides are receding. The coastline sector in front of the two active river mouths between Turshian and Khobar creeks (Keti-Bandar subzone) is found to be encroached by the seawater at a rate reaching up to 50 m/yr, consistent with the findings of [24], mainly due to the reduced sediment discharge and freshwater inflow. The SW monsoon winds further intensify the land loss processes particularly in the deltaic creeks in this zone. High tidal volume inundates the creek banks and erodes the adjoining coastal lands, which consequently increases the cross-sectional area of the seaward entrance of creeks. The erosion of creek banks destroys woody vegetation (mangroves), particularly near the open coastlines where tidal action is strong (Figure 12). The meandering of small creeks and channels due to increased SWI cut-off from regular floodwater and additional sediment deposits are also some of the causes of natural degradation of mangroves [58,59]. Mangroves are known to protect the coastline from erosion by

IR West Zone
Based on the intensity of the advance/retreat curve plot (Figure 11), it is observed that the IR West zone exhibits a very high rate of erosion towards inland with average LRR-based erosion (14.174 ± 0.55 m/yr) increasing from west to east along the IDR coastal stretch. This trend is observed to be correlated with the flat topography. The IR West zone is an embayed coast of tidal mudflats, lagoons, barriers islands, mangrove canopy, and various major creeks. The sector between Waddi/Khuddi and Dabbo creeks in the Mirpur-Sakro subzone is found to be highly erosional and receding at a rate reaching up to 70 m/yr during the last three decades. This result also agrees with earlier observations from Giosan et al. [24] and Hashmi and Ahmad [57], who showed that this sector is eroding at very high and similar rates during the post-damming period, mainly due to a reduction in the sediment influx to the Indus River and SWI. Towards the farther east, the coastal sector between Dabbo and Chhan creeks covering the Ghorabari administrative unit also seems to be receding at an average rate of 15 m/yr, significantly less compared to the subzone Mirpur Sakro However, the maximum erosion rate in Ghorabari coastline segment is found to be increased (up to 50 m/yr in most parts).
The coastline segment from Chhan to Hejamaro creeks in subzone Keti-Bandar shows a similar pattern. The results show an interesting and mixed behavior of accretion and erosion at the mouths of Chhan, Hejambaro, and Turshian creeks. The western (left) sides of these creeks are found to be accreting, whereas the eastern and northeastern sides are receding. The coastline sector in front of the two active river mouths between Turshian and Khobar creeks (Keti-Bandar subzone) is found to be encroached by the seawater at a rate reaching up to 50 m/yr, consistent with the findings of [24], mainly due to the reduced sediment discharge and freshwater inflow. The SW monsoon winds further intensify the land loss processes particularly in the deltaic creeks in this zone. High tidal volume inundates the creek banks and erodes the adjoining coastal lands, which consequently increases the cross-sectional area of the seaward entrance of creeks. The erosion of creek banks destroys woody vegetation (mangroves), particularly near the open coastlines where tidal action is strong (Figure 12). The meandering of small creeks and channels due to increased SWI cut-off from regular floodwater and additional sediment deposits are also some of the causes of natural degradation of mangroves [58,59]. Mangroves are known to protect the coastline from erosion by stabilizing the sediments through their prop root systems and dissipating the storm's energy, thus reducing vulnerability to the climate-induced hazard and increasing coastal resilience.
Remote Sens. 2020, 12, x FOR PEER REVIEW 16 of 24 stabilizing the sediments through their prop root systems and dissipating the storm's energy, thus reducing vulnerability to the climate-induced hazard and increasing coastal resilience.   reducing vulnerability to the climate-induced hazard and increasing coastal resilience.

IR East Zone
Erosive actions are found to be more pronounced in the eastern (right) side of the IDR (i.e., 19.96 ± 0.65 m/yr) as compared to the IR West zone, and the coastline is receding landward (curve plot, Figure 13). This result is in line with the findings of Waqas et al. [27]. Overall, mean erosion is found to be decreasing away from the Indus River outfall location at Khobar Creek towards Sir Creek in the direction of the littoral drift (long-shore sediment transport) from West to East, as shown in Figure 13. Deltaic soils in this zone are mostly marshy and mudflats and can easily be eroded by waves and other oceanic forces of the Arabian Sea. As a result, transects in subzone Kharo-Chhann near the IR outfall have been eroded during the last three decades at an average rate reaching up to 23.42 m/yr-the highest among all the subzones. Coastline recession assessment in this zone is consistent with 23 m/yr erosion rate estimated by the WWF-Pakistan [60] between 1952 and 2006. The rate of sediment influx and the ground displacements significantly affects the position of the coastline, particularly in the subsiding basins [61]. In the case of IDR and the areas in its vicinity, all the aforementioned factors have varied together greatly with respect to the relative SLR causing larger coastline regression and transgression [24,26,27,30,62].
A gentler topographic profile is observed along the entire IDR coastline (i.e., 0%-1.1% slope for Zone II, and 0%-0.3% slope for Zone III). For complete longitudinal profile, see Figure S3. The results indicate that the IDR is a low-lying coastal area. Based on the analysis of predefined categories of the slope position index (SPI), it is observed that most of the inland areas are at the same level as the Arabian Sea, hence, potentially letting the seawater intrude the surface water resources and salinize groundwater aquifer in this low-lying deltaic region. Slope percentage is below 1.1% for almost the entire IDR, hence, making the study region vulnerable to impacts of SLR, land erosion, and coastal flooding. Without appropriate measures, rising sea level, inundation/flooding, erosion, and SWI are likely to affect coastal productivity and economic activities in the study region.
The areas at the east of Wari Creek mouth to Sir Creek in subzone Shah-Bandar have experienced erosion at the rate of 18.41 m/yr on average. The overall erosion rate of this deltaic coast in this subzone from 1989 to 2018 is less than the average rate reported by Giosan et al. [24] for the same segment during 1978-2000. This might be due to the accumulation of the sediments on this segment due to the along-shore wave dispersal system, which transports the sediments from Indus River outfall and mouth towards the east up to Sir Creek and Rann of Kutch [63]. For this reason, average erosion rate in Shah Bandar taluka is found to be smaller than the Kharo-Chhann. Using the transect-based coastline change rate analysis method, a creek-wise status of the coastal erosion hazard was carried out. Details on minimum, maximum, and average rate of coastline change (m/yr) are given in Table 4. Transects on the left and right banks of the Waddi, Khal, Pitiani, Dabbo, Chhann, Khobar, Ghara-Chhann, Kajhar, and Sir creeks along the IDR are found to be eroding. Whereas the transects along the Turshian (one of the IR mouth), Khar, Nar, and Wari creek-banks showed a mixed sign of accretion and erosion. The highest rates of erosion are found at the transects along the banks of creeks in the IR East zone of IDR. IR West zone is protected with the dense mangrove canopy (visually verified from the comparison of the corresponding Landsat images used in this study), whereas the presence of medium and sparse mangrove canopy on the eastern side of the IR East zone has exposed it to the intense mechanical action of the Arabian Sea waves alongshore, particularly during the SW monsoon period. Moreover, there are no rivers or streams transporting the water and sediments from upstream to the coast, thereby exposing this side of the IDR to high wave action. Reduction in the mangrove forest cover and weakening of the eco-shield provided by the vegetation along with disastrous reduction in the freshwater inflow and sediment supply to the IDR [25] has led to the widening of the creeks and increased SWI [25,26].

Analysis of Influencing Factors of Coastline Change
Coastal sustainability is sensitive to coastal erosion due to its impacts on human-natural systems. Impacts of coastal erosion are likely to increase in the IDR, mainly due to the modification of the fluvial regime of the Indus River caused by terrestrial processes, such as reduction in the flow of freshwater and sediment supply to the IDR, interception of the littoral sediment transport due to

Analysis of Influencing Factors of Coastline Change
Coastal sustainability is sensitive to coastal erosion due to its impacts on human-natural systems. Impacts of coastal erosion are likely to increase in the IDR, mainly due to the modification of the fluvial regime of the Indus River caused by terrestrial processes, such as reduction in the flow of freshwater and sediment supply to the IDR, interception of the littoral sediment transport due to human interventions, intensive wave activities, and tidal inundation. As a result, the IR East zone of the IDR seems to have already lost lots of land to open sea. The rate of erosion seems to have significantly increased at most of the transects along the coast, which is also supported by other localized observations from previous studies conducted by the NIO and the WWF-Pakistan. This situation is further aggravated due to the looming climate-change-induced SLR and anthropogenic land-use change.
Considering the severity of the ongoing erosion in the study area, the effects of the climate-change-induced increment in the sea level would be manifold. It is estimated that the global sea level has risen from about 2.5 mm/yr in the 1990s to about 3.4 mm/yr today [64]. Statistical analysis of mean relative sea level at the Karachi tide gauge station for the period 1916-2015, obtained from the PSMSL archive, depicted an increasing trend of MSL along the Sindh coast. The Mann-Kendal and Sen's slope estimation test result indicates that the sea level has been rising quite rapidly (i.e., 1.9 mm/yr, p-value >95%, Figure 14). This is further complemented by the linear regression analysis of the PSMSL sea level time-series (confidence 75%) ( Figure 14). The increasing trend in the sea level along the Pakistan coastline is higher than that reported by Farah and Meynell [65], Figure 14. This situation can lead to more substantial land encroachment by the seawater, accelerating erosion and resulting in higher investment costs for protective measures. This increasing SLR trend coupled with ongoing coastal erosion during last the three decades could be a potential threat to coastal communities and ecosystems, and the recession might get worse under the future SLR. This situation can result in more coastal areas being subject to inundation [33] if proper measures are not taken. Reports from NIO, Pakistan state that the reclamation, given the climatic and geological sea level rise, has already changed the pattern of energy dissipation along the Clifton beach (subzone-Karachi East), causing inundation of roads in southwest monsoon periods during spring tides [66].
localized observations from previous studies conducted by the NIO and the WWF-Pakistan. This situation is further aggravated due to the looming climate-change-induced SLR and anthropogenic land-use change.
Considering the severity of the ongoing erosion in the study area, the effects of the climatechange-induced increment in the sea level would be manifold. It is estimated that the global sea level has risen from about 2.5 mm/yr in the 1990s to about 3.4 mm/yr today [64]. Statistical analysis of mean relative sea level at the Karachi tide gauge station for the period 1916-2015, obtained from the PSMSL archive, depicted an increasing trend of MSL along the Sindh coast. The Mann-Kendal and Sen's slope estimation test result indicates that the sea level has been rising quite rapidly (i.e., 1.9 mm/yr, p-value >95%, Figure 14). This is further complemented by the linear regression analysis of the PSMSL sea level time-series (confidence 75%) ( Figure 14). The increasing trend in the sea level along the Pakistan coastline is higher than that reported by Farah and Meynell [65], Figure 14. This situation can lead to more substantial land encroachment by the seawater, accelerating erosion and resulting in higher investment costs for protective measures. This increasing SLR trend coupled with ongoing coastal erosion during last the three decades could be a potential threat to coastal communities and ecosystems, and the recession might get worse under the future SLR. This situation can result in more coastal areas being subject to inundation [33] if proper measures are not taken. Reports from NIO, Pakistan state that the reclamation, given the climatic and geological sea level rise, has already changed the pattern of energy dissipation along the Clifton beach (subzone-Karachi East), causing inundation of roads in southwest monsoon periods during spring tides [66]. Although the constant evolution of the coastline in the study region is controlled by the multiplicity of terrestrial and marine factors, human interventions are also a major factor affecting the balance of the deltaic region. Concerns have also been raised on the implications of coastal development projects such as Bodha Island City TK's Bodha Islands City Project (https://thomaskramer.com/tks-bodha-island-city-project/) and Zulfiqarabad (Zulfiqarabad Project-Impact on Environment (https://www.dawn.com/news/702023)) on the coastal dynamics of the study area. Results from this study can be utilized in the planning and design of these projects, considering the erosion in the development process. Based on the evidence on coastal erosion hazard in the region presented in this study, we emphasize that the design and implementation of such new Although the constant evolution of the coastline in the study region is controlled by the multiplicity of terrestrial and marine factors, human interventions are also a major factor affecting the balance of the deltaic region. Concerns have also been raised on the implications of coastal development projects such as Bodha Island City TK's Bodha Islands City Project (https://thomaskramer.com/ tks-bodha-island-city-project/) and Zulfiqarabad (Zulfiqarabad Project-Impact on Environment (https://www.dawn.com/news/702023)) on the coastal dynamics of the study area. Results from this study can be utilized in the planning and design of these projects, considering the erosion in the development process. Based on the evidence on coastal erosion hazard in the region presented in this study, we emphasize that the design and implementation of such new coastal projects in a highly erosive coastal region should consider the multidisciplinary approach in order to ensure the project success and minimize the adverse impacts to the environment and coastal systems.
We observe the erosion to be increasing eastward (Figure 7) while the slope of the land elevation (percent) is decreasing eastward (Figure S3), threatening the eastern parts of Karachi city where one of the major deep seaports of the country, Port Qasim, is located. The flatter the topographic slope, the higher the vulnerability to erosion and SWI. The topographic slope-profile analysis shows that subzone Karachi West in the far west has a steep (upward) gradient (0%-15.5%), which manifested itself as sharp hike in percentage rise in the profile graph as compared to Karachi South and Karachi East subzones (0%-3.5%) towards the east ( Figure S4). Due to the extreme land-erosion in this region, the fertile lands have already been converted seabed, which is evident in the SPI map, as most of the inland areas are at the same level as the Arabian Sea ( Figure S4). This situation is likely to adversely affect the livelihood of the inhabitants. However, further evaluation in this regard is a prerequisite for appropriate decision making. Our study can act as a road map for such further evaluations, as it pin-points the areas of different levels of erosion along the Sindh coast ( Figure 9, Figure 11, and Figure 13).

Implications for Coastal Sustainability in Pakistan
Virtually less aggradation or accelerated progradation/recession along with the increasing trend of the SLR has already made the entire IDR as well as the eastern end of the Karachi coastline vulnerable to coastal hazards. As evident from the results, the region is highly susceptible to erosion hazard, which might lead to intense SWI and further inland-flooding. In order to keep the seawater off the coasts along Karachi coastline, the natural coastal barriers such as barrier islands need to be protected/strengthened and other artificial barriers need to be constructed where required. To further reduce risk posed to the humans and environment in the IDR (IR West zone and IR East zone), we suggest nonengineered measures through building-with-nature such as plantation of salt-tolerant mangroves in order to curb the coastal recession and SWI, leading to coastal sustainability in the study area [67]. Restoration of the lost mangroves and conservation of the remaining are integral for stabilizing the coastline [60]. Special seeding techniques might be required for revegetation of the mangroves and their survival. While they reduce coastal erosion due to their soil binding characteristics through their U-shaped prop root system, they also protect the coastal communities from coastal natural hazards (i.e., tropical cyclones and storms), increasing the overall coastal resilience [68]. Current knowledge about the coastline and coastal infrastructure protection seems adequate to derive general guidelines on the role of the coastal vegetation for shielding coastal assets from coastline recession/progradation. However, the effectiveness of the protective role of the vegetation species in combating erosion is site-specific [69]. Adoption of hybrid (mixed engineered-natural) structural measures can also provide increased benefits to coastal systems in the study region. In this way, the present study facilitates promoting risk-informed planning, place-based management, and coastal resilience for long-term sustainability. Taking ecosystem-based measures to reduce the coastal erosion will also assist in achieving sustainable development goals (SDGs) 6, 14, 13, and 15, due to the provisioning of several services from coastal natural systems [70]. Since changes in the coastline morphology alter the flooding hazard, and since long-term coastal management depends on the changing coastline position, simultaneous occurrence of coastal erosion and flood hazard can be jointly studied for comprehensive coastal hazard assessment in Pakistan.
The current study provides a basis for such assessments (i.e., the regions where the coastal erosion is higher).

Limitations and the Way Forward
We use the SRTM DEM to analyze the slope at the coast. However, other information derived from DEM, such as coastline curvature and ground subsidence in the coastal regions, can also be used to estimate the contribution of sea level rise to the coastline changes, and thus separate the global climate drivers from the local anthropogenic (coastal infrastructure and use) and other climatic factors (e.g., torrential rainfalls). However, this may require additional resources, computational capabilities, and employment of different, mixed approaches. A preliminary case study on the ground displacement was carried out by Kanwal et al. [56] for Karachi city. This can be done on a large scale and combined with the erosion phenomena for a comprehensive coastline change quantification as well as vulnerability mapping of coastal communities to compound hazards. This could be a potential research question for future studies. Moreover, the shoreline changes along the barrier islands (BIs) could also be evaluated using our framework to further support the coastal sustainability and resilience enhancement, as the BIs are known as the important defense against coastal storms and provide multiple services to support socio-ecological systems.

Conclusions
This study evaluates the spatial heterogeneities in the large-scale morphological changes along the Pakistan coastline on zonal and subzonal (administrative unit smaller than district level) scales during the past three decades . The rates of the coastal erosion are quantified at transect level (total 1232, with a 500-m spacing between transects) and reported as average statistics for each zone and subzone along the Sindh coastline of Pakistan. This multilevel spatial information of coastal erosion from this study is important for coastal planners and decision-makers for effective coastal planning, taking appropriate actions, and adaptation measures to mitigate the erosion process for coastal sustainability, particularly in the regions which are highly vulnerable to coastal erosion. This first-of-its-kind assessment is essential as well as useful for integrated management and planning of the coastal areas in Pakistan. We find that the Sindh coastline, particularly the IDR coastline segment, is highly erosional, which is critical in the context of coastal infrastructure and communities.
Overall, the present trend of coastline retreat will continue to modify the eastern end of Karachi (including Korangi Creek and outlying BIs) due to both human intervention and marine processes if appropriate measures are not taken. Relevant agencies and the policymakers should reinforce the prevention measures on a priority basis (particularly in high erosion areas) so that the inland SWI could be curtailed in the IDR, which is home to thousands of migratory birds, a source of livelihood for millions, and a cultural heritage site. While the results from this study are important for effective coastal planning and implementation of appropriate measures to reduce coastal erosion and the associated risks in Pakistan, the study acts as a road map to scale-down the regions with higher coastal erosion for further in-depth assessments.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/1/8/s1, Figure S1: Spectral responses of coastline configuration. Dotted lines show the location of the instantaneous coastline whereas the solid lines mark the high tide coastline, delineated along the sea-land boundary of Sindh Coast (left to right). Figure S2. Comparison of coastline change rates (m/year) obtained by different statistical methods is presented on right, (a) EPR vs LRR; (b) EPR vs LMS; (c) LRR vs LMS for the overall Sindh coast. Note: EPR, end point rate; LMS, least median of squares; LRR, linear regression rate. Figure S3. Plot of percentage rise measured for the longitudinal profile of Sindh Coast (unit: km). Slope (percentage rise) varies from 15.5% in the west near Ras Mauri to below 0.5% near Sir Creek. Slope percentage is below 1.1% for almost the entire IDR region. Figure S4. Slope percentage map of the Sindh coastal region. Black line marks the areas of coastline for which profile graph analysis is done and give below. Table S1: Global EPR and LMS statistics of all the study area coastline change rate (m/yr) during (1989-2018) period.
Author Contributions: This research was supervised by X.D. S.K. designed the research experiments, analyzed the results, and wrote the initial draft. S.K., X.D. and M.S. contributed towards the discussion of results and writing, review, and editing of the manuscript; S.A. provided review of the draft of the manuscript. All authors have read and agreed to the published version of the manuscript.