Coastline Vulnerability Assessment through Landsat and Cubesats in a Coastal Mega City

Abstract

According to the Intergovernmental Panel on Climate Change (IPCC), global mean sea levels may rise from 0.43 m to 0.84 m by the end of the 21st century. This poses a significant threat to coastal cities around the world. The shoreline of Karachi (a coastal mega city located in Southern Pakistan) is vulnerable mainly due to anthropogenic activities near the coast. Therefore, the present study investigates rates and susceptibility to shoreline change using a 76-year multi-temporal dataset (1942 to 2018) through the Digital Shoreline Analysis System (DSAS). Historical shoreline positions were extracted from the topographic sheets (1:250,000) of 1942 and 1966, the medium spatial resolution (30 m) multi-sensor Landsat images of 1976, 1990, 2002, 2011, and a high spatial resolution (3 m) Planet Scope image from 2018, along the 100 km coast of Karachi. The shoreline was divided into two zones, namely eastern (25 km) and western (29 km) zones, to track changes in development, movement, and dynamics of the shoreline position. The analysis revealed that 95% of transects drawn for the eastern zone underwent accretion (i.e., land reclamation) with a mean rate of 14 m/year indicating that the eastern zone faced rapid shoreline progression, with the highest rates due to the development of coastal areas for urban settlement. Similarly, 74% of transects drawn for the western zone experienced erosion (i.e., land loss) with a mean rate of −1.15 m/year indicating the weathering and erosion of rocky and sandy beaches by marine erosion. Among the 25 km length of the eastern zone, 94% (23.5 km) of the shoreline was found to be highly vulnerable, while the western zone showed much more stable conditions due to anthropogenic inactivity. Seasonal hydrodynamic analysis revealed approximately a 3% increase in the average wave height during the summer monsoon season and a 1% increase for the winter monsoon season during the post-land reclamation era. Coastal protection and management along the Sindh coastal zone should be adopted to defend against natural wave erosion and the government must take measures to stop illegal sea encroachments.

1. Introduction

A shoreline is defined as the boundary between land and water. The position of the shoreline is dynamic both spatially and temporally, due to hydrological, geological, climatic, and economic developments in coastal areas [1,2]. Therefore, shoreline indicators are used for a consistent representation of a shoreline. The most common shoreline indicators are the tidal datum (e.g., a specific elevation at the land–ocean boundary) and visually discernible features (e.g., a revetment structure, an erosion scarp, a previous high-tide high-water level, low-tide low-water level, or dune vegetation line, etc.) [3,4,5,6,7,8]. The shoreline is an indicator of the ecological health of coastal areas [3,9,10,11]. In past decades, the development of mega-infrastructure projects in coastal areas for urbanization/industrialization has placed more stress on the coastal ecosystem and changed ocean hydrodynamics [12,13,14,15,16]. Unplanned coastal development coupled with the natural action of ocean processes has led to the increasing vulnerability of low lying coastal environments [16,17,18]. This also induces spatiotemporal changes in the shoreline. Short-term geomorphological changes in the shoreline are caused by extreme geological, climatic, and oceanic events (i.e., earthquakes, tsunami, seasonal variation in waves, tides, and storms conditions), hence, such changes are less predictable, while long-term changes to the shoreline are caused by relative changes in astronomical, meteorological, and regional climatic variations (i.e., tides, waves, sea-level rise, and storm surges), and are somewhat predictable [6,19]. Both types of shoreline change are important for understanding trends in coastal sustainability for different times and spaces [3].

Anthropogenic activities coupled with global warming intensify ocean controlling factors that modify the coastal environment [20,21]. Morphological shoreline changes are directly associated with wave height, tide level, wind speed, periodic storms, and sea level [15,22]. Specifically, an accelerated rise in sea level and cyclones pose a considerable threat to populations living within a 100 km vicinity of the shoreline, as well as their economic assets. It is considered that more than 40% of the world’s population is living in coastal cities [18]. For the monitoring of coastal zones, a common method for the assessment of spatiotemporal variability of the shoreline consists of the extraction of the shoreline position using multi-temporal data sources [5,6,22,23]. The primary data sources used for shoreline extraction and evaluation are historical photographs, coastal maps/charts, aerial and beach surveys, global positioning system based shoreline surveys, and images acquired through satellite sensors [15,24,25]. Each source has its own measurement and positional uncertainty [1,4,7,26]. Historical archives of satellite remote sensing open source data with a high spatial resolution of 3 m such as from Planet Scope to a low spatial resolution of 1000 m from the Moderate Resolution Imaging Spectroradiometer (MODIS) have revolutionized science. These datasets have proved their potential for the detection of historically changing trends in extensive coastal land masses, with an insignificant uncertainty in the case of high spatial resolution, and a relatively high uncertainty in cases of moderate to low spatial resolution [4,12,15,16,27,28,29,30].

Different methods, including image classification [4,31], band ratio [15],principal component analysis [32], overlay operation [30], and density slicing [23] have been used for the delineation, mapping, and estimation of the variability of proxy-based shoreline positions through satellite imagery, and have differing levels of measurement accuracy [4,5]. The most commonly used indices for the delineation of the shoreline include the Normalized Difference Water Index (NDWI) [33], the Automated Water Extraction Index (AWEI) [34], and the Modified Normalized Difference Water Index (MNDWI) [35]. Several studies have been carried out globally, of both short-term and long-term shoreline changes using different image analysis techniques. For example, a recent study by Nassar et al. (2018) [36] used the archive of Landsat multi-sensor images from 1989 to 2016 to assess the rates of change by the accretion and erosion of the coastal area of North Sinai that were caused by changes in hydrodynamics and alongshore currents. Cenci et al. (2018) [37] performed a multi-proxy analysis of two different shorelines (i.e., in high and low energy areas) and assessed the rate of shoreline change to highlight the uncertainty in coastal risk management. Bheeroo et al. (2016) [38] carried out a risk assessment of the north-western side of the Mauritius coastline and addressed changes in the coastal environment with reference to high energy wave action and storms. Kermani et al. (2016) [39] used aerial photographs and high resolution satellite images of the Bay of Jijel in Eastern Algeria to investigate changes in morpho-dynamics of the shoreline. Ozturk and Sesli (2015) [15] used historical multi-temporal satellite data to assess changes in the coastal environment in the lagoons of Kizilirmak Delta due to ocean controlling factors (e.g., wind and wave erosion). They also addressed the rate of shoreline change and shrinkage of the lagoon series. Kuleli et al. (2011) [40] assessed the rates of shoreline change on the Ramsar beaches of Turkey using Landsat images from 1975 to 2009 and found that most parts of the wetland areas were under stress with a significant withdrawal rate.

Pakistan is considered the seventh most vulnerable country to climate change induced events. These include four major floods and six historical cyclones, which have damaged the socio-economy of low lying areas within the Sindh coastal region [41]. The coastal belt of the Indus Delta Region is continuously changing due to physical processes such as tidal actions, waves, wind speed, sea level rise, and anthropogenic factors, for example land reclamation and modification [12,16,42,43]. Rising sea levels are contributing to the erosion of less resistant soil in the region. It has been claimed that the mean sea level rose from 1.7 mm/year between 1900 and 2010 and 3.2 mm/year globally over the past decade [44]. Similarly, the Intergovernmental Panel on Climate Change (IPCC) predicts that the global mean sea-level will rise from 0.43 m to 0.84 m in the current scenario and 0.61 m to 1.10 m in the worst case scenario by 2100 [45]. It has been reported that the sea level has already risen from 1.1 mm/year to 1.8 mm/year regionally over the past decade [42]. If sea levels rise above the predicted rate, shorelines in the region will experience erosion in both the long- and short-term. Unfortunately, there has been no study to quantitatively estimate future shoreline changes along the coastline of Karachi, Pakistan. Therefore, this study aims to (i) merge multi-temporal datasets i.e., topographic maps and medium (30 m) to fine spatial resolution (3 m) satellite imagery, (ii) investigate a change in the position of the Karachi shoreline from 1942 to 2018, and (ii) assess the rate of systematic land loss and/or gain due to oceanic processes and anthropogenic activities along the coastline.

2. Study Area

Karachi is considered Pakistan’s largest city, hosting about 7% of the country’s population, making it the fifth largest coastal city in the world [46]. The coast of the Karachi city is about 100 km between the Gharo creek on the east and the river hub on the west, and constitutes a number of tourist beaches [43,47]. The coastline of Karachi is considered an economic hub for Pakistan as 90% of seaborne trade is carried out through two of its international ports i.e., the Port of Karachi and the Port Qasim [43] (Figure 1). The elite class of Pakistan’s population enjoy life along the Karachi coast, and a number of tourist spots have been developed, which has increased its residential property values [16,47]. The study area, the coast of Karachi, is divided into two zones i.e., eastern and western zones. The eastern zone starts from the Defense Housing Authority (DHA) phase 8 to the South Asian terminal on the west and is 25 km long (Figure 1), while the western zone from Manora beach on the east to the Engro beach huts on the west (Figure 1) is 29 km long. The western zone is mostly sandy but toward its western part the nearshore region becomes steeper, with irregular rocky outcrops associated with the underlying geology. During the last decade, several housing societies including the DHA Phase 8 and EMMAR project have encroached on the coastline in the Karachi metropolitan and Crescent Bay areas. The coastline of the eastern zone constitutes a mostly loose stone structure and sandy beaches.

The Sindh coastal zone is mixed wave dominated, with continuous accretion and erosion due to waves and tides of up to 3 m in height, modifying the coastal environment [16]. Strong westerly winds prevail during the summer monsoon and influence ocean circulation with the dominant direction of surface water flow [42]. The pre- and post-monsoon periods have long been associated with cyclone development in the Arabian Sea, which increase erosion and move the shoreline landward. The alongshore sediment transportation is from west to east due to a longshore drift in a clockwise direction. Water depth within 1 km offshore in the study area is shallow (3 m–10 m) with a maximum depth of about 10 m at Hawk’s Bay. Seawater surface temperatures in the near shore water off the coast of Karachi range from 24 °C to 28 °C in summer and 20 °C to 24 °C in winter [43].

3. Data Used

3.1. Historical Topographic Maps

In order to use a long-term data record for the monitoring and evaluation of the shoreline, two topographic maps by the Army Map Service (GDPE), Corps of Engineers, at a 1:250,000 scale for the years 1942 and 1966 (AMS U502 NG 42-13, Edition 1), were used.

3.2. Satellite Data

Two types of satellite datasets were used to assess morphological changes in the shoreline positions i.e., data from Landsat at medium spatial resolution (and Planet Scope at a finer spatial resolution). Landsat collection 1 Level-1 images of Multispectral Scanner (MSS) sensor and Level-2 images of Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Operational Land Imager (OLI) sensors were acquired on different dates (

Table A1. List of satellite images used in this study.

Sr. #	Satellite	Sensor	Scene Identifier	Path/Row	Acquisition Date
1	PlanetScope	Optical	20180422_053509_1033_3B_ AnalyticMS_SR	53/509	22 Apr 2018
2	PlanetScope	Optical	20180422_053510_1033_3B_AnalyticMS_SR	53/510	22 Apr 2018
3	PlanetScope	Optical	20180422_060440_0f32_3B_AnalyticMS_SR	60/440	22 Apr 2018
4	PlanetScope	Optical	20180422_060441_0f32_3B_AnalyticMS_SR	60/441	22 Apr 2018
4	Landsat-5	TM	LT05_L1TP_152043_20110211_20161010_01_T1	152/43b	11 Feb 2011
5	Landsat-7	ETM+	LE07_L1TP_152043_20021211_20170127_01_T1	152/43b	11 Dec 2002
6	Landsat-5	TM	LT05_L1TP_152043_19900217_20170131_01_T1	152/43b	17 Feb 1990
7	Landsat-2	MSS	LM02_L1TP_163043_19761214_20180425_01_T2	163/43a	14 Dec 1976

Note: The superscripts a and b represent the path/row for WRS-1 and WRS-2, respectively.

) from the United States Geological Survey (USGS) Earth Explorer (http://earthexplorer.usgs.gov/) and USGS Earth Resources Observation and Science (EROS), as well as the Center Science Processing Architecture (ESPA) On Demand Interface (https://espa.cr.usgs.gov/), respectively. The finer spatial resolution data from Planet Scope (PS) was acquired for the year 2018. Four continuous strips of PS (Table A1) scenes covering the study area and Level 3B product (orthorectified and surface reflectance image product) of PS were acquired.

3.3. Wind-Wave Data

Winds are the driving source of waves and the transformation of wind energy results in strong ocean currents [42,48,49]. Wave action is considered the most important factor for a change in a coastal environment during the summer and winter monsoon seasons, as prevailing winds from the south-west and northeast intensify attacks on the coastline. Long-term offshore daily-hourly wind-wave data collected from National Oceanic and Atmospheric Administration (NOAA) Wave-Watch Ⅲ Model [50] were used to evaluate the influence of wind-wave climatology on the coastline of Karachi before and after the mega land reclamation project of 2008 to 2017.

3.4. Tidal Data

Tides significantly affect wave dominant shorelines [49,51] as they induce longshore drift, which can nourish or erode the coast. The variations of tides at Karachi can be used to analyze the impact of tidal currents along the coastline and identify the coastal waterlines. In this study, tide height was used as a tool in the image selection criteria to extract the true shoreline position. Port Muhammad Bin Qasim is located about 5 km east of the coast of Sindh where the tide height has been monitored by the National Institute of Oceanography (NIO) Karachi, Pakistan from 1976 to 2018. Historical tide height records show that maximum and minimum tide heights recorded between 1976 and 2018 at the station were 3.8 m and a little less than 0.1 m.

4. Methodology

4.1. Pre-Processing of Data

4.1.1. Rectification of Topographic Maps

The geo-referencing of scanned topographic maps was made in-house using the Universal Transverse Mercator (UTM) coordinate system (WGS1984-UTM-Zone-42N; EPSG-32642 OTF) to project over the study area so as to estimate the shoreline position. The topographic sheets were then rectified using 35 Ground Control Points (GCP) and the NIO shoreline profile along with the 2016 image to reduce measurement uncertainty in the shoreline position [4,26,52].

4.1.2. Atmospheric Correction of Satellite Imagery

The pre-processing of images from different sensors in order to make it consistent with other images by normalizing them is very important in shoreline change detection studies [4,7,31,53,54]. The image based Dark Object Subtraction (DOS) model was used to convert the digital number recorded by the MSS sensor to generate a surface reflectance product [55,56,57]. Then, the surface reflectance products of L2 MSS images with a 60 m spatial resolution were resampled at 30 m to conform to the spatial resolution of other Landsat sensors (TM, ETM+, and OLI) [16,56].

The Scan Line Corrector (SLC) error encountered on all images acquired after 31 May 2003 from ETM+ caused about a 22% loss of data [57]. Therefore, a correction of the SLC error was carried out using the “Fill no data tool” under the raster tools available in Quantum Geographical Information System (QGIS 2.8.8) software, though some degree of error remains [4,26,52].

4.2. Satellite Data Selection Criteria

Images for the same months of different years for several decades enabled us to investigate the variability of shoreline morphology with nominal inherited measurement and positional uncertainty [15,43,58,59]. Therefore, we only selected Landsat (L2, L5, L7, and L8) and PS images that met the following conditions: (1) 0% cloud cover over the study area, (2) non-flooding months (December to April), and (3) tide height of <1 m within a time window of ±1 h of the image acquisition time. Non-flooding months provided information about the stability of the ocean condition because tropical storms mainly prevail during the months of May to June and September to November [16]. Furthermore, low tide height provided stable outer boundaries of the coastline [4,5,26,59]. This image selection criteria (Figure 2) resulted in a limited number of images available, including four images from Landsat sensors and one from PS which was prioritized over the OLI images due to a higher resolution (Table A1). The image selection criteria in this study was helpful in reducing or removing the inherited seasonal uncertainty in shoreline change analysis.

4.3. Shoreline Identification and Extraction

The basic requirement for detection of the true shoreline position is to determine the shoreline indicator i.e., the high tide line, the high-water line, dune line, or vegetation line [3,4,5,7,26,60]. These indicators show a variability in the shoreline position for the same time and space [3,4,5]. For the present study, detection of a high water line (HWL) under stable oceanic conditions (low wave-tide condition) from satellite images was calculated by band ratioing. The HWL shoreline indicator for images acquired during non-flooding months at a low tide condition reduced the overall positional uncertainty to a minimum as HWL migrates to contemporary low water line (LWL) or mean water line (MWL) to draw satisfactory continuity of the shoreline positions during different time periods [26,27]. The best results for shoreline identification were obtained from NDWI (Equation (1), which is used to monitor changes in water content) by evaluating the results of all indices (NDWI, MNDWI, and AWEI) [33,34,35] visually with true color satellite images. Mud-flats were also considered to represent the shoreline [48].

$$NDWI=\frac{Green-NIR}{Green+NIR}$$

(1)

The land and water mask were then developed by image thresholding [61]. Raster binary images were then converted to linear vectors by using the raster to vector conversion tool available in QGIS 2.8.8 software. Minor edits were made to the shoreline position by visual interpretation at a scale of 1:20,000 for mapping the shoreline position with more accuracy than the pixel size [62,63]. The maximum measurement uncertainty for the MSS sensor was presumed as ±37 m (about one pixel of resampled MSS image) while for TM and ETM+ sensors it was ±34 m (about one pixel of TM/ETM+) and for PS ±3.4 m (about one pixel of PS). Historical shoreline positions for 1942 and 1966 prior to the available satellite images were acquired from rectified GeoRef topographic sheets through digitization using GIS state of the art tools at a scale of 1/8000, having a maximum measurement uncertainty of about ±37 m [26,52,62]. The MWL for topographic maps could be located to within 0.150 mm on a map scale or to within 37 m on the ground for a map of scale 1:250,000. The comparison of MWL between topographic sheets and HWL for satellite images was possible due to their near equivalence [52]. The best estimate of maximum annualized average shoreline uncertainty for the study was considered ±0.43 m. Here it should be noted that the uncertainty 0.43 m/year is an average value for change in rate i.e., not constant for all transects of the local area. For example, in the case of Hawks Bay (in the western zone), only 105 transects out of the total 268 transects encounter accretion and erosion less than 0.07 m/year, which is less than the annualized uncertainty rate of 0.43. Furthermore, it needs to be clarified that the value of 0.07 m/year (by the EPR method) reflect the averaged value of shoreline change (for all transects drawn including erosional and accretional) over the period of 76 years for the Hawks Bay area (

Table 1. Average shoreline change rate of the eastern and western zone.

Eastern Zone
	DHA ZB	DHA GC	DHA ZC	DHA ZD	DHA EXT	S V	SJC	SAT
TT	86	104	65	236	149	257	162	135
T ID	1–86	87–190	19–255	256–491	492–640	641–897	898–1059	1060–1194
ShoreT	RS	RS	RS	RS	SGB	SB	SedB	RS
EPR	14.3	22	22.4	21	8.7	5.9	12	8.38
LMS	2.51	5.05	6.98	16.98	11.54	7.2	4.78	7.01
LRR	16.29	24.96	23.95	25.8	11.26	5.8	9.95	5.57
SCE	1.73	1.73	1.77	1.89	0.90	0.57	1.01	0.82
NSM	1.63	1.68	1.71	1.59	0.67	0.45	0.91	0.64
Western Zone
	MB	SS	KP	HB	JG	SG	ARG	HAG	EBH
TT	342	291	62	268	77	46	40	96	100
T ID	1–342	343–633	634–695	696–963	964–1040	1041–1086	1087–1126	1127–1222	1223–1322
ShoreT	SGB	SB	RB	SB	SB	SB	SB	SB	RB
EPR	−0.61	−0.77	0.18	0.07	−1.55	−0.95	−0.73	−0.75	−0.87
LMS	−0.07	−0.19	0.12	−0.01	−1.7	−0.54	−0.82	−0.28	−1.34
LRR	−0.32	−0.72	0.07	−0.22	−1.79	−0.65	−0.12	−0.74	−1.2
SCE	0.23	0.10	0.10	0.11	0.15	0.11	0.14	0.17	0.13
NSM	−0.05	−0.06	0.01	0.01	−0.12	−0.07	−0.06	0.17	−0.07

Notes: (1) TT = Total Transects drawn, T ID = Transect Identity, ShoreT = Shoreline Type; (2) RS = Riprap Structer, SB = Sand Beach, SGB = Sand and Gravel Beach, SedB = Sediment Bank, RB = Rocky Beach; (3) Eastern Zone: DHA = Defense Housing Authority, ZB = Zone B, GC = Golf Club, ZC = Zone C, ZD = Zone D, Ext = Extension, SV = Sea-view, SJC = Shireen Jinnah Colony, SAT = South Asian Terminal; Western Zone: MB = Manora Beach, SS = Sandspit, KP = Kaka Pir, HB = Hawks Bay, JG = Jamali Goath, SG = Somar Goath, ARG = Abdul Rehman Goath, EBH = Engro Beach Huts; (4) EPR, LMS, LRR = Rate of shoreline change (m/year) and NSM, SCE = Shoreline movement (Km).

), and do not reflect a change within a pixel size. These low rates in change represent very low or moderate anthropogenic activity compared to the eastern zone where there are more anthropogenic activities (Section 5).

4.4. Digital Shoreline Analysis System (DSAS)

The Digital Shoreline Analysis System (DSAS) version 4.4 is now provided by the Environmental Systems Research Institute (ESRI) as an add-in utility for the ArcGIS version 10.5 service pack, which was developed by the USGS under the Coastal and Marine Geology Program [64]. The statistical algorithm in DSAS allows the measurement of shoreline change statistics for multiple historical shorelines at each transect drawn at a user specified interval and length, relative to the baseline. The baseline is drawn at a user specified distance from the source shoreline, which remains parallel and offset to the reference shoreline orientation. Then transects were cast, which must intersect each shoreline position to calculate the change in shoreline positions. The DSAS executes statistical operations based on measured differences between historic shoreline positions to measure the rate of change in the following terms:

• SCE: Shoreline change envelopes (distance measurement);
• NSM: Net shoreline movement (distance measurement);
• EPR: End point rate (point rate);
• LRR: Linear regression rate (regression statistics);
• WLR: Weighted linear regression (regression statistics);
• LMS: Least median of squares (advanced statistics).

The extracted shorelines (i.e., from the historical topographic sheets and multi-temporal images) and the relative discrete baselines were placed in a personal geo-database prepared using Arc Catalog and prerequisite attribute fields (i.e., date, shape length, cast direction, etc.) were added to them. The movement of shoreline recorded at each transect drawn relative to the orientation of baseline, was considered as a landward movement if there was a negative (–) value recorded at the transect drawn, while positive (+) values indicated seaward movement (accretion) at that specific transect [64]. In this study, NSM, SCE, EPR, LMS, and LRR metrics were used to determine changes in the shoreline during the study period. Here SCE and NSM measured the net change in historical shoreline positions (distance) while EPR was used to measure the net rate of change in the shoreline [15]. The LMS and LRR were used to calculate the median and mean susceptibility rate of the shoreline change by fitting a regression line for each transect [64]. The SCE measured the distance between the two farthest shoreline positions and the distance value from the SCE method was always positive, as it does not record the accretion and erosion in shoreline positions. NSM calculated the distance of the youngest from oldest shoreline position, which represents the change in positive and negative values [40,64]. SCE was independent of the date factor of shoreline position while NSM was associated with the date factor for these two shoreline positions, and both estimated the distance, but not the rate. The EPR was used to estimate the rate of erosion and accretion per year at each transect [52]. In LMS, the median value of the squared residual was used for fitting a regression line shoreline. Its value was obtained by dividing the total coastal change distance, or NSM, by the total number of years at each transect [52]. In LMS, the median value of the squared residual was used for fitting a regression line relative to different positions of shoreline at each transect, which exhibited the rate of shoreline change. Similarly, the LRR method used the mean value of residuals to fit a line for each transect in which the slope of the line predicted the rate of shoreline change. A great variation between these two methods was observed because in the LRR method, each input offset data point had an equal influence on the estimation of the slope’s change rate, while in the LMS method each offset data point had less influence on the slope as the median value was considered the rate. The LRR rate of shoreline change was susceptible to effects of outlier values because it could underestimate or overestimate shoreline change rate values observed by the EPR method [58].

4.5. Shoreline Change Analysis

Different parts of the study area showed different characteristics of coastal change [7,65]. Therefore, we divided the study area into two major zones as eastern and western zones of a length of 25 km and 29 km, respectively. These major zones were further divided by their local area name to better understand shoreline change at a micro level. The eastern zone had local areas named DHA Zone B, DHA Golf Club, DHA Zone C, DHA Zone D, DHA Extension, Sea-view, Shireen Jinnah Colony, and South Asian Terminal (Figure 3). The shorelines of DHA Zone B, DHA Golf Club, DHA Zone C, DHA Zone D, and South Asian Terminal have loose stone structures while DHA Extension has a sand and gravel beach (Table 1). The shoreline of Sea-view has a tidal flat sand beach and Shireen Jinnah Colony has a sediment bank (Table 1). Similarly, the western zone has local names Manora Beach, Sandspit, Kaka Pir, Hawks Bay, Jamali Goth, Somar Goth, Abdul Rehman Goth, Haji Ali Goth, and Engro Beach Huts (Figure 4). The shoreline of Sandspit, Hawks Bay, Jamali Goth, Somar Goth, Abdul Rehman Goth, and Haji Ali Goth have sandy beaches while Kaka Pir and Engro Beach Huts have rocky beaches (Table 1). The shoreline of Manora Beach constitutes sand and gravel beaches (Table 1).

A baseline was created for the two zones offshore and onshore of the coast, which remained parallel to the shorelines, in order to draw transects for determining shoreline changes. The baseline of the eastern zone (Figure 3) constitutes segments of baseline and a baseline of the western zone is divided into three segments because of the discrete shoreline orientation to draw transects of length 3500 m and 1000 m respectively from the baseline as it intersects all the shorelines. A total of 2516 transects were drawn for the study area at 20 m intervals along the baseline (1194 for the eastern zone, and 1322 for the western Zone) (Figure 3 and Figure 4). This transect interval was small enough comparatively to the satellite image resolution, but transect intervals below than this would not provide a better understanding of shoreline change.

5. Results

An analysis of the shoreline change rate revealed that both zones exhibited both types of change, with erosional trends as negative values and accretional trends as positive values. The eastern zone exhibited accretional action in the shoreline position throughout the study period while the western zone mostly showed both accretion and erosion. Therefore, the results of the present study identified areas with accretion and erosion in the shoreline position at each transect drawn (Table 1 and

Table 2. Maximum and mean actions (accretion and erosion) at eastern and western zones.

	Rate	Eastern Zone			Western Zone
Accretion		Max	Location		Max	Location
	EPR	32	South Asian Terminal		3.88	Haji Ali Goth
	LMS	27	DHA zone D		4.37	Haji Ali Goth
	LRR	35.49	DHA zone D		5.66	Haji Ali Goth
	NSM	2.43	South Asian Terminal		0.30	Haji Ali Goth
Erosion		Max	Location		Max	Location
	EPR	−1.85	South Asian Terminal		−3.96	Manora Beach
	LMS	−14	Shireen Jinnah colony		−3.98	Manora Beach
	LRR	−2.21	South Asian Terminal		−3.27	Manora Beach
	NSM	−0.14	South Asian Terminal		−0.30	Manora Beach
Accretional Transects		Mean	Count	Percent %	Mean	Count	%
	EPR	14.05	1130/1194	94.6	1.37	339/1322	25.6
	LMS	9.99	1084/1194	90.8	0.58	523/1322	40
	LRR	15.31	1130/1194	94.6	0.98	296/1322	22.4
	NSM	1.06	1130/1194	94.6	0.08	338/1322	25.6
Erosional Transects		Mean	Count	Percent %	Mean	Count	%
	EPR	−1.48	64/1194	5.4	−1.15	983/1322	74.4
	LMS	−3.23	111/1194	9.3	−0.85	799/1322	60.4
	LRR	−1.89	64/1194	5.4	−0.99	1026/1322	77.6
	NSM	−0.11	64/1194	5.4	−0.08	984/1322	74.4

Note: EPR, LMS, LRR = Rate of shoreline change (m/year) and NSM = Shoreline movement (Km).

). Table 1 represents the types of shoreline, transect-based change rates from selected metrics (EPR, LMS, LRR, SCE, and NSM) for both zones with respect to their local areas, and a total number of transects drawn for each local area. Table 2 shows maximum accretion and erosion rates recorded for transects in both zones, and the percentage of transects that recorded erosion or accretion.

5.1. Rate of Shoreline Change for the Eastern Zone

The maximum average accretion rate in this zone was observed at DHA Phase 8 Zone C with a rate of 22 m/year at a transect identity (TID) from 191–255 (Table 1). The lowest average rate of accretion was observed at Sea-view Clifton beach with an average rate of 6 m/year (TID 641–897) (Table 1). The lowest to highest EPR rate in the eastern zone, from −1.85 to 32 m/year was observed at South Asian terminal (Table 2). Using the EPR method in the eastern zone, 1130 out of 1194 transects (94.6%) experienced accretion, while 64 (5.4%) exhibited erosion. The mean values of accretion and erosion from EPR were 14 m/year and –1.48 m/year for the total accretional and erosional transects in the eastern zone (Table 2).

The LMS was the second metric used for estimating shoreline change rates. The eastern zone showed a high variation of accretional shoreline change rates with a mean maximum of 17 m/year at DHA Phase8 zone D and a mean lowest change of 2.5 m/year at DHA Phase8 zone B (Table 1). The transects with the greatest erosional LMS rate of –14 m/year and greatest accretional rate of 27 m/year were estimated at Shireen Jinnah Colony and DHA phase8 zone D respectively (Figure 5, Table 2). The LMS rate showed that 90.7% of transects were accretional with an average rate of 10 m/year while 9.3% of them were erosional with a mean rate of –3.23 m/year. Similarly, observed values of the LRR for eastern zone showed that the DHA golf club exhibited a maximum mean accretional rate of 25 m/year and South Asian terminal had the lowest mean rate of 5.6 m/year (Table 1). The greatest significant LRR erosional rate of -2.21m/year and accretional rate of 35 m/year (Table 2) were estimated for the transects of DHA Phase8 zone D and South Asian terminal respectively (Figure 5). Furthermore, the overall LRR records show that 94.6% of transects drawn in eastern zone experienced shoreline accretion, with mean rate of 15.3 m/year while only 5% of transects recorded erosion, with mean rate of -1.89 m/year (Table 2).

The estimation of the shoreline changes in the study area for the eastern zone constitutes on the predicted values of shoreline change envelope (SCE) and net shoreline movement (NSM). The results of SCE reveal that a mean highest accretional SCE value of 1.9 km was obtained for the DHA phase 8 zone D (TID from 256-491) (Table 1). However, it cannot be concluded that the shoreline moved seaward. Similarly, the lowest mean accretional SCE value of 0.57 km was observed for the Sea View Clifton beach (Table 1). The highest and lowest SCE values of 2.4 km and 0.15 km were obtained for the South Asian terminal (TID from 1060-1194) (Table 2). Furthermore, the NSM values recorded in the analysis indicate that the maximum mean accretion distance between the oldest shoreline position of 1942 and the recent shoreline position of 2018 was 1.7 km, obtained for DHA phase 8 zone C (TID from 191–255) (Table 1). Similarly, the lowest mean NSM value of 0.4 km was obtained for Sea View Clifton Beach (TID from 641–897) (Table 1). The highest accretional and erosional values of NSM were obtained for South Asian terminal ranging from 2.43 km to -0.14 km respectively (Table 2, Figure 6). Overall results for NSM shows that 95% of the eastern zone transects experienced accretion with a mean seaward movement of the shoreline of 1.06 km, while 5.4% of transects in the eastern zone show erosion, with a mean retreat of -0.11 km (Table 2) (Figure 6).

5.2. Rate of Shoreline Change for Western Zone

The western Zone (Figure 4) constitutes rocky and sandy beaches, and exhibited both erosional and accretional trends in the shoreline (Figure 7 and Figure 8). Hawks Bay, Haji Ali Goth (French Beach) and Manora Beach are recreational beaches with beach huts and water sports parks with thousands of visitors daily. The maximum mean accretion rate of 0.18 m/year by EPR was observed for the Kaka Pir area, having TID ranges from 634–695 (Table 1, Figure 7). The lowest mean erosional rate of −1.5 m/year by EPR (Table 1) was observed at Jamali Goth area (TID from 964–1040) (Table 1). The highest rate of accretion and erosion shoreline change was recorded at Manora Beach with a value of 3.86 m/year and −3.96 m/year (Table 2). The EPR method shows that 26% of the 1322 transects 1322 experienced accretion, with a mean linear rate of 1.37 m/year (Table 2). Similarly, erosive action was dominant in the western zone, with 1322, or 74% of transects drawn showing erosion, at a mean rate of −1.15 m/year.

LMS values for the western zone showed that the highest mean accretional rate occurred at Kaka Pir area at a rate of of 0.12 m/year (Table 1) and highest mean erosional rate of −1.7 m/year for Jamali Goth (Table 1). The highest accretional rate of shoreline change from the LMS method was observed over Haji Ali Goth, with a rate of 4.37 m/year (Table 2). Similarly, the highest erosional rate of change by LMS was obtained for Manora beach with a value of −3.94 m/year (Table 2) (Figure 7). The overall results for the LMS method indicated that 40% of transects drawn encountered accretion, with mean accretional rate of 0.5 m/year and 60% of them recorded erosion with the mean erosional rate of −0.85 m/year.

Moreover, the values of LRR showed that the highest mean accretional rate of 0.07 m/year was observed for Kaka Pir and the highest mean erosional rate of −1.79 m/year for Jamali Goth in the zone (Table 1). Similarly, the highest accretional and erosional rate by LRR method for each transect drawn shows that Haji Ali Goth had a highest positive rate of 5.66 m/year while Manora beach with a highest erosive rate of −3.27 m/year (Table 2). The LRR showed accretion in 22% of transects in the western zone with mean accretion rate of 0.98 m/year, and erosion in 78% at a rate of −1 m/year (Table 2).

Movement of the shoreline position was estimated by the SCE and NSM methods (Figure 8). The highest mean distance change by SCE shows value of 0.23 km obtained for the Manora beach (TID from 1-342) which shows the highest variability in shoreline position during the study period (Table 1). Similarly, the lowest mean distance change by SCE method shows 0.10 km value for Sandspit beach area (TID from 343-633) (Table 1). The highest variability in shoreline position at any transect drawn in western zone by SCE method was obtained for the Manora beach with a value of 0.769 km and minimum distance variability in shoreline position for Haji Ali Goth observed with a value of 0.016 km (Figure 8) (Table 2). The estimated NSM values show displacement of the shoreline on both landward and seaward sides. The NSM result indicate that highest mean accretion occurred at Haji Ali Goth (TID from 1127–1222) with a value of 0.17 km (Table 1). Similarly, the highest mean landward shifting (regression) of the shoreline encountered at Somar Goth (TID from 1041–1086) with a value of −0.07 km (Table 1). The highest net shoreline shift towards sea (accretion) encountered at Haji Ali Goth with a value of 0.29 km and highest landward shoreline shifting (erosion) was observed at Manora beach with a value of 0.3 km (Figure 8) (Table 2). The NSM results show that 26% of transects drawn in the western zone moved sea, with mean net shoreline change of 0.08 km, while 74% of transects showing landward movement, with a mean retreat of −0.08 km (Table 2).

6. Discussion

This study examines variability in historical shoreline position along the coastal city of Karachi, using historical topographic maps (1942 & 1966), Landsat satellite imagery (1976, 1990, 2002, and 2011) and Planet Scope satellite imagery (2018). Historical shoreline positions were extracted using NDWI and shoreline changes were observed using the DSAS approach. The study provides in-depth fine spatio-temporal shoreline change analysis and contributes to exploring the causes of such changes. Results show significant accretionary change in shoreline position in the eastern zone (Figure 3) (Table 2), where large residential projects have encroached on the sea in the past decade. On the other hand, the western zone (Figure 4) (Table 2) has the highly stable shoreline. The results also demonstrate that the shoreline changes in the western zone were related to the changes in hydrodynamics due to land reclamation projects in the eastern zone, associated with the local longshore current. Similar findings have been reported by recent studies conducted in other parts of the Indus delta region where large land reclamation projects have increased coastal vulnerability by altering the ocean hydrodynamics. For instance, Waqas et al. (2019) [16] found that 38% of the barrier islands in the Indus delta were vulnerable to oceanic factors. Similarly, Salik et al. (2015) [41] reported the vulnerability of the communities living in the lower Indus deltaic plain due to climatic and anthropogenic impacts on the coast.

6.1. Qualitative Rate of Shoreline Change

The shoreline positions in both zones have experienced accretion and erosion in different ways. The length of shorelines in the eastern and western zones have increased by 2.6% and 16.6% respectively, over a period of 76 years. The maximum accretional and erosional rates of shoreline recorded at each transects by EPR and NSM methods help to quantify the length of shoreline vulnerable to extreme ocean events (

Table 3. Shoreline classification based on EPR [39].

Rate of Shoreline Change (x m/y)	Shoreline Classification	Affected Shoreline of Easter Zone (%)	Affected Shoreline of Western zone (%)
x ≤ −2	Very High Erosion	0.0	9.9
−1 < x ≤ −2	High Erosion	5.1	27.1
0 < x < −1	Moderate Erosion	0.3	37.4
x = 0	Stable Shoreline	0	0
0 ≤ x < 1	Moderate Accretion	0.0	13.7
1≤ x <2	High Accretion	0.8	8.4
x ≥ 2	Very High Accretion	93.8	3.6

Note: Affected shoreline length reflects the percentage (%) of total shoreline length of each zone.

) while SCE can help to track the overall movement in the shoreline position. This study revealed that 5.4% of transects drawn in the eastern zone (25 km), representing 1.28 km of the whole shoreline length, experienced erosion (Table 2) and of this, 1.22 km at South Asian terminal underwent high erosion rates (Table 3). Similarly, 94.6% of transects drawn in the eastern zone, corresponding to a 22.6 km shoreline length, showed accretion (Table 2), of which 22.4 km exhibited very high accretion rates (Table 3) at DHA Zone B, DHA Golf Club, DHA Zone C, DHA Extension, Sea-view, Shireen Jinnah Colony, and South Asian Terminal. The DHA Zone B, DHA Golf Club, DHA Zone C, and DHA Extension corresponding to a combined shoreline length of 12.8 km showed high accretional rates (Table 3) of shoreline change because of massive land reclamation projects for urbanization and industrialization (Figure 9). Similarly, Sea-view is a recreational beach with a low substrate slope and had high accretion rates (Table 3). Table 3 indicates that 4.94 km of the beach length was nourished by trapping, dredging, and constructional materials from offshore drift. This beach regularly receives large numbers of visitors from all over Pakistan (Figure 9). Furthermore, Shireen Jinnah Colony had high accretion rates (Table 3) due to soft sediments deposited by longshore drift after construction of breakwaters preventing seawater and high energy waves from entering the residential areas (Figure 9) (Section 6.2).

The study also revealed that during this period, 74.4% of transects drawn in the western zone (29 km) corresponding to 19.7 km of the total coastline experienced erosion (Table 2), and of this, 9.9 km experienced moderate erosion (Table 3) at Manora Beach, Sandspit, Kaka Pir, Hawks Bay, Jamali Goth, Somar Goth, Abdul Rehman Goth, Haji Ali Goth, and Engro Beach Huts (Figure 9). A total of 2.5 km out of the 5.3 km shoreline length of Hawks Bay and 4 km of the 5.82 km of Sandspit shows moderate erosion action (Table 1 and Table 3) because they have almost flat substrate slopes that are being eroded by waves and current action during the south-west monsoon season (Section 6.2). Similarly, 1.82 km out of the shoreline length of the 6.64 km western part of Manora beach and 0.7 km out of the 1.54 km of Jamali Goth showed severe erosion from strong wind-driven longshore currents in the western part of study area (Figure 9). Furthermore, 25.6% of transects drawn in the western zone corresponding to 6.8 km of the shoreline experienced accretion, of which 3.62 km showed moderate accretion (Table 3). These involved some parts of Manora beach, Hawks Bay, Kaka Pir, and Haji Ali Goth. High to very high accretion was encountered at the northeastern part of Manora beach due to soft sediments deposited by longshore drift from the construction of the deep-water container port, corresponding to 1.38 km to 0.88 km of shoreline (Figure 9).

6.2. Hydrodynamics of the Sindh Coastal Zone

Pakistan’s coastal hydrodynamics are controlled by the reversal of the monsoonal system. Significant oceanic factors such as wind and wave climatology influence the coastal environment of Karachi because these changes are associated with the local wind driven longshore current [16,51,66,67]. To analyze the wind and wave climate influencing the study area seasonally before and after the major encroachment by urban development (Figure 10 and Figure 11), the long-term offshore daily-hourly data from the National Oceanic and Atmospheric Administration (NOAA) was divided into two periods: Pre-encroachment (2008–2010) and post-encroachment (2011–2017). During the pre-summer monsoon season (March–May) of the pre-encroachment period, the wave height (WH) results showed an absence of swell for 35% of the time, while waves exceeding the 2 m height were recorded for only about 2% of the whole data (Figure 12). During this season, winds were northwesterly, with an average wind speed (WS) of 4.4 knots having an average wave period (WP) of 13 s, and only 0.5% of the time did WS exceed 9 knots, generating large wavelets (Figure 13). However, during the summer (Jun–Sep) i.e., the southwest (SW), monsoon season winds were southwesterly (Figure 10) with an average WS of 9 knots, generating strong waves with a mean significant WH of 3 m, and an average WP of 11 s.

Surprisingly, 22% of the time WH exceeded 4 m, and occasionally reached 6 m due to a rising tide, generated by heavy SW winds of WS above 12 knots (Figure 13). During the SW monsoon, high-energy southwesterly waves strikes the coast, accelerating erosive action in the western zone and sediment deposition in the eastern zone from longshore drift. During the retreat of the monsoon season (Oct–Nov), northeasterly winds with an average WS of 5 knots generated waves with an average WH lower than 1 m and a mean WP of 11 s. During this season, 86% of the observed time WS was observed to be between 3–9 knots, which generated small waves with a significant WH lower than 2 m (Figure 12). Furthermore, during the northeastern (NE) monsoon transition season (Dec–Feb), northeasterly winds (Figure 10) with a low mean WS of 5 knots generated small wavelets of an average 1 m wave height, of 1 m, which allowed the settling of suspended sediments. During this season, the NE wind remained lower than 9 knots for 99% of the time (Figure 13) and WH remained less than 2 m for 98% of the time (Figure 12).

During the pre-summer monsoon (NW) season (Mar–May) of the post-encroachment (2011–2017) period (Figure 11), the average WS was more than 4 knots and only 14% of the observed time the WS remained between 6–9 knots (Figure 12). The persistent NW winds generated a small swell of which about 97% of the swell’s significant WH remained below 2 m (Figure 12) with an average WP of about 13 s. Similarly, average WS of about 9 knots during the SW monsoon season (Jun–Sep) (Figure 11) generated strong wave swells with a mean wave height of more than 3 m having a mean WP of 11 s. During the SW season, 26.6% of the time WH exceeded more than 4 m generated by WS reached to 22 knots (Figure 13) and a maximum 7.7 m of WH recorded during the post-encroachment period (Figure 12). During the retreating period of the monsoon season (Oct–Nov), wind prevailed from the NE side with an average WS of about 5 knots, which generated a swell wave of an average significant wave height of 1 m. The average WP of 13 s during the season showed the calmness of the sea where a total 98% of the time WH recorded below 2m, which helps in the settlements of the sediments in the bottom (Figure 12). During the winter monsoon season (NE) (Dec–Feb), winds with an average speed of more than 5 knots generated swells with an average significant WH of 1 m. Consistent and normal wind generated 98% of swell wave less than 3 m (Figure 12) with an average WP of about 11 s.

Generally, the wave height and its energy in the study area directly or indirectly subject of the wind speed, wave period, surface area over which the wind blows, and bathymetry of the area. The wave packets remain in action until the energy of the waves is absorbed by coastal region. Therefore, shoreline vulnerability of Karachi will probably be accelerated by high wave-height generated by high-speed winds at a high frequency (Figure 12 and Figure 13). The pre- and post-hydrodynamic study revealed that about a 2.95% average WH increased significantly during the SW summer monsoon season and a 0.97% average for the NW winter monsoon season indicating the vulnerability of the shoreline position in the future. The sediments usually deposited during the stable ocean condition when the energy of waves is low and the deposited soft sediments in the shelves are driven away from the coast by the action of longshore drift and high wave energy [43]. The wave and wind are important factors for changing the nourishment and erosion of the coastal environment and development of urban areas along the Sindh coastal zone [15,68,69].

7. Conclusions

This study presented tracks changes along a coastal mega city Karachi, Pakistan by employing historical topographic sheets dated back to 1942 and images from different missions of Landsat and Cubesats, until 2018. The total shoreline length was divided in two eastern zone and western zones. The mean accretional End Point Rate (EPR) for the eastern zone was observed to be 14 m/year, which indicated that 95% of the total transects of the eastern zone encountered accretion. This implies that the eastern zone remained under extensive land reclamation activities. Similarly, the EPR rate of change for the western zone revealed a mean accretional rate of 1.37 m/year, and with that 26% of the total transects underwent positive change. This indicates a small positive change rate for a quarter of the study area, while 74% of the transects in the western zone experienced mean an erosional EPR value of −1.15 m/year, indicating that a major proportion of the study area faced erosional activities. Furthermore, a total of 23.5 km of the shoreline length of eastern zone was designated as very high accretional while 2.9 km of western zone shoreline was found prone to very high erosional activity. The pre- and post-land reclamation hydrodynamics revealed that an average wave height increased by 2.95% during the south-west summer monsoon season and 0.97% during the NW winter monsoon season, indicating potential future threats to the shoreline position. Most of the western part of the shoreline is still undisturbed by human intervention. The shoreline position at a place is mostly determined by sand supply for nourishment, geomorphology, and oceanic forcing factors coupled with sea level change. Coastal land development via encroaching sea is a two-way process and triggers disturbance in the ocean hydrodynamics, which may create future stress on communities living in low-lying coastal areas. There are potential chances of land subsidence on reclaimed coastal areas following extreme events. Our study shows the potential of GIS and remote sensing techniques for comprehensive coastal risk assessment and our method would be very useful for the entire delta or other deltaic areas. Furthermore, it was observed that low to medium spatial resolution data sets (i.e., topographic sheets, Landsat MSS, TM, ETM+, and OLI) have higher uncertainty than the higher spatial resolution (3 m) Planet Scope imagery. The availability of Planet Scope imagery for educational and research purpose has opened new doors for the exploration of new techniques to preserve the environment. It is recommended that proper coastal protection and coastal management along ocean facing coastal areas should be adopted to defend against the erosive action of the ocean. The master plan for Karachi coastal areas should be revised to be friendlier to the coastal environment, with more emphasis on sustainable coastal development and coastal safety. In future, we will employ coastal vulnerability tools to quantify physical and socio-economical vulnerabilities of the Sindh coastal region to examine the effect of seawater intrusion due to sea level rise.