Analysis of the Deformation Behavior and Sinkhole Risk in Kerdabad, Iran Using the PS-InSAR Method

: The unexpected collapse of land surface due to subsidence is one of the most signiﬁcant geohazards that threatens human life and infrastructure. Kabudrahang and Famenin are two Iranian plains experiencing several sinkholes due to the characteristics of the underground soil layers and extreme groundwater depletion. In this study, space-based Synthetic Aperture Radar images are used to investigate the ground displacement behavior to examine the feasibility of Sentinel-1 data in detecting precursory deformation proceeding before the sinkhole formation. The selected sinkhole occurred in August 2018 in the vicinity of Kerdabad village in Hamedan province with a 40 m diameter and depth of ~40 m. Time series of the European constellation Sentinel-1 data, spanning from January 2015 to August 2018, is analyzed, and the results revealed a 3 cm annual subsidence (–3 cm/year) along with the line-of-sight direction. Time-series analysis demonstrated that the driving mechanism of the sinkhole formation had a gradual process. Displacement of persistent scatterers (PSs) near the cave area had an acceleration by approaching the sinkhole formation date. In contrast, other areas that are far from the cave area show linear subsidence behavior over time. Additionally, the one-kilometer deformation proﬁle over the cave area indicates a high subsidence rate precisely at the location where the sinkhole was formed later on 20 August 2018.


Introduction
Collapsing of the surface layer that is mainly caused by karst processes [1], chemical dissolution of carbonate rocks [2], or suffusion processes [3] are the primary catalysts of sinkhole formation. From a geohazard point of view, a sinkhole, according to its potential risks and unpredictability, can pose significant hazards for people, properties, and agricultural activities [4]. Different factors, such as karst processes and overexploitation of groundwater, can make the ground susceptible to a sinkhole.
Dissolution of soluble minerals such as halite, anhydrite, gypsum, calcite, and dolomite has been identified as the major reason for karst formation. Gypsum is several times more soluble than carbonate rocks [5]. Therefore, karst in gypsum rocks evolves at a faster pace and often causes severe problems to the built environment, frequently fostered by anthropogenic effects [6]. Overexploitation of groundwater creates small inner cavities, and penetration of surface water into the subsurface layers, which is a result of high solubility, is another main reason for sinkhole formation [7].
According to different perspectives regarding geological properties and the formation rate, sinkholes are classified into various categories [8]. In the areas with limestone outcrops, We assume that precursory deformation may occur before sinkhole formation, so we have exploited the SAR data archive before the sinkhole formation date. Several factors have made the study of displacement behavior in small areas challenging, such as sinkholes with radar images and InSAR techniques. These problems include low spatial resolution (for Sentinel-1~20 × 5 m) and extensive vegetation coverage in the sinkhole area, which leads to the limitation of PS selection in the low coherence area. Additionally, due to the lack of prior knowledge and unpredictable location of the cave, we need to examine the spatial deformation behavior in a sizeable area. For this purpose, the Sentinel-1 time-series data and the spatiotemporal deformation behavior have been investigated.
The study area and its geological characteristics are described in Section 2. Materials and Methodology are discussed in Section 3. Finally, Sections 4 and 5 are dedicated to the Results and Discussion, and Conclusions, respectively.

Study Area and Geological Setting
The study area is located in the northeast of Hamedan and approximately encompasses the area between longitudes 48 • 44 to 48 • 49 E and latitudes of 35 • 5 to 35 • 9 N WGS84 ( Figure 1). The area is confined to Kabudrahang, Razan, and Ghahavand plains from northern, northeastern, and southeastern parts, respectively. This region is a smooth plain having an overall north-to-south slope that is intersected by the Hamedan-Tehran Road, on the side where the Mofatteh power plant is located. The Kabudrahang and Famenin are two adjacent plains, covering 1348 and 800 square kilometers, respectively, considered the lowest plains topographically in the semi-arid Hamadan state in west Iran, facing groundwater depletion and sinkhole formation [50].

Materials and Methods
We acquired 76 Sentinel-1 images, all in ascending direction. Processing of all the images (76 SLC images) leads to a decrease in the amplitude dispersion index value and a low value of ADI, which causes a reduction in the PSs density or complete loss of PSs in the cave area. So, all images were divided into two groups, the first group represents the images taken from January 2015 to December 2016, and the second group includes the pe- Hamedan state consists of the central plains of Kabudrahang, Razan, Famenin, Ghahavand-Komijan, Hamedan Bahar, and Markazi plains. The Kabudrahang and Famenin plains are located on the northern border of the Sanandaj-Sirjan zone and the northern part of the central Zagros thrust belt [51]. On the surface of outcrops or heights adjacent to Kabudrahang and Famenin plains, faults are described as follows (see Figure 1a):

1.
Faults located in the highlands of Pashalu mountain, situated in the northeast of the area with a generally east-west direction with a length of 3 to 28 km. 2.
Faults located in the Achini mountain highlands with a general north-south direction of 5 to 20 km.

3.
Potential buried faults of the Kabudrahang-Famenin plains having a general westsoutheast direction of 25 km.
Drillings in the Famenin and Kabudrahang fields show that the bedrock is made from Oligo-Miocene limestone, and the karst structure is predominant [52]; the bedrock in the area is highly prone to dissolving. Additionally, a significant amount of dissolved free carbon dioxide has been detected in the wells, representing a high CO 2 pressure associated with the wells [52,53]. Investigating the gathered information and samples taken from Ghahavand and Famenin reveals that underground water containing CO 2 gas has impacted the region's bedrock by creating cavities, flaws, and cracks through dissolving a high proportion of available calcium carbonate [53]. The transformation process into karst happened inside the bedrock so that its porosity is about 20 to 25% [54].
The sinkhole studied in this research was formed in the last days of August 2018. It is worth mentioning that the occurrence date is remarkable, especially when there is a strong relationship between the groundwater level and land subsidence or sinkhole formation during summers [39]. In addition, overexploitation influences water levels in different seasons, and proportionally, most sinkholes can be found near wells excavated on limestone bedrock [51]. One of the major sources of water withdrawal in the area is the Mofatteh power plant situated on the Hamedan-Tehran Road, located in the southeastern part of the plain since 1994. It should be noted that the active power plant in the study area requires water for cooling and steam turbine; for this purpose, 50 wells with a depth of 80 to 120 m in the southeastern part of the power plant were drilled [55]. Deep wells with a pumping rate of 500 l/s have significantly lowered the underground water level of the southern and eastern regions of the Mofatteh power plant. Statistics indicate an average consumption of 15 million cubic meters of water for the power plant, so that in this area, the water drop in nine years (1992 to 2001) was about 50 m [52]. Sinkholes in the state of Hamedan consist of 50% subsidence, 43% collapse, and 7% solution sinkholes [56]. In total, most of the sinkholes in the Kabudrahang and Famenin plains are subsidence and cover sinkholes [56].
The sinkhole under analysis in this research is a subsidence sinkhole located 2 km southeast of Kerdabad, Hamedan, and was formed on 22 August 2018. The Kerdabad region's bedrock is made of carbonate, and the sinkhole is formed by the rupture of the soil on the lower karst limestone. The solubility of calcium carbonate in solutions containing carbon dioxide is an essential factor in the dissolution of carbonate bedrock. Moreover, since the bedrock of Kerdabad is karst, this sinkhole is formed with a depth and width of 40 m ( Figure 1 depicts the location of the study area). Characteristics of major sinkholes in the Hamedan state are presented in Table 1.

Materials and Methods
We acquired 76 Sentinel-1 images, all in ascending direction. Processing of all the images (76 SLC images) leads to a decrease in the amplitude dispersion index value and a low value of ADI, which causes a reduction in the PSs density or complete loss of PSs in the cave area. So, all images were divided into two groups, the first group represents the images taken from January 2015 to December 2016, and the second group includes the period of January 2017 to August 2018. Therefore, the reason for dividing the database into two separate groups is to process fewer dataset images to obtain PS candidates with high ADI values in the cave area. In addition, the overlap of the two groups has been considered to compare their displacements during the same months; therefore, both groups start from January. Additionally, the data of the second group is closer to the date of sinkhole formation (see Table 2). In the first instance, using a co-registered stack total of 28 and 46, interferograms were generated from first-and second-group S-1A datasets using a star graph (see Figure 2). The methodology implemented in this study is based on the standard PS-InSAR method [57]. This technique selects image pixels with excellent phase stability over the SAR data acquisitions period, called persistent or permanent scatterers (PS) [57], and measures surface deformation over months or years. These pixels have strong backscattering to the SAR sensor to generate a reliable deformation velocity map. The powerful PS-InSAR method overcomes the shortcomings of the conventional differential SAR interferometry (DInSAR) method and makes it possible to create a deformation time-series at millimeter precision. Thus, this is a suitable method for investigating ground displacement behavior over time because stable pixels (PSs) over time are features of sub-pixels that are not affected by temporal and spatial decorrelation [58].
The temporal and spatial separation of image pairs is called the temporal and spatial baseline, respectively. Using the baselines of each image, we can set the spatiotemporal space interferometric configuration to use for selecting the master image(s) and generation interferometric image pairs (interferograms) [59].
The master image is selected approximately in the middle of interferometric configuration to minimize the spatiotemporal decorrelation [60]. The chosen master image is used for co-registration of all SLC images (slave images) into the single master stack and forma-tion of the high coherence interferograms [61]. Finally, by a single master configuration for N number of SAR acquisitions, we can generate N − 1 interferograms. Each pixel in the interferograms is made of several terms, which are the orbital phase, flat earth effect, and topography phase component that can simply be removed by orbital information, simulated phase image related to imaging geometry, and external digital elevation model, respectively. Finally, the remaining terms as in the following equation [58] where the first term is the residual topographic phase due to external DEM inaccuracy (h), the second indicates linear deformation velocity between acquisitions (v), ϕ i atm is the atmospheric phase effect, and ϕ i noise denotes decorrelation errors such as baseline estimation error, temporal, and geometrical decorrelations. The other parameters are as follows: It can be observed that to estimate the deformation value, it is necessary to remove all other terms except the deformation term. For this reason, coherent pixels should be selected to overcome some shortcomings of conventional InSAR, detailed information of which can be found in [57]. As the atmosphere is correlated in the spatial domain, this term can be ignored in the adjacent pixels in the spatial graph [62].
The amplitude stability index (ASI), which corresponds to a dispersion of amplitude (DA), is used for selecting pixels that have a low phase standard deviation as PS candidates. This step can be implemented by thresholding on a reflectivity map or amplitude stability index map. A threshold of~DA < 0.25 (ASI > 0.75) is recommended by Ferretti et al. [57], but in different situations, other values may be chosen. The PSCs as nodes establish a spatial network using Delaunay triangulation [57], which can be used for estimating unknowns on the edges (between nodes) using the maximization of the periodogram function [58]. Then, using at least one pixel reference point on a non-deformed area integrates all the estimations to achieve the absolute value at each PSC. Later, based on spatially and temporally correlated components [62], the atmospheric phase can be roughly estimated. Then, increasing the number of PSCs can be obtained by a lower threshold of ASI, and the parameters are re-estimated for whole final PSs. Estimations are assessed through analyzing temporal coherence to obtain guaranteed results. Finally, the deformation values of the PS points are used to generate the velocity map.
In this study, SARPROZ software, which utilizes the standard PSI method, was used for estimating deformation in PS-InSAR [63]. The PS-InSAR approach has been selected instead of SBAS to gain acceptable numbers of PSs in the sinkhole area. The applied workflow of processing steps is shown in Figure 2.
The mentioned PSInSAR method was implemented on ascending Sentinel-1 data in the IW imaging mode. In this study, the ground deformation trend was estimated using SARPROZ software [63]. The master image was selected based on the maximum correlation between temporal and spatial baseline and doppler centroid frequency for both groups [60]. First-and second-group slave images were co-registered concerning the master images acquired on 17 February 2016 and 14 November 2017, respectively. Using a singlemaster baseline configuration (star graph), which guarantees high coherence of whole interferograms, for the first and second groups, 28 and 46 interferograms were created, respectively (see Figure 3). For eliminating the main part of a topographic phase, the 1-arcsecond shuttle radar topography mission (SRTM) digital elevation model was used. The mentioned PSInSAR method was implemented on ascending Sentinel-1 data in the IW imaging mode. In this study, the ground deformation trend was estimated using SARPROZ software [63]. The master image was selected based on the maximum correlation between temporal and spatial baseline and doppler centroid frequency for both groups [60]. First-and second-group slave images were co-registered concerning the master images acquired on 17 February 2016 and 14 November 2017, respectively. Using a single-master baseline configuration (star graph), which guarantees high coherence of whole interferograms, for the first and second groups, 28 and 46 interferograms were created, respectively (see Figure 3). For eliminating the main part of a topographic phase, the 1-arcsecond shuttle radar topography mission (SRTM) digital elevation model was used. After removing the effect of the flat Earth and topographic phase, PS candidates were selected using the ASI criteria, and then the spatial network was formed by Delaunay  The mentioned PSInSAR method was implemented on ascending Sentinel-1 data in the IW imaging mode. In this study, the ground deformation trend was estimated using SARPROZ software [63]. The master image was selected based on the maximum correlation between temporal and spatial baseline and doppler centroid frequency for both groups [60]. First-and second-group slave images were co-registered concerning the master images acquired on 17 February 2016 and 14 November 2017, respectively. Using a single-master baseline configuration (star graph), which guarantees high coherence of whole interferograms, for the first and second groups, 28 and 46 interferograms were created, respectively (see Figure 3). For eliminating the main part of a topographic phase, the 1-arcsecond shuttle radar topography mission (SRTM) digital elevation model was used. After removing the effect of the flat Earth and topographic phase, PS candidates were selected using the ASI criteria, and then the spatial network was formed by Delaunay After removing the effect of the flat Earth and topographic phase, PS candidates were selected using the ASI criteria, and then the spatial network was formed by Delaunay triangulation. In our analysis, PS points were chosen to have a DA lower than 0.35 (ASI > 0.65). This resulted in 23712 PSC points for the first group and 24204 PSC points for the second group. The atmospheric phase component was estimated in the location of the PSCs, which were selected based on the ASI.
For each connection (edge) among points, the DEM error and amount of displacement via optimized periodogram function are estimated. A single reference point is taken based on the selected persistent scatterer candidates (PSCs), which are not affected by subsidence. The calculated parameters are subtracted from the observed interferometry phase, and only atmospheric effects and noise remained. Due to the atmospheric phase's spatial correlation, applying a low pass filter will allow us to separate this term from noise. Finally, the deformation velocity of PSs was used to generate velocity deformation maps.

Results and Discussion
The linear model assumption was exerted for the first and second groups' results to estimate the ground deformation and generate a displacement map. The velocity maps from the analyses are shown in Figure 4a,b. The displacement rates in the line of sight (LOS) direction were found in the range between −30 mm/year and +10 mm/year and −27 mm/year and +10 mm/year for groups 1 and 2, respectively. The annual velocity map indicates that we have the highest deformation at the southeast of the study area. The black circle represents the cave location (see Figure 4). As can be seen in the annual velocity map of the second group, the subsidence rate increases over time near the cave area (see Figure 4b).
displacement via optimized periodogram function are estimated. A single reference point is taken based on the selected persistent scatterer candidates (PSCs), which are not affected by subsidence. The calculated parameters are subtracted from the observed interferometry phase, and only atmospheric effects and noise remained. Due to the atmospheric phase's spatial correlation, applying a low pass filter will allow us to separate this term from noise. Finally, the deformation velocity of PSs was used to generate velocity deformation maps.

Results and Discussion
The linear model assumption was exerted for the first and second groups' results to estimate the ground deformation and generate a displacement map. The velocity maps from the analyses are shown in Figure 4a,b. The displacement rates in the line of sight (LOS) direction were found in the range between −30 mm/year and +10 mm/year and −27 mm/year and +10 mm/year for groups 1 and 2, respectively. The annual velocity map indicates that we have the highest deformation at the southeast of the study area. The black circle represents the cave location (see Figure 4). As can be seen in the annual velocity map of the second group, the subsidence rate increases over time near the cave area (see Figure  4b).  In the following, several subjects, such as spatial and temporal deformation analysis, are explored. At the first step in Section 4.1, the deformation histogram velocities for the first and second groups are discussed. In the next part, spatiotemporal deformation analysis is investigated for both groups. Following this, to determine the uniqueness of the displacement behavior in the sinkhole prone area, the time-series behavior of the PSs on the cave region is compared with two other regions. Finally, the levels of piezometric wells and their relationship with the nearest PS time-series deformation are discussed in Section 4.4.

Progressive Subsidence
Velocity displacement maps for each group (first and second) reveal progressive subsidence in the study area. This can be determined by investigating peaks of histograms of PSs velocities (see Figure 5). One can conclude that the black curve, which is for the first group displacement rates (from 5 January 2015 to 25 December 2017), mostly includes displacements between +5 and −5 mm/year (histograms peak is in the range of zero deformation). In comparison, the second group's diagram (red curves), which is drawn from 01/18/2017 to 08/17/2018, has tended to the left, including higher subsidence regions. This means the subsidence zone has increased in the area before sinkhole formation.
Remote Sens. 2021, 13, 2696 9 of 18 sidence in the study area. This can be determined by investigating peaks of histograms of PSs velocities (see Figure 5). One can conclude that the black curve, which is for the first group displacement rates (from 5 January 2015 to 25 December 2017), mostly includes displacements between +5 and −5 mm/year (histograms peak is in the range of zero deformation). In comparison, the second group's diagram (red curves), which is drawn from 01/18/2017 to 08/17/2018, has tended to the left, including higher subsidence regions. This means the subsidence zone has increased in the area before sinkhole formation.

Investigating Temporal and Spatial Deformation Behavior of the Cave Area
This section investigates the temporal and spatial behavior of ground deformation in the cave area. In this regard, sinkhole prone area PSs have to be selected and plot their time-series deformation trend. This leads to a better understanding of PSs temporal deformation behavior. Accordingly, eight PSs, which include boundaries inside and outside the cave area, were selected (see Figure 6).

Investigating Temporal and Spatial Deformation Behavior of the Cave Area
This section investigates the temporal and spatial behavior of ground deformation in the cave area. In this regard, sinkhole prone area PSs have to be selected and plot their timeseries deformation trend. This leads to a better understanding of PSs temporal deformation behavior. Accordingly, eight PSs, which include boundaries inside and outside the cave area, were selected (see Figure 6). Each PS diagram is divided into black and red parts (see Figure 7). Separation points for these diagrams are located on 26 March 2018 and 12 March 2016 for the second and first groups, respectively. To select the separation point, a third-order polynomial model was fitted, and its curvature is chosen as the separation point where the slope changes. Then, for showing their trends, each part of the diagram was interpolated with the linear method. As can be observed (Figure 7b), points outside the cave (the PSs number 1, 2, 3, 7, and 8) have the same slopes at each part of the diagrams (black and red). In other words, there is no difference in the interpolated displacement trend over time in the two parts. However, for the points inside the cave area (PSs number 4, 5, and 6), the second part of the diagram receives a steeper slope, which indicates a different behavior compared to prior sinkhole formation. Each PS diagram is divided into black and red parts (see Figure 7). Separation points for these diagrams are located on 26 March 2018 and 12 March 2016 for the second and first groups, respectively. To select the separation point, a third-order polynomial model was fitted, and its curvature is chosen as the separation point where the slope changes. Then, for showing their trends, each part of the diagram was interpolated with the linear method. As can be observed (Figure 7b), points outside the cave (the PSs number 1, 2, 3, 7, and 8) have the same slopes at each part of the diagrams (black and red). In other words, there is no difference in the interpolated displacement trend over time in the two parts. However, for the points inside the cave area (PSs number 4, 5, and 6), the second part of the diagram receives a steeper slope, which indicates a different behavior compared to prior sinkhole formation. Additionally, the time-series deformation trend for the first group (January 2015 to December 2017) has been investigated. As shown (Figure 7a), there is no apparent change in the diagram slope at the black and red parts. For this part, the behavior observed in the PSs numbers 4, 5, and 6 have stable deformation behavior; consequently, the cave area is affected by deformation only before sinkhole formation (See Figure 7b, PSs numbers 4, 5, and 6).
Moreover, to study the local spatial deformation behavior over the cave area, a onekilometer profile (ij dashed line in Figure 6) is drawn for the first and second groups (see Figure 8). According to Figure 8a, there is local subsidence of 40 mm in the sinkhole area for the second group (the displacement at the sinkhole's area has increased), while for the first group (Figure 8b), the whole diagram fluctuated between +10 and -10 and there is no specific displacement behavior. In other words, the proximity of the near-sinkhole PSs shows homogeneously constant velocities and deformation trends (Figure 8b), but the Additionally, the time-series deformation trend for the first group (January 2015 to December 2017) has been investigated. As shown (Figure 7a), there is no apparent change in the diagram slope at the black and red parts. For this part, the behavior observed in the PSs numbers 4, 5, and 6 have stable deformation behavior; consequently, the cave area is affected by deformation only before sinkhole formation (See Figure 7b, PSs numbers 4, 5, and 6).
Moreover, to study the local spatial deformation behavior over the cave area, a onekilometer profile (ij dashed line in Figure 6) is drawn for the first and second groups (see Figure 8). According to Figure 8a, there is local subsidence of 40 mm in the sinkhole area for the second group (the displacement at the sinkhole's area has increased), while for the first group (Figure 8b), the whole diagram fluctuated between +10 and −10 and there is no specific displacement behavior. In other words, the proximity of the near-sinkhole PSs shows homogeneously constant velocities and deformation trends (Figure 8b), but the deformation observed in sinkhole-prone land PSs shows a unique pattern similar to a sink (Figure 8a). deformation observed in sinkhole-prone land PSs shows a unique pattern similar to a sink (Figure 8a).

Comparison of the Displacement Behavior of the Cave Area with the Highest Subsidence Region and the Area Closest to the Cave
For investigating the temporal behavior in the second group (from January 2017 to August 2018), three regions were chosen (see Figure 9). The properties of regions are described in Table 3. Table 3. Characteristics of three regions whose time-series behavior was analyzed.

Regions
Description Region 1 Nearest subsidence area to the cave Region 2 cave area Region 3 Near to Mofatteh power plant In each region, four PSs were chosen as representatives for the time-series displacement ( Figure 10). Diagram depicts an increasing slope of trend only in the cave area

Comparison of the Displacement Behavior of the Cave Area with the Highest Subsidence Region and the Area Closest to the Cave
For investigating the temporal behavior in the second group (from January 2017 to August 2018), three regions were chosen (see Figure 9). The properties of regions are described in Table 3. deformation observed in sinkhole-prone land PSs shows a unique pattern similar to a sink (Figure 8a).

Comparison of the Displacement Behavior of the Cave Area with the Highest Subsidence Region and the Area Closest to the Cave
For investigating the temporal behavior in the second group (from January 2017 to August 2018), three regions were chosen (see Figure 9). The properties of regions are described in Table 3. Table 3. Characteristics of three regions whose time-series behavior was analyzed.

Regions
Description Region 1 Nearest subsidence area to the cave Region 2 cave area Region 3 Near to Mofatteh power plant In each region, four PSs were chosen as representatives for the time-series displacement ( Figure 10). Diagram depicts an increasing slope of trend only in the cave area  Table 3. Characteristics of three regions whose time-series behavior was analyzed.

Region 1
Nearest subsidence area to the cave Region 2 cave area Region 3 Near to Mofatteh power plant In each region, four PSs were chosen as representatives for the time-series displacement ( Figure 10). Diagram depicts an increasing slope of trend only in the cave area (Figure 10b, region 2). As shown in Figure 10a,c, the trend slope at red and black colored parts are equal for regions 1 and 3, but the subsidence rate increases for the cave area at the red dots, which is close to the time of the collapse. ( Figure 10b, region 2). As shown in Figure 10a,c, the trend slope at red and black colored parts are equal for regions 1 and 3, but the subsidence rate increases for the cave area at the red dots, which is close to the time of the collapse.

Deformation Behavior and Piezometric Wells
The groundwater level is one of the most critical factors influencing subsidence. Groundwater, due to its volume and pressure, can keep the substrates stable or cause instability. Accordingly, overexploitation leads to a disturbance in the underground layers.
The drop in the region's water level due to overexploitation and subsequent compaction in the water table is the primary explanation for subsidence in this area. Excessive water withdrawal in the study area can be seen in the Khalatabad and Noabad piezometric wells (see Figure 11). As can be seen, the water level of the Khalatabad piezometric well has decreased sharply since 2013 and the Noabad piezometric well since 2007.

Deformation Behavior and Piezometric Wells
The groundwater level is one of the most critical factors influencing subsidence. Groundwater, due to its volume and pressure, can keep the substrates stable or cause instability. Accordingly, overexploitation leads to a disturbance in the underground layers.
The drop in the region's water level due to overexploitation and subsequent compaction in the water table is the primary explanation for subsidence in this area. Excessive water withdrawal in the study area can be seen in the Khalatabad and Noabad piezometric wells (see Figure 11). As can be seen, the water level of the Khalatabad piezometric well has decreased sharply since 2013 and the Noabad piezometric well since 2007.
As can be observed in Figure 9, the Khalatabad and Noabad wells are the nearest ones to the deformation velocity map coverage, with a distance of 4 and 1.7 km from the cave, respectively (blue circles depict locations of these wells in Figure 9). According to PS deformation behavior near the Khalatabad well (see Figure 12, Khalatabad a1 and b1) and the water level fluctuations in the piezometric well, it can be seen that displacement changes in different sections are consistent with the piezometric water level fluctuations that are also true for Noabad well (Figure 12, Noabad a2 and b2). It should be noted that the break-in Noabad diagram is due to a lack of data. To summarize, taking a look at different parts of the diagram, one could note a proportional behavior between subsidence and piezometric changes in wells. To investigate this relationship quantitatively, the data have been interpolated using a linear model f(x) = p 1 *x + p 2 (Figure 12), and the slope of each diagram can be expressed by the value of p1. As can be seen, there is a positive slope in the trend of the Khalatabad piezometric well, and this positive slope has revealed itself in the uplift of the time series. In contrast, in the Noabad well, the general trend is negative with the p1 value of −0.2399. Correspondingly, the time-series behavior of ground displacement shows it as the subsidence.

Deformation Behavior and Piezometric Wells
The groundwater level is one of the most critical factors influencing subsidence. Groundwater, due to its volume and pressure, can keep the substrates stable or cause instability. Accordingly, overexploitation leads to a disturbance in the underground layers.
The drop in the region's water level due to overexploitation and subsequent compaction in the water table is the primary explanation for subsidence in this area. Excessive water withdrawal in the study area can be seen in the Khalatabad and Noabad piezometric wells (see Figure 11). As can be seen, the water level of the Khalatabad piezometric well has decreased sharply since 2013 and the Noabad piezometric well since 2007.  Groundwater level decline is one of the main causes of subsidence and sinkhole formation in this area. As shown in Figure 11, over a period of 20 years, there has been a severe decrease in water levels of 41 m for Khalatabad and 85 m for Hamedan. This is one of the causes of subsidence and sinkhole occurrence in the region. The empty karst spaces created by groundwater suction provide a large amount of room for sediment suffusion (see Figure 11, Khalatabad and Noabad stratigraphic columns), but it is not the only factor. The karst bedrock is covered with a thick overburden of soil and fine-grained sediments (alluvial deposits) to a depth of about 70 to 80 m. In this soil mass, due to the suction of sediments as a consequence of groundwater over-exploitation, a series of internal caves have been created. These caves, when the water table reaches the contact of sediment and karst, cause an imbalance of the retaining forces and initiate sinkhole formation. These factors must occur together, so there are many aquifers where groundwater levels declined, but no sinkholes have formed.

Uncertainty Assessment
Due to the absence of a GPS station or leveling in the study area, it was impossible to evaluate our velocity and time series displacement against external observations. So, the spatial statistical uncertainty of the velocities was calculated to determine the error of the estimations. Assuming the PSs velocity of a selected region is zero (no displacement), the error of the velocity estimates in this non-subsided area is based on the Equations (2) and (3) [65].
where MEAN, STD, and N represent mean velocity values, the standard deviations of the velocities, and the number of PS points in the selected region, respectively. The village of Vian in the southwest of the study area has little movement, so this area was chosen as a stable area without displacement. The average velocity is 0.85 mm/year, and 1778 points were selected in this area. Figure 13 shows the velocities histogram of this area. The error value calculated in this area based on Equations (2) and (3) is equal to 0.8 mm per year. different parts of the diagram, one could note a proportional behavior between subsidence and piezometric changes in wells. To investigate this relationship quantitatively, the data have been interpolated using a linear model f(x)=p1*x+p2 (Figure 12), and the slope of each diagram can be expressed by the value of p1. As can be seen, there is a positive slope in the trend of the Khalatabad piezometric well, and this positive slope has revealed itself in the uplift of the time series. In contrast, in the Noabad well, the general trend is negative with the p1 value of -0.2399. Correspondingly, the time-series behavior of ground displacement shows it as the subsidence.  Groundwater level decline is one of the main causes of subsidence and sinkhole formation in this area. As shown in Figure 11, over a period of 20 years, there has been a severe decrease in water levels of 41 m for Khalatabad and 85 m for Hamedan. This is one of the causes of subsidence and sinkhole occurrence in the region. The empty karst spaces created by groundwater suction provide a large amount of room for sediment suffusion (see Figure 11, Khalatabad and Noabad stratigraphic columns), but it is not the only factor. The karst bedrock is covered with a thick overburden of soil and fine-grained sediments (alluvial deposits) to a depth of about 70 to 80 m. In this soil mass, due to the suction of sediments as a consequence of groundwater over-exploitation, a series of internal caves have been created. These caves, when the water table reaches the contact of sediment and karst, cause an imbalance of the retaining forces and initiate sinkhole formation. These factors must occur together, so there are many aquifers where groundwater levels declined, but no sinkholes have formed.

Uncertainty Assessment
Due to the absence of a GPS station or leveling in the study area, it was impossible to evaluate our velocity and time series displacement against external observations. So, the spatial statistical uncertainty of the velocities was calculated to determine the error of the estimations. Assuming the PSs velocity of a selected region is zero (no displacement), the error of the velocity estimates in this non-subsided area is based on the Equations (2) where MEAN, STD, and N represent mean velocity values, the standard deviations of the velocities, and the number of PS points in the selected region, respectively.
The village of Vian in the southwest of the study area has little movement, so this area was chosen as a stable area without displacement. The average velocity is 0.85 mm/year, and 1778 points were selected in this area. Figure 13 shows the velocities histogram of this area. The error value calculated in this area based on Equations (2)

Conclusion
This research attempted to estimate the deformation of land prone to sinkholes using the PS-InSAR technique and to evaluate the feasibility of Sentinel-1 data in the detection of precursory signal before the sinkhole formation. The sinkhole formed on 22 August 2018 was selected for this purpose. We used the Sentinel-1 IW imaging mode data stack between 5 January 2015 and 17 August 2018 (a few days before sinkhole formation). The data were divided into two groups: (1) images that were taken between January 2015 and December 2016, and (2) images that were taken between January 2017 and August 2018. These groups were analyzed separately. In total, approximately −3 cm annual subsidence is observed on the entire period of the data stack.
By investigating and analyzing the temporal and spatial behavior of the cave area, a noteworthy deformation behavior before sinkhole formation was detected; i.e., a localized subsidence is observable exactly at the location where the sinkhole later formed on 20 August 2018. Additionally, the subsidence velocity of the PSs located on the cave was increased with time; therefore, local deformation can be detected before sinkhole formation. Time-series deformation behavior reveals steeper subsidence in the cave area at least three months before occurrence. According to the second data group results, three regions were chosen for further investigations. The results imply that noteworthy behavior and subsidence patterns can only be detected in the cave region. The analyses performed in this research detected precursors for small phenomena such as sinkholes. Additionally, the introduced method can be generalized to phenomena such as landslides, collapses, settlements, or any other phenomenon with gradual deformation. Since the studied sinkhole is located in the vicinity of a power plant, identifying the precursor helps design an early warning system for proper and timely management. However, the authors seek to use high-resolution SAR images in different bands (X and L bands) to find a more significant number of PSs for better analysis. Although the X-band provides a higher resolution, decorrelation effects due to vegetation are still the main limitations of the sinkhole study that make the retrieval of meaningful information more difficult. It is worth mentioning that the InSAR method only can study surface deformation, so other technologies and techniques must be used in the design of an EWS.