Combined analysis of PS-InSAR and hypsometry integral (HI) for comparing seismic vulnerability and assessment of various regions of Pakistan

InSAR-based deformation analysis and the geomorphic hypsometric integral (HI) technique are powerful tools for assessing the susceptibility and comparison of seismic sites to earthquakes. Therefore, this paper mainly focuses on surface deformation analysis associated with the Mw 5.0 earthquake (2019) in Mach and Quetta, Balochistan, Pakistan. Sentinel-1 IW data was used to perform PS-InSAR time series analysis. SRTM DEM of 30 m spatial resolution was utilized for the geomorphic Hypsometry Integral (HI) method. The obtained results of the Interferogram indicate the changes in velocity and vertical displacement during pre-seismic, co-seismic, and post-seismic activity. Integral values were calculated using Hypsometry curves delineating the future probability and comparison of vulnerable seismological sites in Mach, Quetta, Ghazaband, Chamman and surroundings of Balochistan region. The combined results of HI and PS-InSAR revealed that Mach and Quetta regions are in between two lines known as the mature stages. Class 1_moderate (0.35 ≤ HI ≤ 0.52); with an integral value of HIMach = 0.398 and HIQuetta = 0.435 with a modest seismic forthcoming rate in future and susceptible to both erosion/uplifting with a vertical displacement rate more than existing ± 55 mm/year. Class 2_high (HI ˃ 0.53) with the younger and more tectonically active region surrounded by Chaman fault, which possesses a future susceptible tendency towards subsidence more than an existing velocity rate ~ 8 mm/year and Ghazaband fault towards uplifting more than 5–6 mm/year. No region of the study area was found at Monadnock: class 3_Low (HI ˂ 0.35) stabilized condition, all sites are unstable and tectonically active. Therefore, obtained results through combined PS-InSAR and HI techniques can be used for the identification of most vulnerable seismic sites and can ascertain future safe metropolitan planning.


Study area
The region of Balochistan has experienced many earthquakes in the past years with both small and high magnitudes. The most severe earthquake's magnitude was observed near Sharigh with an Mw 6.8 on 24 August 1931. Also, it was followed by the major city of Balochistan Mach, Mw 7.3 (27 August 1931), and these earthquakes triggered severe destruction and loss of lives since the significance in economic and cultural life is high in that region and relatively due to the great seismic activity of these regions. The Mach and Quetta city were selected as the target region for the earthquake displacement analysis and evaluated the seismic landforms characteristics. As the deformation related to uplifting and erosion continues in this region, potentially moderateto-large magnitude earthquakes are common and cause severe damage. The variety of deadly and destructive hazards associated with earthquakes impact a severe threat to renovation efforts and the economic development and growth of this region. Moreover, the Arabian Plate has caused subduction underneath the Eurasian Plate in this region with a 4 cm/year total rate. It is also linked with a sedimentary accretionary chunk of deposits.
The study area in Fig. 1 demonstrates the variable style of deformation (erosion and uplifting) which provides discernment into the kinematics of seismic landscape and describe a framework that can help to constrain deformational and displacement trend along the Chaman fault belt that extends 900 km to the north-eastern edge of the Makran accretionary in the south 22 .

Materials and methods
Datasets. PS-InSAR based time series analysis were carried out to monitor the earthquake induced deformation. We have used 34 images of Sentinel-1 IW data ranging March 2018 to May 2019 for pre, co, and post-seismic events, respectively are listed in Table 1. The IW-Single Look Complex product has an advantage as it moderates the signal-to-noise ratio (SNR) and scalloping effect, resulting in a uniform image quality 23 . As the sentinel-1 acquires at C-band products brings some significant advantages concerning other sensors with open accessibility and with a revisit time of 12 days 24 . The sentinel-1 acquires at C-band products brings some significant advantages concerning other sensors with open accessibility and a revisit time of 6 days. Sentinel-1 works under a pre-determined operational mode system to produce long-term consistency in data archives built especially for the long-term time series analysis. It is also one of those first five missions that European Space Agency (ESA) is emerging for Copernicus's initiative, an open data source providing free SAR-Sentinel-1 datasets 25 . SARPROZ and MATLAB platforms were used for the processing of SAR data for the ground deformation analysis and the InSAR time-series analysis. Shuttle Radar Topography Mission (SRTM) dataset with the resolution of 1 arcsecond (approximately 30 m) was used for the elevation profiles analysis and the hypsometric curves extraction in affected area. MATLAB and ArcGIS Package (ArcMap for 2D mapping and ArcScene for 3D mapping) were further used for the processing of DEM-based watersheds analysis and the geomorphic hypsometry integral  Table 1. Sentinel-1 wide swath datasets acquisition from 5 March 2018 to 23 May 2019 of (Pre, co, and post) seismic event along with track, temporal and baseline information.
The temporal baseline (day) is calculated in Table 1 to master image (25 September 2018), and the acquisitions permit a choice of pair images with the suitable baseline information and time interval.
Methods. The processing chain consists of Sentinel-1 and DEM data processing. Sentinel-1 interferometric pairs have the following components like the pre-processing sentinel-1 raw data; baseline estimation; co-registration; SPS (Sparse point Selection) using local maxima algorithm; Interferogram (single-multi) processing; Multi-Image InSAR processing; APS (atmospheric Phase screen processing); Sparse Point processing, coherence based cumulative displacement analysis and Time-series based deformation analysis. In a second chain, SRTM (1 arc-second, 30 m resolution) is used for the elevation classification based on elevation profiles and calculated hypsometric curves for each affected area. The elevation classification and area accumulation are used for the linear regression and the hypsometry integral E ≈ HI values extractions. The last landform type was estimated through HI values and multispectral 3D seismic topographic site indicated using the multispectral Landsat 8 data. An overview of the carried-out processing is shown in Fig. 2.
As the perpendicular baseline of the acquired images exhibits the coherence difference between the images. The smaller the existing perpendicular baseline between the master and slave images maintain the better overlap and demonstrate the effects of coherence among the images. The coherence has indirect proportionality using the perpendicular baseline, as with the decrease of perpendicular baseline, the coherence increases 26 . It is also recommended that to use pairs of acquired images with less than 600 m baseline values for attaining well coherence and overlaps between the salve and master datasets 27 . Hence, this study used the pair of sentinel-1 SAR images using normal minimum baselines ranging between 82.1 and − 90 m while ensuring better overlap. The pairs of sentinel-1 images below and above this range were excluded from the resultant SAR data processing. The graph exhibits the temporal/normal baseline information of datasets from 5 March 2018 to 23 May 2019 of (Pre, co, and post) seismic event of the study area as shown in Fig. 3.
After the baseline calculation, the Co-registration process was performed. Co-registration in sentinel-1 data processing was carried out by identifying the consistent pixels in both master and slave datasets. As in interferometric processing coherence is used for estimation 28 . Moreover, to attain every single imagery for the time series analysis, the co-registration process ensures that the same pixel resolution and grid gone through this process 29 . Afterward, the resultant interferograms showed the difference in phase between both co-registered master and slave images. The difference in variance in the form of colour fringes pattern on the interferogram are visible. Moreover, the idea about the generation of the interferogram is adopted from 30 Eq. (1): www.nature.com/scientificreports/ In Eq. (1), int (i.k) indicating the interferogram taken between two images (i.k) , img (j) a notation indicating the complexity of master and img * (k) is for the complex conjugate of the slave imagery. With the increases in terrain value, the phase value stays wrapped between the range from 0 and 2π radians, and cyclic get repeated. This phase information of the interferogram has already been unwrapped. It then regulates the absolute phase value taken from the relative phase value of the interferogram because of adding and multiplying integer of 2π to the edges color fringe pattern. The information about the phase unwrapping phase has given by R. Burgmann 31 in Eq. (2): where, ∅ unw is the unwrapped phase, ∅ wrapped is the phase wrapped, and the value of n is taking any positive integer. Whereas the process of interferogram phase unwrapping is used to filter and remove the effect of noise contributions in a resolution cell that contributed from the different sources. The Goldstein filtering technique is used for this processing. Goldstein's filtering method is used to remove the noise present in the resultant products.
PSI (permanent scatter interferometry) approach. Previously, Interferometric Synthetic Aperture Radar-PSI technique has been utilized to show deformation movement in the vulnerable zones using models, which in results showed the high Line of Sight (LOS) deformation velocity and land subsidence prediction at high susceptible regions 32,33 . In our study, for the selection of PS points, the amplitude stability index thresholding method is used. A threshold of 0.75 was chosen for permanent scatters as these scatters remain coherent throughout the period. While these permanent scatters (PS) have constant amplitude values in SAR datasets acquired from March 2018 to May 2019. The simplified Permanent Scatter approach fully utilized the coherence values from the temporal baseline of the images mainly in mountainous areas and gives the best results in high coherence areas with the rapid decrease in coherence value with time 34 . The temporal coherence of Permanent scatter (PS) is shown in Fig. 4.
The main procedure consists of the following steps: (2) ∅ unw = ∅ wrapped + 2nπ , www.nature.com/scientificreports/ 2D sparse LS phase unwrapping using interferograms of the co-seismic event with and without Goldstein 15 × 15 window size. ScanSAR interferometric-based analysis for the identification of (pre, co, and post) seismic events with maximum vertical displacement values. The processing chain of PSI consists of multiple segments; Removal and estimation of the atmospheric component phase can be done using a set of Spatio-temporal filtering techniques [35][36][37][38][39] , Renovating the phases into displacements values and generating the time series-based deformation maps using the atmosphere attenuation free interferograms through the APS removal method. Afterward, the identification of the PSC (Permanent scatter candidates) and removing the corresponding APS (atmospheric phase screen) attenuation, the estimation of correlation between the PSC, and the map  Hypsometric integral analysis. The hypsometric curve refers to the division of heights over an area of land.
To estimate the progressive status of landforms vulnerability which is strongly associated with the extent of uplifting and erosion that has taken place in a given region 36 . The method of Hypsometric curves interconnected and correlated to tectonic evolution and geomorphic status of the basins in terms of their progressions and placements 37 . The most valuable feature of the hypsometric curve is the geomorphic status of the basins, and its diverse dimensions can be related to each other with resembled relative area and height ratios which are the impulsive functions of the total height and total area. However, the hypsometric integral curve values are fully independent of inconsistencies invariance, dimensions, and deviations 37 . The N. Strahler has classified the types of landforms using hypsometric integral values with the typical basin dissection stages. For this purpose, the hypsometric percentage method has been used to plot the hypsometric curves. It determines the relationship involving various variables plotted against one another. In curves, the y-axis shows the relative to total elevation (h/H) values, and the x-axis represents the relative to total area (a/A) information, as shown in Fig. 5. Additionally, the relative elevation is defined as the height of a given contour (h) ratio from the even elevation to the basin's maximum elevation (H). The relative area values are taken as a proportion value of the upper contour zone (a) to the total surface outlet area of the basin (A). The resultant correlation has a range from zero to one where the lowest value of zero is at the lowest point of the drainage basin (h/H = 0), and the highest value of 1 is at the uppermost point of the drainage basin (h/H = 1) 38 . The ordinate represents the relative to total elevation (h/H) data ratio. The abscissa represents the total to relative area (a/A) ratio as shown in Fig. 5 where concave up and convex down curves exhibit youngest, young, old, and oldest basins lithology type.

Approximation of HI (hypsometric integral).
This approximation of integration is adopted from the hypsometric curve that gives the hypsometric integral (HI) values equal to the ratio of elevation and relief (E) as suggested by Pike and Richard 39 .
Mathematically, it is defined as in the Eq. (3): In Eq. (3), E refers to elevation-relief ratio, HI = hypsometric integral value, E mean defines the mean elevation; E min is for minimum elevation, and E max taking as maximum elevation. The HI value is orderly precise by the basin geometry and the area and relief ratio of the basin 40,41 and inversely associated with the total value of the steepness relief of the basin, plus channel gradients. The hypsometric integral values quantify the phases of geologic lithology expansion and erosional status. The highest range of integral corresponds to the stage of youthfulness and less eroded landform areas. The lowest integral values demonstrate the landscape is towards its maturity and old stage of the basin. The integral value is taken as an indicator of percentage in the residue of the current paralleled volume and the original volume of the area basin and a complete cycle of erosion of the basin. This entire phase or cycle of basin erosion is assembled into three classes, each demonstrating the three distinct types of the basin stages to its geomorphic cycle. The first one is known as the Monadnock stage with the integral value (HI ≤ 0.35), where the basin is completely alleviated and stabilized. The second stage is the mature stage

Combined analysis of PS-InSAR and HI.
The results obtained from the PS-SAR method were compared to the yield of measured geomorphic hypsometry integral curves values on the deflected streams networks. The combined SAR (remote sensing) and HI (geographical information system) results have provided evidence of the presence of active faults. Moreover, the integral values in the study areas show differential erosional rates, and InSAR-derived results show the deformation trend in these regions with their annual velocity rates. The hypsometry technique and PS-InSAR results will help to provide information about seismically active locations and deformation status in the future.

PS-InSAR-derived results. PS-InSAR technique has been used for displacement analysis and given
results were examined in a selected pair of sentinel-1 IW data taken for (pre, co, and post) seismic events from March 2018 to May 2019, as shown in Fig. 6. The study area from top to bottom along the vertical direction consists of sand and hard rock landform. The estimated results show that there is a constantly increasing trend in uplifting, which is observed just before the earthquake areas, as shown in red rectangles of Fig. 6. During the co-seismic event, the rate of uplifting gradually decreased and shifted temporarily to subsidence represented in red regions in red rectangular areas. Afterward, in a post-earthquake situation, areas seemed to recuperate the pre-earthquake uplifting trend in both Mach-A and Quetta-B city because of the strong reactional forces.
It is clear to note that in Fig. 6, there was an increase in the uplifting trend found before the earthquake (positive values), as shown in Fig. 6a. During the earthquake, it is clearly shown in co-seismic differential interferograms Fig. 6b-d that the rate of uplifting decreased, and in the meantime, the trend shifted towards subsidence (negative values) as shown in. Generally, there is a clear shift in the vertical displacement values between 12.0 and + 2.0 mm in along with the zero-reference line, and the topography in these areas mainly consists of thick stones. In the post-seismic interferograms Fig. 6e,f, it clearly shows that the displacement area has retained its original position, such as in pre-seismic events, because of strong reactional forces.
Time series analysis of surface deformation. Time series-based surface deformation trends measured over an area of Mach and Quetta city that has undergone earth activities are shown in Fig. 7. In this situation, the PSI simplified approach permits obtaining off the highly measured density points over the area of interest. The earthquake occurred on 16 March 2019 in the Mach and Quetta region; for this purpose, 34 pair images   www.nature.com/scientificreports/ Figure 7 illustrates the linear deformation time series analysis over areas associated with earthquakes in Balochistan using full PSI approach. The LOS (line-of-sight) linear velocity information calculated from the PSI time series datasets for the indicated area that is characterized by earthquake activities using descending track orbits. As the deformation is demonstrated in the LOS direction, the positive velocity (color bar) represents the ground movement toward the satellite. The non-positive (negative) color bar (yellow to red colors) demonstrates the motion away from the satellite. Figure  Hypsometry integral values calculation for seismic susceptibility. There are approximately nine hypsometric curves extracted at the surrounded areas of Mw 5.0 earthquakes in Mach and Quetta city to evaluate the seismic landscape characteristics and status of landforms that are strongly associated with the extent of erosion and uplifting. The gentle and steepening of slopes throughout the basins escorted by a more rapid rate change in elevation to a change of horizontal cross-sectional area. The hypsometry curves of Mach, Sibi areas exhibit steep part to the elevation in the beginning and then relatively gentle slopes belts with the integral values of HI Mach = 0.398 and HI Sibi = 0.4. This shows that the basin is at equilibrium stage where the rate of differential erosional stages is equal to the neo-tectonic level which are expressed in Fig. 8.    Topographic elevation profiles EP. Elevation Profile was extracted using the Digital Elevation Model (DEM) and Landsat 8 data to identify the 3D topography visualization of landscape characteristics. At the same time, the elevation profile between Chaman and Ghazaband fault shows the steepening of slopes from 2100 to 2258 m and the mid-section accompanied by the total elevation from 1420 to 1650 m by a more rapid rate in the decrease of elevation as shown in Fig. 9a. The elevation profile between Mach (810 m) and Sibi (129 m), at first thought, suppose that at steeper parts, the curve is towards Mach and ultimately coincides with the belts of relatively least gentle slopes curve of Sibi area with 130 m elevations as shown in Fig. 9b. Whereas in Fig. 9c, the elevation profile between Quetta (1665 m) to Mastung (1660 m) shows two big peaks in the middle of profile at 2101 m and second at 2020 m. Moreover, in Fig. 9d of Ziarat (2668 m) to Harnai (911 m), peak-valley trend could be seen from Ziarat-2668 to 2125 m. Afterwards, big dip comes in centre and slope becomes less steep from 1250 to 911 m until Harnai. In contrast, reports of major earthquakes in the surrounding cities, e.g., Quetta and Mach, showed widespread damage to cities and ground disruption. Therefore, this study implies the earthquake assessment and vulnerability analysis using geomorphological examination, Sentinel-1 based PS InSAR technique, HI, and hypsometric curves (HCs) were used to interpret the recent, youthful, and active surface deformation areas and their influence on the local topography of Balochistan.

Combined PS-InSAR and hypsometry integral (HI). The obtained results from the combined InSAR
and Hypsometry Integral approach suggested that the area of Quetta (B) and Mach (A) in Fig. 7. It experienced the subsidence with a total velocity rate − 4 to − 6 mm/year with the mature Hypsometric Integral values where the rate of uplift is nearly equal to the rate of erosion as shown Fig. 8. This region will experience a modest seismic forthcoming rate in the future and be susceptible to both erosion/uplifting with a velocity rate of more than − 6.0 mm/year. Previously this region experienced a velocity rate defined as − 5.8 mm/year. This velocity Using GIS techniques, the obtained InSAR deformation results were combined and compared to the measured geomorphic hypsometry integral curves values. Hypsometry Integral values of Quetta, Mach, Chaman fault refer towards the youthful, active tectonic regional area with an Integral values HI Mach = 0.398 and HI Quetta = 0.435 (0.35 ≤ HI ≤ 0.52); and HI ˃ 0.53 more youthful susceptible class region with S and convex-up shape. These results were previously noted in the Neogene to Quaternary sedimentary Valley, which provides a good technique for evaluating fault activity and related erosion processes 45 . These combined PS-InSAR and GIS results are evidence of the presence of active tectonics in this region and responsible for bringing out the surface deformational activities in these regions. The approach of active tectonic presence in the region using geomorphic HI values was previously distinctly identified by other researchers 46 . To compare Fig. 10A shows the time series surface www.nature.com/scientificreports/ deformation map. Figure 10B,C shows schematic cross-section of 3D topographic characteristics revised from a previous study 47 and elevation profile to describe the Lithology in the neighbouring areas. Figure 10a shows  Fig. 10f which clearly shows that region is as active as Ghazaband.
In this study, we have also proven that using the Hypsometric technique with the integral value for evaluating future seismic fault activity with the value of 0.377 towards subsidence. These hypsometry curves of Mach, Ziarat and Sibi areas are seen towards steep part to the elevation in the beginning and then relatively towards gentle slopes as shown in Fig. 8. In addition, the integral values in the corresponding regions of HI Mach = 0.398 and HI Sibi = 0.4 shows the rate of differential erosional rate is equal to the neotectonic level, and the deformation trend in these regions can exceed from the current rate up to ± 2.1 to ± 3.5 mm annual velocity rate in future. This hypsometry technique 2 , along with InSAR data, provided future seismic activity sites with deformation status and exposed that the integral values calculated by elevation-relief ratio method was accurate, less weight, and easy to estimate using GIS environment.
Evaluation of results. The selection of study areas in Balochistan, Pakistan, was done from the reports of numerous international organizations that have been in process over a legitimately long time as shown in Fig. 11. Moreover, the GPS rates of the surface deformation were also retrieved in the Quetta Valley by the University of Balochistan, the Balochistan University IT Department (QTIT), the Balochistan University Geology Department (QTAG) with the data acquired using the four GPS monitoring sites in the valley. All the stations in the valley www.nature.com/scientificreports/ were not continuously operating and were recorded at an interval of two months or greater. Therefore, the ground-based GPS rates showed the deformation at around − 7.7 mm/year to identify the subsidence in the valley, which is closer to the displacement rate of − 6.0 mm/year, calculated in our study. However, the previous studies could identify the deformation up to − 5.1 mm/year rate. For better and more accurate results, we require more working stations data. The earthquakes data used in this research paper was retrieved from the United States Geological Survey (USGS PDE), and Pakistan Meteorological Department (PMD) earthquake data catalogues as shown in Fig. 11. For identifying the total assessment in the affected earthquake region, persistent scatter interferometry (PSI) based approach has been used and inspected the total surface displacement before and after the events with help of Sentinel-1 time-series data.
The hypsometry curve technique is used to evaluate the area of seismic susceptibility sites. It is a handy tool for the identification of the most vulnerable sites. In the end, this combined approach allows us to identify the seismic areas related to young and old active tectonics sites based on the deflected streams data through google earth imagery.
The Chaman Fault (or fault group) is shown by the red color line boundary in Fig. 11, which stretches towards sub-N-S until it meets another fault named Herat in the north. Chaman Fault is a left-lateral strike-slip that has not unconfined a high earthquake activity that took place around this region in the recent past. The researcher speculates whether a substantial amount of accumulated stress has prepared this region for many more major earthquakes determined by the Ambraseys and Bilham 48 . Moreover, the Ghazaband fault area is also shown by The evidence from the assessment and analysis in the recent research led to several large earthquakes that have rattled the area in the past and contributory epochs, as shown in Fig. 12. The major cause of these devastating earthquakes lies in the tectonic and fault systems (Chaman fault; Ghazaband fault) surrounding the region. The earthquake assessment and vulnerability sites were carried out in area sources, and in terms of fault system sources 49 . The obtained results through hypsometry curves integral values can be used for the identification of seismic sites and can ascertain the engineer's future urban safe planning sites. The verification of PS InSAR analysis, all the results were firstly compared with the ones obtained from the previously conducted studies to estimate the surface deformation in the Quetta and Mach region as shown in Fig. 13.
Comparison of surface displacement results in the Quetta region: previous and recent studies. The previously conducted studies to estimate the surface deformation in the Quetta region have shown the vertical displacement   Fig. 13A to identify the impact of tectonic processes and the pattern of earthquakes in this region 42 . The previous researcher also used the small baseline technique by combining differential interferograms to calculate surface deformation velocities by phase profiles after reducing atmospheric and orbital errors 27 . On the same descending path, the deformation trend with high coherent PS points spanning from 2018 to 2019 can be seen with a total vertical displacement velocity value − 6.0 mm/year subsidence on the right side of Fig. 13B. Positive values can be considered as a shift eastward or uplift. This movement is ascribed mainly to uplift for the following reasons: (1) Tectonics of the region is endangered to thrust and north-south strike-slip faults and overriding valleys. We could not resolve the north-south strike-slip Chaman Fault movement because the strike-slip movement is perpendicular to the Line of Sight (LOS) direction. (2) The GPS observations obtained from the stations are deployed in previous studies in Quetta Valley during the year of 2006-2016 that exhibit that observed deformation is mainly towards the direction of subsidence 42 , which is also proven in the recent InSAR data processing (2018-2019) during the deformation analysis of Mach-Quetta earthquakes. This relative leading uplifting could be due to pressure increase in the Ghazaband fault along with the Quetta city and responsible for earthquake strike there as shown in Fig. 14. However, locations of uplifting and subsidence during the pre-seismic period were not stable in the study areas. After assembling the earthquake's catalogue and time-series deformation analysis were performed, it has estimated the primary uplift/subsistence of earthquakes Mw 5.0 in Mach and Quetta. Mach and Quetta earthquakes have epicentres along with the Chaman and Ghazaband fault, respectively. Temporal variations in vertical displacement in Mach and Quetta cities may cause by one of the significant surrounded faults, Chaman and Ghazaband. Figure 15 shows the graphical representation of PS-InSAR time-series displacement analysis of Mach, Quetta, Ghazaband and Chaman in Fig. 15a-d respectively with its velocity, temporal coherence and max displacement values. In their comparisons, only Ghazaband fault gives the positive values for velocity and max displacement which explains the instability and more active seismic site with respect to others.

Discussion
This section analysed and summarized the comparison results obtained through combined Persistent Scatter Interferometry (PSI) technique and Hypsometry Integral (HI) analysis in various regions of Pakistan. Most studies focused on analyzing earthquake-induced deformation using sentinel-1 based RADAR interferometry approaches 50-52 that is why we used Mach and Quetta earthquakes for our study. This study is an attempt to www.nature.com/scientificreports/ evaluate and compare seismic vulnerability through combined PSI and HI techniques. For this purpose, openly accessed data was used, and the multi-temporal PS-InSAR method was applied for the displacement analysis 51 . The high temporal resolution enabled us to acquire radar images before and after the earthquake event. The data was processed and monitored with corresponding hypsometry integral values to identify the future tendency of seismic sites 53,54 . The processing of Sentinel-1(open access) satellite data is done using the open-source software SARPROZ provided by ESA and processing tool by Periz using MATLAB programming. Moreover, the Geographic Information Systems (GIS) techniques employed to process digital elevation model (DEM) data to extract the valuable mapping representation of all the different parameters and collected data, enabling the spatial representation of the various geomorphic findings 5 . The result obtained from time-series deformation analysis indicates that the Mach-A city has observed the deformation as uplift of 7.56 mm/year and subsidence of − 53.84 mm/year. Quetta-B city experienced the deformation from − 150 mm/year to − 27.96 mm/year, as shown in Fig. 16. Additionally, the simplified approach is used in these affected areas based on the most protuberant coherent PSC's points. A similar situation was observed during the Pasni earthquake, and the total surface deformation was estimated using RADAR interferometry.
Firstly, the results obtained from PS-InSAR and then combined with the Hypsometric Integral values. The areas of Quetta and Mach experienced the subsidence with a total velocity rate from − 4 to − 6 mm/year with the mature Hypsometric Integral values as in Fig. 17. The rate of uplift is nearly equal to the rate of erosion, as discussed in the results section. This region can experience a modest seismic forthcoming rate in the future and is susceptible to both erosion/uplifting with a velocity of more than − 6.0 mm/year. Chen also introduced this hypsometric technique of tectonic regimes in terms of stream-gradient. The hypsometric analysis was used to illustrate the relative activities sites in different tectonic regimes of Western Foothills 55 . Previous studies indicate the − 5.8 mm/year velocity in this region resulted from the tectonic processes and subsequently eroded towards subsidence 42 .
The   The rate of change in regional surface deformation (along the y-axis) was estimated by the coherent scattered PS points (along the x-axis). www.nature.com/scientificreports/ that has not been eroded. In Fig. 18a, hypsometric curve of Quetta consists of a steeper slope in the basin with the Integral values of HI Mach = 0.43578. A deeper trench with the rapid fall of elevation observed at the end of the basin indicates an active tectonic landscape. The hypsometry curve of Mach was supposed to be steep at the beginning with the Integral values of HI Mach = 0.398, then gentle in the middle of the basin, and a deeper trench with rapid fall in elevation was observed in the last. It shows that the basin is at an equilibrium stage where the rate of differential erosion is equal to the neo-tectonic level. Figure 18b shows the stage of equilibrium having with stable land characteristics as defined by Strahler 43 and Fig. 18c shows Geomorphic cycle by Perez-Pena 56 which describe the geomorphic properties of the site and their comparison. The changes in hypsometric describe prospectively.
Hypsometry Integral values of Quetta, Mach, Chaman fault have been started towards the active tectonic regions with Integral values HI Mach = 0.398 and HI Quetta = 0.435 (0.35 ≤ HI ≤ 0.52); and HI ˃ 0.53 was the more susceptible class region with S and convex shape. These results were previously noted in the Neogene to Quaternary sedimentary Valley, which has proven a good technique for evaluating fault activity and related erosion processes 45 . The combination of PS-InSAR (remote sensing) and GIS (geographical Information system) results obtained through this study are evidence for active tectonics in this region. It is more likely to say that tectonic movement is responsible for causing surface deformation in the study areas. The Higher values of Hypsometry Integrals (HI) explain that it has not been eroded significantly. It is possibly formed by active tectonics and may recommend younger landscape characteristics, especially in the areas of Quetta, Chaman, and Ghazaband Faults. According to the given results of HI, we have assumed that the lower portion of the HI is perhaps related to tectonic uplifting along a Ghazaband fault, as shown in Fig. 19. Given results of the PS-InSAR and Hypsometric analysis indicates the seismic landscape vulnerability. Figs. 20 and 21. These are the indication of the active element of the Chaman and Ghazaband fault system and these faults systems surround the populated cities of Quetta and Mach, Balochistan that provides an example of the interface between co-evolving and growing structures landscapes. Drainage patterns and landform/landscape examined in the Google earth system have been used to examine the lateral propagation of thrust in different regions. These systems along the fault lines are developing a semi-dendritic pattern with less defined and migrating or deflected stream patterns Fig. 20a,c. In the mapped area, the stream patterns are interchanging irrespective of the direction of strike-strip Chaman Fault.

Evidence of faults through fluvial incision deflection. The evidence of Faults through fluvial incision deflection is shown in
The dominant control on this drainage pattern is the actively evolving structures associated with the other factors together with, bedrock lithology, land cover/land use and climate are constant in the study area and cannot explain the change in stream pattern. These streams commonly flow along with the lateral position of the fault. In these steeper terrain areas of Balochistan erosion rates are already determined using the Hypsometry Integral values (HI) and the fluvial incision has produced a variable degree of deflection. Deflected Fluvial networks along mountainsides with dendritic patterns along the active faults are also widespread in Fig. 21a,c of the Ghazaband Fault. In Fig. 21b, Time series surface deformation [mm/year] of Ghazaband with blue portion denotes movement towards radar LOS.
Their development and motion appeared to be organized by active tectonics, and therefore they hold the information about active landscape evolution. This fault line landscape evolution was also identified by the   www.nature.com/scientificreports/ uplifts and this dominant uplifting rate could be due to pioneer pressure increase in the Ghazaband fault and responsible for bringing out the Earthquake in the cities of Mach and Quetta.

Conclusion
An earthquake of Mw 5.0 stokes on March 2019 in the Balochistan region. This earthquake caused significant surface deformation in Mach and Quetta, Balochistan. This study mainly focuses on the relative susceptibility and assessment of land deformation caused by Mw 5.0 earthquakes in Mach and Quetta, Balochistan, Pakistan. Moreover, combined approach of PS-InSAR and Hypsometric techniques has been used to indicate the future tendency and comparison of seismic sites. The objective of this research has been achieved by highlighting the surface displacement in tectonically active areas that are robustly influenced by deadly earthquake activities and comparison of surrounding seismic sites were also done. The estimated deformation compared the seismic vulnerability sites using the combined techniques of RS and GIS. The multi-temporal PS-InSAR method is applied for displacement analysis before and after the seismic activity in the study areas. The Mach region shows vertical displacement of 2 mm/day and − 9.5 mm/day during the co-seismic event. The PS time series analysis shows the displacement along the Line of Sight (LOS) in the pre-seismic event of Mach from − 10 to − 37 mm/year. Moreover, the landmass of Quetta city experienced two major earthquakes of Mw 4.2 in December 2018 and Mw 5.0 and in March 2019, respectively. Total vertical displacement of − 10 cm/day and − 12 cm/day was observed due to these earthquakes. Our research strongly observes that Mach and Quetta regions are between two lines known as the mature stage: class 1_moderate (0.35 ≤ HI ≤ 0.52); with integral HI Mach = 0.398 and HI Quetta = 0.435 with a modest seismic approaching rate in the future. Whereas, surrounding areas like Sibi, Harnai, Ziarat and Mastang are comparable less active than Mach, Quetta, Chaman and Ghazaband. These areas are also vulnerable to erosion/uplifting because of vertical displacement rates of more than ± 55 mm/year. However, class 2_high (HI ˃ 0.53) with the younger and more tectonically active region is surrounded by Chaman fault, which possesses future vulnerability regarding subsidence more than the existing velocity rate ~ 8 mm/year and Ghazaband fault showed unequalled uplifting, which is more than 5-6 mm/year rate. Chaman, Ghazaband, Quetta and Mach are found to be most seismic vulnerable and active regions where the probability of earthquakes are greater than other surrounding regions. Our study limitations were all the stations in the valley were not continuously operating and were recorded at an interval of two months or greater. Therefore, the ground-based GPS rates showed the deformation at around − 7.7 mm/year to identify the subsidence in the valley which is closer to the displacement rate of − 6.0 mm/year calculated in our study. However, the previous studies could identify www.nature.com/scientificreports/ the deformation up to − 5.1 mm/year rate. We can overcome this by placing more ground-based continuously working GPS stations to get the results more accurate. In the future, we will endeavour to investigate the pressure between the plates and assess the seismically vulnerable sites during surface rupture.