*4.1. Result of PS-InSAR*

The ground deformation rate along the line of sight (LOS) direction during the acquisition time interval was obtained based on the PS-InSAR technique (see Section 3.1). SARPROZ (SAR PROcessing tool by periZ) [51] software was used to implement the PS-InSAR technique in the current research paper. The star configuration of the S-1A SLC time-series data stack is shown in Figure 4. The master image for the dataset was selected based on maximizing the stack coherence [68]. All slave images were co-registered with respect to the master image (2019/09/11). The shuttle radar topography mission (SRTM) DEM was used to remove topographic-related phase components from the interferometric phases. After the selection of 7881 PSCs, a spatial network was created by Delaunay Triangulation, connecting each point to the other. The unknown LOS velocity and DEM error were calculated along with the connections by maximizing a periodogram. All the obtained

parameters were then integrated into the absolute values with respect to a reference point, which had no subsidence rate. The atmospheric phase for PSCs can be resampled on the uniform image grid as the APS. With having APS compensated for all slave images, the unknowns were estimated again for more PS points, selected by a lower threshold on ADI to obtain a more dense subsidence map. Figure 5a shows the LOS deformation map. According to Figure 5a, the maximum velocity was about −175 mm/year, which occurred in the southern agricultural part of the region of interest, where more ground-water extraction was observed.

Further, the LOS deformation rates should be decomposed into horizontal and vertical deformation components as it is the inherently vertical movement of the Earth's surface with a slight horizontal displacement. It has been proved that the horizontal deformation is a very small portion of motion compared to the vertical deformation [67,69]. Hence, the LOS deformation could be assumed as negligible and converted simply into the vertical deformation rates using the cosine of the incidence angle of the radar signal. The interpolated land subsidence inventory map was designed based on the vertical velocity deformation map for Shahryar County, depicted in Figure 5b. ADI was more successful in identifying PS points in man-made areas with stable targets than agricultural areas [70]. Since vegetated regions are the main land cover in our case study, an interpolation was applied to the obtained vertical map to extend the deformation information for the whole study area. The inverse distance weighted (IDW) interpolation was used with a weighting power of 2 and neighboring radius of 12 for calculating the vertical velocity interpolation.

**Figure 4.** The perpendicular baseline graph for the time-series data stack. The dots represents the images, and the edges denote the interferograms.

**Figure 5.** (**a**) The annual velocity map based on the PSInSAR technique. The map is superimposed on the Google Earth (GE) imagery; (**b**) IDW interpolation raster of the velocity map.
