Hazard Potential in Southern Pakistan: A Study on the Subsidence and Neotectonics of Karachi and Surrounding Areas

: Coastal communities in deltaic regions worldwide are subject to subsidence through a combination of natural and anthropogenic processes. The city of Karachi in southern Pakistan is situated along the diffuse western boundary of the tectonically active Indian Plate, making it more susceptible to natural subsidence processes from plate motion-related deformational events such as earthquakes and faulting. Karachi has a dense population of over 16 million people, and determining the rate of subsidence and extent of neotectonic activity is crucial for mitigating seismic hazards. Excessive abstraction of groundwater and extensive groundwater use in irrigation are some of the anthropogenic contributions to subsidence in the area. A combination of the lack of historical data and few previous studies of the area make it difﬁcult to determine the rate and extent of deformation in this region. We present an analysis of subsidence and neotectonic activity in Karachi and its surrounding areas using Interferometric Synthetic Aperture Radar (InSAR) timeseries techniques. The InSAR results for satellite LOS velocity change in both ascending and descending Sentinel-1 tracks indicate subsidence in key residential and industrial areas. Further decomposition into two dimensions (east–west and vertical) quantiﬁes subsidence in these areas up to 1.7 cm per year. Furthermore, InSAR data suggest the presence of an active north–east dipping listric normal fault in North Karachi that is conﬁrmed in the shallow subsurface by a 2D seismic line. Subsidence is known to cause the reactivation of faults, which increases the risk of damage to infrastructure


Introduction
Coastal communities in deltaic regions worldwide are subject to subsidence through a combination of natural and anthropogenic processes [1][2][3][4]. Natural processes include the dewatering/compaction of sediments [2,5,6], faulting, and load-induced crustal downwarping [2]. These processes cause slow and steady deformation over time due to their interaction with other natural processes that offset the subsidence rate, such as sedimentation and the production of new soils from organic decay [2]. Conversely, anthropogenic processes that cause subsidence, such as the extraction of groundwater leading to the compaction of aquifers and induced fault motion through mining or fluid extraction, tend to cause relatively faster and temporally variant rates of subsidence [2,3,7].
A significant anthropogenic cause of subsidence is excessive groundwater abstraction. Fluid extraction rates outpacing recharge leads to compaction within the aquifer [3,[8][9][10]. Combining these natural and anthropogenic processes significantly increases the risk of subsidence in densely populated coastal areas.
The city of Karachi in southeast Pakistan is one such coastal city in a deltaic environment. Situated along the diffuse western boundary of the tectonically active Indian Plate, it is more susceptible to natural subsidence processes from plate motion-related deformational events such as earthquakes and faulting. West of the city, the triple junction between the Indian, Eurasian, and Arabian plates lies just off the coast of Pakistan. This kind of deltaic setting near a tectonically active junction is also found at the San Francisco Bay area and Los Angeles, California, where the seismic hazard is characterized as high [3,11]. As such, it can be argued that the lack of existing data and research conducted in Karachi leaves this area in an extremely vulnerable position and susceptible to accelerating subsidence rates.
Due to the lack of regulatory legislation in place for the pumping of groundwater, extensive groundwater use in irrigation, and the Indus Basin being one of the most stressed aquifers in the world [12], water shortages are becoming increasingly frequent in Karachi. With all these factors in consideration, it would be safe to assume that the city of Karachi is on track to follow the same subsidence trends as the Houston-Galveston area before the establishment of the Groundwater Management and Subsidence Districts [8,13,14]. Understanding spatial and temporal variations in deformation may be the key for assessment of the current and potential future hazards related to subsidence and tectonics.
The mobile fold and thrust belt that underlies the city and extends out towards Hyderabad, termed the 'Karachi Arc' [15], is considered to be currently active and is suggested to have been for a long time. The earliest documented major earthquake (8 or greater on the Modified Mercalli Scale) occurred in 893 AD and is believed to be responsible for the destruction of the port city of Debal, which was located near modern-day Karachi [11,15,16]. Although there have been reports of shaking from earthquakes in the city over the past few decades, the city has never experienced severe damage due to an earthquake. Despite the evidence of recurrent low to moderate grade seismic activity [15], the time scale required to classify the seismic cycle of major earthquakes varies from 80 to 1000 years depending on the specific area within the region [11,16]. The question posed by past research [11], whether Karachi happens to be in an aseismic zone or whether it is simply a lack of information on the city's history that would suggest so, remains unanswered, but there is evidence of nearby towns destroyed by earthquakes in the past [11,17]. Considering the given evidence and the lack of research conducted in the area, it is imperative that an innovative approach to risk assessment be adopted for the region.
We present an analysis of the subsidence and neotectonic activity of Karachi and its surrounding areas ( Figure 1) using Interferometric Synthetic Aperture Radar (InSAR) timeseries analysis techniques. We estimate the rate of subsidence to be about 1.7 cm per year in key areas of Karachi such as Defense, Korangi, and North Karachi. For the first time, we present evidence of an active fault in the city using a new full resolution InSAR approach along with inverse fault modeling.

Geological Background
The city of Karachi is located at the southern end of the Kirthar Mountain Rangea fold and thrust belt that formed as a result of the Late Cretaceous collision between the Indian and Eurasian plates [18][19][20]. This range marks the diffuse western boundary of the Indian plate and has predominantly undergone transpressional deformation since its initial collision with the Afghan block of the Eurasian plate [18,21,22] (Figure 1). The city itself is built on what is known as the Karachi Arc, a large fold and thrust belt which shows thin-skinned deformation and eastward movement [15,18,23,24], with the detachment surface being located in the Goru shale sequences, about 2-4 km below sea level [15,18,23]. Northeast of the city, the Karachi arc is characterized by asymmetric, doubly plunging anticline-syncline pairs trending north-south in parallel to en echelon arrangement [18] (Figure 2). Historical seismicity data suggests that the Karachi Arc is still active [15,16]. The eastward creep may be due to active subsidence in the underlying Hyderabad graben [15,25]. Northeast of the city, the Karachi arc is characterized by asymmetric, doubly plunging anticline-syncline pairs trending north-south in parallel to en echelon arrangement [18] ( Figure 2). Historical seismicity data suggests that the Karachi Arc is still active [15,16]. The eastward creep may be due to active subsidence in the underlying Hyderabad graben [15,25].

Figure 2.
A map of the study area depicting key tectonic features such as anticlinesyncline pairs, faults, and recent earthquake hypocenters, along with InSAR coverage footprints. Earthquakes data source: Pakistan Meteorological Department Karachi, as presented by [14]. Faults and folds traced from Sindh geological map published by Geological Survey of Pakistan.

Methods
We use the Interferometric synthetic aperture radar Scientific Computing Environment (ISCE) developed by NASA's JPL to process 276 descending track and 134 ascending track Sentinel-1 Single Look Complex (SLC) images available from 2015 to 2020 with a 1 Figure 2. A map of the study area depicting key tectonic features such as anticlinesyncline pairs, faults, and recent earthquake hypocenters, along with InSAR coverage footprints. Earthquakes data source: Pakistan Meteorological Department Karachi, as presented by [14]. Faults and folds traced from Sindh geological map published by Geological Survey of Pakistan.

Methods
We use the Interferometric synthetic aperture radar Scientific Computing Environment (ISCE) developed by NASA's JPL to process 276 descending track and 134 ascending track Sentinel-1 Single Look Complex (SLC) images available from 2015 to 2020 with a 1 arc-second SRTM digital elevation model into a stack of coregistered SLCs. This is accomplished by using the Sentinel-1 Tops Stack Processor [26]. The stack of coregistered SLCs is then used as the input for the Miami Phase Linking in Python (MiaplPy) software to perform full resolution InSAR processing using non-linear phase inversion [27,28]. The resulting files are then imported and processed using MintPy [29] for phase unwrapping error corrections and timeseries analysis. A simplified diagram of the InSAR workflow can be seen in (Figure 3). is then used as the input for the Miami Phase Linking in Python (MiaplPy) sof perform full resolution InSAR processing using non-linear phase inversion [27, resulting files are then imported and processed using MintPy [29] for phase unw error corrections and timeseries analysis. A simplified diagram of the InSAR w can be seen in (Figure 3).

ISCE
We developed the stack using the 10 nearest neighbor connections for each tion and used the Network-based Enhanced Spectral Diversity (NESD) approach form coregistration [26].

Miaplpy
The SLC and geometry files output by ISCE are loaded into MiaplPy at ful tion. The SLC stack is then subdivided into patches (200 × 200 pixels) and proc parallel using a Beowulf cluster with 16 workers. The patches are organized in stacks of 10 images each for phase linking. A full network non-linear phase linkin formed by the sequential Eigenvalue decomposition-based Maximum Likelihoo terferometric phase (EMI) method to estimate the wrapped phase timeseries [30,3 the following equation [27,31]: where ˆ is the estimated complex coherence matrix and ∘ represents the Ha product. The desired solution for this equation is the eigenvector ( ) that corresp the minimum eigenvalue ( ) [27]. After the patches are concatenated, interferograms are generated using a si erence and unwrapped using the Statistical Cost, Network Flow Algorithm for Ph wrapping (SNAPHU) [32]. The unwrapped interferograms are loaded into MintP rect for phase unwrapping errors before the network of interferograms is inverte timeseries of deformation. We use MiaplPy to perform least square inversion of wrapped interferogram network and convert the estimated unwrapped timese range-change timeseries. We then perform timeseries error corrections and geo final results in MintPy.

Mintpy
Once the stack is loaded into MintPy, a reference pixel is chosen in the far deformation with high temporal coherence (≥0.85). Temporal coherence ( ) is ca to assess the quality of each pixel in the raw phase timeseries [29] using the initi value ( ) and the estimated phase value ( ) in the following equation [27,33]

ISCE
We developed the stack using the 10 nearest neighbor connections for each acquisition and used the Network-based Enhanced Spectral Diversity (NESD) approach to perform coregistration [26].

Miaplpy
The SLC and geometry files output by ISCE are loaded into MiaplPy at full resolution. The SLC stack is then subdivided into patches (200 × 200 pixels) and processed in parallel using a Beowulf cluster with 16 workers. The patches are organized into mini-stacks of 10 images each for phase linking. A full network non-linear phase linking is performed by the sequential Eigenvalue decomposition-based Maximum Likelihood of Interferometric phase (EMI) method to estimate the wrapped phase timeseries [30,31] using the following equation [27,31]: whereΓ is the estimated complex coherence matrix and • represents the Hadamard product. The desired solution for this equation is the eigenvector (ν) that corresponds to the minimum eigenvalue (ζ min ) [27]. After the patches are concatenated, interferograms are generated using a single reference and unwrapped using the Statistical Cost, Network Flow Algorithm for Phase Unwrapping (SNAPHU) [32]. The unwrapped interferograms are loaded into MintPy to correct for phase unwrapping errors before the network of interferograms is inverted into a timeseries of deformation. We use MiaplPy to perform least square inversion of the unwrapped interferogram network and convert the estimated unwrapped timeseries to a range-change timeseries. We then perform timeseries error corrections and geocode the final results in MintPy.

Mintpy
Once the stack is loaded into MintPy, a reference pixel is chosen in the far field of deformation with high temporal coherence (≥0.85). Temporal coherence (γ) is calculated to assess the quality of each pixel in the raw phase timeseries [29] using the initial phase value (θ n ) and the estimated phase value (ϕ n ) in the following equation [27,33]: where N is the number of SAR images in the stack, i represents individual images, and nk is a wrapped phase interferogram generated using images acquired at time n and k. Temporal coherence is used as a reliability measure or statistical "goodness of fit" [27,33]. Unwrapping errors are then determined and corrected using the bridging method [29]. After network inversion in MiaplPy, we apply a temporal coherence threshold of 0.7 to mask out anomalous pixels. A quadratic phase ramp is estimated and removed from the reliable pixels in the data in order to correct for residual tropospheric and ionospheric delays [29]. This step is followed by correction of the topographic phase residual [29].
To improve the quality of the data, the estimated residual phase is used to determine noise in the timeseries. The root mean square error is calculated on the residual phase and used to identify and remove noisy acquisitions through the following equation from [29]: N], Ω represents the reliable pixels chosen from temporal coherence masking,φ i resid represents the residual phase at time i, and λ represents the radar wavelength. The median absolute deviation is calculated and SAR acquisitions with an RMS higher than three median absolute deviations are considered noisy and excluded from further processing. Finally, the average velocity is estimated from the timeseries to determine the rate of deformation using the following equation [29]: where v is the velocity, t i represents the time at the i th acquisition, c is an unknown offset constant, andφ i disp is the displacement timeseries. The ascending and descending track velocities are resampled and subset to the same spatial resolution and coverage. The resampled datasets are then used as inputs for decomposition in two dimensions (horizontal and vertical) using Mintpy. The LOS velocity can be decomposed into three individual directional components of displacement. These are north-south (v n ), east-west (v e ), and vertical (v u ) [34,35]: where v los is the line-of-sight displacement velocity, α is the azimuth angle of the satellite heading, θ is the radar incidence angle at the ground surface, and δ los is measurement error due to imprecise satellite orbit geometry, tropospheric delay, topographic phase contribution, etc. This equation can be applied to both ascending track (v asc los ) and descending track (v desc los ) results to form an equation with three unknowns for each track: v asc los = (sin θ sin α) asc v n − (sin θ cos α) asc v e + cos θ asc v u + δ los Assuming negligible contributions from both the measurement error δ los and northsouth displacement component v n , Equations (6) and (7) can be used to determine the displacements in the horizontal (east-west) and vertical directions using Equations (8) and (9), respectively [36]: v e = cos θ asc v desc los − cos θ desc v asc los sin θ asc cos α asc cos θ desc − cos θ asc sin θ desc cos α desc (8) v u = sin θ desc cos α desc v asc los − sin θ asc cos α asc v desc los cos θ asc sin θ desc cos α desc − sin θ asc cos α asc cos θ desc (9)

Results
The deformation time series and LOS velocity change in both ascending and descending tracks (Figure 4) indicate movement away from the satellite in the following key areas: Defence Housing Authority (DHA), Korangi, Landhi, Leely Town, and North Karachi. There are apparent uplift regions due to the selection of a reference point with high temporal coherence in the far field of deformation. While it is possible to select a reference point in the areas that appear to be uplifting, these areas were determined to have a slightly lower temporal coherence, and they are in proximity to actively deforming regions. As such, the current reference point (denoted on the figure with a black star) was selected for reliability and stability. There are apparent uplift regions due to the selection of a reference point with high temporal coherence in the far field of deformation. While it is possible to select a reference point in the areas that appear to be uplifting, these areas were determined to have a slightly lower temporal coherence, and they are in proximity to actively deforming regions. As such, the current reference point (denoted on the figure with a black star) was selected for reliability and stability.
The LOS velocities indicate a similar trend for these focus areas. North Karachi has a displacement rate of between −0.7 to −0.8 cm/yr and displays a sharp transition boundary between no displacement and movement away from the satellite. While there is no apparent pattern for DHA, the areas of Korangi, Landhi, and Leely Town all show a bullseye pattern of no displacement at the edges and increasing displacement inward with maximum displacement at the center of between −0.7 to −1.8 cm/yr. Further decomposition into two dimensions (east-west and vertical) quantifies the vertical component of displacement in these areas to about 1.7 cm downward per year ( Figure 5). Timeseries of individual pixels representing the average displacement in each of the key areas ( Figure 6) show consistent displacement from 2015 to 2020.

Natural Causes
Proximity to the Indus Delta renders Karachi susceptible to the effects of deltaic subsidence. The eastern Indus delta is influenced by footwall subsidence associated with thrust faulting in the Kachchh mainland, as well as recurrent tectonic activity in the Rann of Kachchh [37,38].
Subsidence in the north may be influenced by neotectonic activity. The line of sight velocity results (Figure 4) and the vertical displacement results ( Figure 5) display a sharp boundary between downward motion and no displacement. This boundary is aligned with the inferred trace of the Allah Bund fault, which is estimated to have a 50-70°N dip with listric normal fault geometry [39]. The inferred fault trace is traced using the geological map of Sindh published by the Geological Survey of Pakistan. Although the displacement boundary is clear in the northwestern portion of the fault, it diminishes to the southeast. This may be due to a change in the fault dip or the fault terminating in the area. It is more likely that the signal is obscured by proximity to the river and the strong subsidence signals in the southeast compared to the relatively stable northwest. Recent earthquakes clustered in the hanging wall block suggest strain along the fault may be contributing to recurrent low-to moderate-grade seismic activity. Although there is a strong indication

Natural Causes
Proximity to the Indus Delta renders Karachi susceptible to the effects of deltaic subsidence. The eastern Indus delta is influenced by footwall subsidence associated with thrust faulting in the Kachchh mainland, as well as recurrent tectonic activity in the Rann of Kachchh [37,38].
Subsidence in the north may be influenced by neotectonic activity. The line of sight velocity results (Figure 4) and the vertical displacement results ( Figure 5) display a sharp boundary between downward motion and no displacement. This boundary is aligned with the inferred trace of the Allah Bund fault, which is estimated to have a 50-70 • N dip with listric normal fault geometry [39]. The inferred fault trace is traced using the geological map of Sindh published by the Geological Survey of Pakistan. Although the displacement boundary is clear in the northwestern portion of the fault, it diminishes to the southeast. This may be due to a change in the fault dip or the fault terminating in the area. It is more likely that the signal is obscured by proximity to the river and the strong subsidence signals in the southeast compared to the relatively stable northwest. Recent earthquakes clustered in the hanging wall block suggest strain along the fault may be contributing to recurrent low-to moderate-grade seismic activity. Although there is a strong indication of neotectonic activity in the area, the fault parameters are poorly constrained. The geological map of Sindh published by the Geological Survey of Pakistan extends the inferred trace of the Allah Bund fault through this area. However, this is not the case for the published geological maps of the Karachi area. There is also conflicting information in the literature.
According to Sarwar (2004), the Allah Bund fault is further south and characterized as a south-dipping thrust fault that terminates within the Rann of Kutchh. There is also an unknown dextral strike-slip fault with east-west orientation indicated that may coincide with this neotectonic activity (see [25], Figure 1). According to the structural analysis of the Mangho Pir anticline in Niamatullah and Imran (2012), the fault may be an extension of, or splay off, the Pir Mangho fault (see [40], Figures 6 and 11). Sarwar and Alizai (2013) further expand on the unknown fault in Sarwar (2004), and suggest it is a wrench fault that extends through the city and terminates at the western coast (see [15], Figure 6). The existence of conflicting information makes it difficult to constrain the fault parameters. However, the abundance of references suggests that it is highly likely there is some neotectonic activity occurring in the area. Analysis of a perpendicular seismic line from the Union Texas 1999 2D seismic reflection survey of Karachi shows clear displacement between reflectors and onlapping bedding planes consistent with listric normal fault geometry (Figure 8). We use the Geodetic Baysian Inversion Software (GBIS) [41] to model the fault parameters using InSAR LOS velocity data from both tracks. Optimal fault model parameters determined a strike of 278° and a slip at depth of about 9.9 cm over 1,000,000 trials. The model results can be Analysis of a perpendicular seismic line from the Union Texas 1999 2D seismic reflection survey of Karachi shows clear displacement between reflectors and onlapping bedding planes consistent with listric normal fault geometry (Figure 8). We use the Geodetic Baysian Inversion Software (GBIS) [41] to model the fault parameters using InSAR LOS velocity data from both tracks. Optimal fault model parameters determined a strike of 278 • and a slip at depth of about 9.9 cm over 1,000,000 trials. The model results can be seen in Table 1.  Fluvial processes may be contributing to subsidence in the central and southern regions of Karachi. In the case of the Korangi and Defence areas, these processes can be attributed to interaction with the Malir River. Subsurface compaction may be enhanced in Fluvial processes may be contributing to subsidence in the central and southern regions of Karachi. In the case of the Korangi and Defence areas, these processes can be attributed to interaction with the Malir River. Subsurface compaction may be enhanced in these areas due to the dewatering of sediments and fluvial erosion. The Defence area is not only subject to similar effects at the outlet of the Malir River but is also exposed to high energy wave action [42].
Timeseries plots of individual pixels representative of average displacement in each of the key areas display a relatively consistent rate of subsidence from 2015 to 2020 ( Figure 5). The highest variations are found during the monsoon season between the months of June and September each year. Historically, this is the time of year when Karachi experiences severe flooding. There have been varying degrees of flood severity between 2017 and 2019, but the city experienced its worst floods in August 2020 since precipitation data collection began in 1931 [43].

Anthropogenic Causes
The construction of dams and barrages upstream has had a significant impact on the geomorphology of the Indus delta by severely limiting the supply of water and sediment downstream, making it a sediment-starved delta. As a result, the now primarily tidedominated delta front is eroding at a faster rate and experiencing further reworking of abandoned streams into intruding tidal channels [38,42].
Localized subsidence within areas of the city can be attributed to multiple causes. One of the primary potential reasons is excessive abstraction of groundwater fueled by rapid urbanization and a growing demand for water in domestic and industrial applications. A steadily increasing shortage in the municipal water supply has encouraged widespread reliance on groundwater abstraction through individual and community wells [44].
There are thousands of suction pumps and tube wells installed at private residences throughout the city, especially in areas like Defence and Nazimabad. However, there is no centralized database with records of these wells to determine a direct correlation between the locations of groundwater wells and rates of subsidence. A regional overview of groundwater in the area can be assembled using the Monthly Mass Grids of Water Equivalent Thickness produced by NASA's Gravity Recovery & Climate Experiment (GRACE) [45][46][47].
Although the data has a coarse spatial resolution of 300 km, the groundwater in the Karachi region has been declining since 2002 ( Figure 9). Recent studies have suggested that this decline is attributed to excessive groundwater abstraction for agricultural, industrial, and residential use due to the inadequate surface water supply from Haleji Lake, Kheenjar Lake, and Hub Dam [24,48,49].
Furthermore, a large portion of Defence is built on land reclaimed between 1984 and 2021 ( Figure 10).

Risk to Infrastructure
Buildings and infrastructure most at risk of damage are not necessarily located in areas where subsidence is greatest. The risk is highest where there is a steep gradient in the deformation velocity [50], as differential subsidence is known to cause surficial faulting and fracturing [51]. We calculate an InSAR deformation gradient to highlight these areas in Karachi (Figure 11). The vertical displacement data is interpolated using an Inverse Distance Weighted approach and resampled to a 100 m resolution continuous raster. A sobel operator is applied to the raster and a 3 × 3 kernel is used to smooth the data. The resulting raster highlights those areas most at risk for damage.
A major cause for concern is that, in some cases, the subsiding areas may be caused in part by damaged infrastructure itself. Specifically, broken sewer and water lines can be attributed to variations in localized subsidence patterns where leaking water may cause liquefaction and subsidence over time [52,53]. This suggests that the areas where differential subsidence is greatest would have the problem further exacerbated by damage to sewer and water systems.     areas where subsidence is greatest. The risk is highest where there is a steep gradient in the deformation velocity [50], as differential subsidence is known to cause surficial faulting and fracturing [51]. We calculate an InSAR deformation gradient to highlight these areas in Karachi ( Figure 11). The vertical displacement data is interpolated using an Inverse Distance Weighted approach and resampled to a 100 m resolution continuous raster. A sobel operator is applied to the raster and a 3 × 3 kernel is used to smooth the data. The resulting raster highlights those areas most at risk for damage. Figure 11. InSAR gradient created using vertical displacement data.
A major cause for concern is that, in some cases, the subsiding areas may be caused in part by damaged infrastructure itself. Specifically, broken sewer and water lines can be attributed to variations in localized subsidence patterns where leaking water may cause liquefaction and subsidence over time [52,53]. This suggests that the areas where differential subsidence is greatest would have the problem further exacerbated by damage to sewer and water systems.
The high InSAR gradient areas in Defence Phase VIII are concentrated along the historical landmass extents ( Figure 9). It has been suggested by [52,54] that this is especially alarming considering the extensive construction of high rise buildings and housing communities on reclaimed land, and that further unmitigated development could have catastrophic results [52].
Considering the density of Karachi's population as well as the density of buildings, there is an alarming potential for catastrophe. As seen in cases like the Brownwood subdivision in the greater Houston area, unmitigated subsidence in coastal cities can lead to the sinking of entire neighborhoods [55]. The InSAR gradient map should help to serve as Figure 11. InSAR gradient created using vertical displacement data.
The high InSAR gradient areas in Defence Phase VIII are concentrated along the historical landmass extents (Figure 9). It has been suggested by [52,54] that this is especially alarming considering the extensive construction of high rise buildings and housing communities on reclaimed land, and that further unmitigated development could have catastrophic results [52].
Considering the density of Karachi's population as well as the density of buildings, there is an alarming potential for catastrophe. As seen in cases like the Brownwood subdivision in the greater Houston area, unmitigated subsidence in coastal cities can lead to the sinking of entire neighborhoods [55]. The InSAR gradient map should help to serve as a guide for city planners to determine which areas require the most attention to address based on the amount of risk.

Limitations and Considerations
There are a few areas in Karachi where subsidence may be occurring but which are not explored in detail in this work. Areas of note include Defence Phase 8, Machar Colony, and Sector 5. DHA Phase 8 and Machar Colony are masked out due to low temporal coherence, likely caused by rapid and sporadic construction. Sector 5 has areas masked out due to the high concentration of vegetation. However, in each of these cases, the same gradient or pattern of deformation as the focus areas of this study can be seen partially formed. Specifically, lower rates of deformation are seen on the outsides, and the few reliable pixels in these areas show increasing deformation rates moving inward. These patterns are apparent in the unmasked data ( Figure 12). However, as temporal coherence is used as a reliability measure, these areas have been masked out but should be taken into consideration in future work. due to the high concentration of vegetation. However, in each of these cases, the same gradient or pattern of deformation as the focus areas of this study can be seen partially formed. Specifically, lower rates of deformation are seen on the outsides, and the few reliable pixels in these areas show increasing deformation rates moving inward. These patterns are apparent in the unmasked data ( Figure 12). However, as temporal coherence is used as a reliability measure, these areas have been masked out but should be taken into consideration in future work.  Moving forward, we recommend the establishment of a dense network of GPS stations encompassing the city. This would help develop a baseline reference for the ground surface and enable more detailed and precise research for subsidence in Karachi. While InSAR is an effective tool for studying subsidence and ground deformation, it is limited to comparative motion between the study area and the reference point. A network of GPS stations can help enhance InSAR studies by providing ground truth data that can be used to normalize the InSAR results and obtain absolute displacement values. This ground truth data can further be enhanced by the installation of extensometers to different depths/lithologic units to gather differential subsurface compaction data. The extensometer data would help to determine which units are compacting the most, and researchers can analyze the data to find the driving causes.
The establishment of a well locations catalogue would be of great benefit to researchers as this would allow them to perform colocation and hotspot analysis. The spatial distribution of groundwater wells plays an important role in subsidence as areas with densely clustered wells have been shown to have higher subsidence rates compared to wells distributed over a larger area [56].

Conclusions
The results of this study suggest that the city of Karachi and surrounding areas are undergoing subsidence due to both anthropogenic and natural causes. Interferometric analysis of Sentinel-1 data from 2015 to 2020 shows vertical motion of 1.7 cm/yr downward, which suggests that the rate of subsidence is steadily increasing. For the first time, we report the presence of an active normal fault that strikes south-east to north-west. InSAR data demonstrate clear displacement along the fault, most notably along a 10 km-long section. This poses an alarming risk to the densely inhabited city of Karachi as there is a lack of studies on the strain accumulation of this fault, leaving the city unprepared to forecast and plan for the potential of a large seismic event.