Monitoring and modeling of a landslide in Kahroud, Iran, by InSAR measurements and slope stability analysis

This work demonstrates the potential of synthetic aperture radar interferometry (InSAR) based on interferometric point target analysis (IPTA) approach for landslide-induced deformation measurement. Twenty-eight Sentinel-1A images acquired between October 14, 2014, and October 27, 2016, are employed in order to measure the surface deformation associated with Kahroud landslide, northern Iran. Due to the lack of geodetic measurements and thanks to capability of processing small baseline highly coherent interferograms in the study area, the conventional small baseline subset (SBAS) time series analysis is exploited for investigation of the IPTA performance. The cross-comparison indicates some unwrapping errors in the IPTA results and removal of a large number of point target candidates due to the remained nuisance terms in the phase measurements, the deviation of the deformation behavior from a simple linear model employed by IPTA and the temporal gap in the available dataset. Two different sliding mechanisms, however, are revealed based on the SBAS time series analysis results. In the second part of the research, the slope stability analysis to characterize the slope dynamics is finally performed based on the deformation field extracted from two ascending and descending single interferograms. Accordingly, the soil shear strength parameters as well as a 3D model of the sliding surface are estimated followed by the calculation of the soil mass involved by the Kahroud landslide.


Introduction
Landslide is a downhill intermittent sliding of the rock, debris or earth occurred usually as a consequence of anthropogenic or natural activities (Bhattarai et al. 2005). Landslides are one of the most destructive and active processes that lead to erosion and changes in the land surface (Guzzetti et al. 2005). Iran, with its predominantly mountainous topography, highly active tectonic plates and different climatic and geological conditions, is prone to landslide hazards. One of the slope instabilities have occurred in the northern heights of 1 3 Kahroud village located in the central Alborz Mountains, northern Iran. According to the previous studies, Kahroud landslide has been active since, at least, 15 years ago (Dehghani 2016). Downhill movement of rock and soil materials due to Kahroud landslide would potentially cause to block Haraz road, Haraz River and tough damage to Kahroud village.
In order to investigate the landslide behavior, one requires to provide extensive spatial and temporal displacement measurements that could be costly and time-consuming, particularly in regional scales. Compared to other field-based geodetic techniques such as leveling and campaign Global Positioning System (GPS) surveys, interferometric synthetic aperture radar (InSAR) as a satellite-based approach is an effective tool for deformation monitoring from regional to local scales which can cost-effectively provide highprecision measurements. Since early 1990s, InSAR employing phase difference of two or more radar images over an area experiencing surface displacements has been recognized and used as an effective and well-established tool in the study of land surface deformation. InSAR-based technology in mountainous regions is widely applied for revealing and monitoring slow-moving landslides phenomena over the world (e.g., Ye et al. 2004). There are, however, significant limitations such as loss of correlation and inaccuracies in digital elevation model (DEM) which make conventional InSAR fails to measure the displacement specifically in a mountainous area with rough topography. Therefore, measuring displacement due to landslides which are happening in the surface slope would be significantly challenging.
Persistent scatterer InSAR (PSI) is proposed to address the decorrelation problem, residual topographic effect and atmospheric signal. Different PSI techniques have been used for landslide detection and its displacement monitoring (Wegmuller et al. 2008;Zhang et al. 2018;Aslan et al. 2020;Bekaert et al. 2020). In spite of the ability of PSI for monitoring the ground surface deformation in areas which are subject to decorrelation effects, some limitations may exist which highly affect the efficiency of the PSI methods. Hence, the performance of the existing methods which are almost freely available to the researchers' society have been carefully evaluated (Colesanti et al. 2003;Colesanti and Wasowski 2006;Dehghani 2016) for mapping and monitoring landslides. Few investigations, however, have been conducted for assessment of those approaches which are commercially accessible such as interferometric point target analysis (IPTA) (Werner et al., 2004).
In this study, the performance of IPTA as one of the PSI commercial techniques developed by Gamma Remote Sensing is investigated for mapping and monitoring the Kahroud landslide behavior. This approach identifies the point target (PT) candidates based on amplitude information of SAR images. Moreover, it employs a pre-defined temporal deformation model for final PS selection and phase unwrapping. In previous study carried out by Dehghani (2016), a modified version of Stanford method for persistent scatterer (StaMPS) was employed for displacement monitoring of Kahroud landslide. The main differences between IPTA and StaMPS is the way they detect the stable pixels and unwrap the interferometric phase. The performance of StaMPS in PS detection in Kahroud landslide was acceptable as expected due to the application of phase analysis in PS identification step. However, the phase unwrapping failed because the Nyquist sampling criterion employed in the temporal unwrapping is not met. This problem was fixed by incorporating a couple of coherent interferograms which might not always be available. The first objective in the present study is to assess the IPTA performance for Kahroud landslide monitoring regarding PT selection as well as displacement measurement. Unfortunately, no geodetic measurements are available in Kahroud for the evaluation of the IPTA results. However, it should be noticed that a significant number of coherent small baseline interferograms can be generated over Kahroud landslide so as to use small baseline subset (SBAS) time 1 3 series analysis whose results can be considered as a reference for evaluation. The SBAS time series analysis approach proposed by Berardino et al. (2002) in which the residual topographic effect and atmospheric signal are carefully estimated is applied. A cross-comparison between the SBAS and IPTA results is finally performed.
In some studies, the InSAR results and geotechnical analysis are combined in order to characterize the landslide behavior (e.g., Refice et al. 2019;Bovenga et al. 2017). For better assessment of the Kahroud landslide, the slope stability is analyzed through the determination of sliding surface associated with Kahroud landslide and estimation of soil shear strength parameters including the soil friction angle and apparent cohesion intercept of soil particles. The shear strength parameters are usually determined based on laboratory or field measurements. In this paper, however, it is tried to extract and estimate these parameters based on InSAR measurements. For this purpose, two Sentinel-1 image pairs acquired from ascending and descending orbit during a specific time period are employed to decompose the line-of-sight (LOS) displacement into the vertical and horizontal directions. The results are then used in a general limit equilibrium method by assuming simple circular sliding surfaces. A circular sliding surface is often recognized as a common mechanism of failure. To do so, a series of analyses are conducted by assuming various radii and centers for the sliding surfaces. A combination of these two parameters provides a large number of possible sliding surfaces, and the most critical one can be identified by finding the one with the minimum factor of safety. Of course, as the slippage has occurred, the minimum factor of safety must be around 1.0, in addition to be consistent with the deformation profile already obtained by the InSAR measurements. To attain such a factor of safety, the soil shear strength parameters have been changed incrementally to achieve this requirement. Different sliding surfaces are finally combined in a Geographic Information System (GIS) to generate the three-dimensional (3D) sliding surface and calculate the soil mass involved with the landslide during the time period of interest.

Study area
Kahroud village is located in a valley in the northeast of Mount Damavand as a part of Alborz Mountain Range, northern Iran. The sliding zone and the location of the Kahroud village is shown in Fig. 1. The upward slope is gently moving toward the village. The location of head and toe area of the landslide is shown in Fig. 1. Figure 2 illustrates the geology map of the Kahroud area. The landslide area dominantly contains sedimentary rocks such as sandstone, siltstone, shale, carbonaceous shale, claystone, quartzite and conglomerate. As shown in Fig. 2, the upward slope above the Kahroud village is mostly subject to the rock fall (Q r.f ). The southern part of the area is mainly covered by lava. The main part of the village is located within the lava composed of trachyandesitic and trachytic flows. Two different anticlines, Barzadeh and Panjab, are also surrounding the study area.
Different earthquakes with the magnitude above Mw7.0 have been reported in the Alborz Mountains. The landslides occurring in the Alborz Mountains have been probably initiated by these earthquakes. The deformation and rock fall caused by Kahroud landslide have imposed significant damages to the Kahroud village, Haraz road and their infrastructures. The Haraz valley along which considerable human activities have been developed is potentially threatened by Kahroud and other similar landslides.

3
In order to interpret the landslide movements which will be further presented, the slope and aspect of the study area are shown in Fig. 3. The slope and aspect maps are extracted from the ALOS DEM with the resolution of 12 m. According to Fig. 3, the dominant slope in the landslide area is between 15° and 30° toward the east and southeast direction. A small part located in the middle area faces toward the south.

Data
In order to measure the deformation caused by Kahroud landslide, 28 Sentinel-1A single look complex (SLC) images acquired from track 130 spanning between October 14, 2014, and October 27, 2016, are applied. Due to the fact that the Kahroud landslide is located in a slope facing east, assuming the right looking geometry of the Sentinel-1 satellite, the ascending orbit is selected in order to reduce the foreshortening effect. All data are freely available from the Copernicus Open Access Hub (https:// scihub. coper nicus. eu/ dhus/#/ home). Figure 4 shows the dataset collection with respect to spatial and temporal baseline which is employed in the IPTA and conventional SBAS time series analysis approaches.  In the IPTA, the interferograms are generated using the image acquired in December 20, 2015, as a single master (Fig. 4a). The master image is selected so as to minimize the temporal and spatial decorrelation in the processed interferograms.
On the other hand, the interferograms characterized by small temporal and spatial baselines are computed to be applied in the conventional SBAS time series analysis (Fig. 4b).
The minimum and maximum amounts of the perpendicular baseline are 20 m and 110 m, respectively, while the temporal baseline lies within the range of 12 days and 130 days. Using these thresholds for the baselines and eliminating the interferograms which are highly influenced by the atmospheric signal and decorrelation noise, a total number of 80 interferograms are remained to be used in the SBAS time series analysis.
In addition to the time series analysis, it is also tried to decompose the LOS landslide deformation into horizontal and vertical components. To this end, two pairs of ascending and descending Sentinel-1 images covering a time period of interest are applied (Table 1). It should be noted that the ascending and descending interferograms which are used for decomposition should cover the same time interval in order to minimize any distortion originating from different time coverage.

Methods
In SAR interferometry, using a significant number of differential interferograms, we are able to generate the deformation time series used for investigating the temporal and spatial behavior of the landslide. The principles of two different time series analyses, IPTA and SBAS, which are employed over the Kahroud landslide are presented in the following sub-sections. This is followed by an introduction to the slope stability analysis served for estimation of the shear strength parameters and the sliding surface. The overall flow chart of the methodology is shown in Fig. 5.

Interferometric point target analysis (IPTA)
The first objective of this paper is to evaluate the performance of the IPTA for monitoring landslide lacking bright scatterers provided by man-made structures. The method starts with a stack of SLC images co-registered to a single-master image. A candidate list of stable points is then generated based on the amplitude information of the returned echo. The first criterion used for an initial selection of the PT candidates includes low temporal variability of the radar cross section (RCS). The temporal variability is calculated as the mean to standard deviation ratio , hereafter called MSR. and are the temporal average and the standard deviation of the backscattering coefficients. In the case of a stable PT, a value of MSR considerably above 1.0 is generally expected. The first group of candidate stable points are thus selected by application of a threshold (e.g., 1.45) to the values of MSR. The second PT candidates list is generated based on the spectral coherence information as the stable PT typically contains a high degree of coherence. The two lists are then merged in order to achieve the final PT candidates. From this step onward, all processing steps will be performed on the selected points in a vector data structure for more efficiency and reduction in the required storage. The differential single-master interferometric phase are next generated for every PT using shuttle radar topography mission (SRTM) DEM with the resolution of 30 m. The interferometric phase is typically a point data stack of complex valued which is wrapped in the initial run. The double difference phase which is generated by making the phase difference between one point and another as a spatial reference indicates the spatial (perpendicular) and temporal baselines dependence. Therefore, a two-dimensional (2D) regression analysis using a phase model as a function of spatial (perpendicular) and temporal baselines is carried out for estimation of the unknown parameters estimation, i.e., residual topographic component and deformation rate. The phase model shows a linear perpendicular baseline dependence for the residual topographic phase which is a result of DEM inaccuracies. Moreover, this model indicates a linear time dependence for deformation rates which differ between one point and the reference one. The slopes of the regression line in the perpendicular and temporal baselines dimensions correspond to the residual topographic and relative deformation rate components, respectively. The problem, however, is that phase measurements are stated in a wrapped form. Therefore, the correct phase ambiguities are found as one part of the optimization problem in the regression solution. The corrections are further made to the DEM and initial linear deformation rate.
Any discrepancy between the phase double difference and the regression model is an indication of nuisance terms mainly due to the atmospheric contribution, nonlinear deformation component and the phase noise. Thanks to different behaviors of these terms in space and time, using a combination of spatial and temporal filtering, we are able to suppress phase noise, estimate the nonlinear deformation and atmospheric phase components (Ferretti et al. 2001). Consequently, the phase time series are corrected for the nuisance terms. The PT list is finally filtered by applying a threshold on the phase standard deviation obtained from the regression analysis. This process is done in an iterative manner to refine the unknown parameters as well as to obtain the final PT list.

Small baseline subset (SBAS) time series analysis
One approach for deformation time series analysis is the small baseline subset (SBAS) in which the processed interferograms are characterized by small spatial and temporal baselines (Berardino et al. 2002;Dehghani et al. 2009). The main idea behind the SBAS time series analysis is to obtain the deformation corresponding to each acquisition date by inverting the interferograms. The first step after flattening and calibration of the interferograms is to estimate the residual topographic effects remaining in the differential interferograms due to inaccuracies of DEM. To do so, a deformation model which can simply considers a constant velocity for each pixel is taken into account. It generally leads to jointly estimation of the main trend of the deformation time series hereafter is called the low-pass (LP) component and the residual topographic effects. The estimated LP component of the deformation and topographic artifacts is then subtracted in the wrapped form from each interferogram in order to reduce the fringe rate in the interferograms so as to improve the de-noising and unwrapping steps. Further steps of time series analysis and interferogram inversion are implemented after adding back the LP component to the newly unwrapped residual interferograms.
It should be noted that the estimated phase is contaminated with the atmospheric signals which are presented as undesired fluctuations in the time series. It can be estimated through a combination of low-pass and high pass filters applied in the spatial and temporal domain, respectively (Ferretti et al. 2001). The atmospheric phase component is finally subtracted from the estimated phase values to obtain the deformation time series as well as the deformation rate.

Slope stability analysis
A complete load-deformation analysis, revealing the possibility of failure and also the time-dependent or permanent deformation, requires the soil mechanical properties governing the deformability and the strength, to be fully prescribed. Despite such analyses are granted by completeness, they are seldom performed to assess the stability of the ground. In contrast, stability analysis requiring only the shear strength parameters and by assuming that the soil behavior can be idealized by a rigid perfectly plastic model has been commonly performed (Drucker and Prager 1952;Chen and Liu 1990;Pietruszczak 2010). The limit equilibrium (LE) technique, beside others like the limit analysis, is one of the most commonly applied methods in stability analysis of slopes prone to sliding. One important application is to estimate the most critical sliding surface and the overall shape of a landslide. Such analyses are also made in order to back calculations. Cornforth (2005) explained the limit equilibrium method and its applications to a series of case studies mainly categorized as landslides. In some cases, the form of the sliding surface and the sliding materials properties are found by back analysis. The universally admitted Mohr-Coulomb failure criterion has been found to be successfully applicable for such problems. Ground surface deformation measured by InSAR can be potentially applied for slope stability analysis from which the sliding surface is estimated. In the limit equilibrium method, there are several techniques which are mainly based on the method of slices and analyses are made to equilibrate the driving and resisting forces, moments or both to achieve some measure of the factor of safety. Interested readers may refer to Cornforth (2005) or more theoretically, Atkinson (2007) or Harr (1966) to get some practical as well as theoretical idea of the limit equilibrium. The whole methodology of the geotechnical modeling used in the present study is shown in Fig. 5.
The slope stability analysis in a time period of interest is performed along every single transect covering the entire landslide area. This requires the amount of deformation projected along the transect direction which can be extracted from the deformation field recognized by InSAR. Two SAR images acquired from two different geometries, i.e., ascending and descending, are used to extract the 2D surface deformation field, assuming low sensitivity of InSAR to the surface deformation along the north-south direction (e.g., Wright et al. 2004). It should be noticed that the stability analysis refers to the period covered by SAR data used for deformation field extraction.
Following the decomposition of LOS deformation into east-west and vertical directions, the deformation along every transect will be achieved using the azimuth angle of its alignment. In the first step of the geotechnical modeling, assuming some initial values for the soil shear strength parameters, the plastic limiting equilibrium analysis is performed in order to estimate an approximation of two-dimensional (2D) sliding surface in the underlying soil layers. In the next step, the factor of safety (F.S.) for each transect (along a presumed sliding surface) which is supposed to be 1.0 is estimated. The F.S. of 1 is an indication of the initiation of the landslide. One may note that higher values shows no movement, and lower values are technically meaningless. Using this assumption and soil shear strength parameters which are now estimated in a trial-anderror procedure, the final geometry of the 2D sliding surface for each transect is found. The F.S. of 1.0 as well as the consistency of the parameters between transects indicate the reliability of the inverse analysis in identifying the correct sliding surface and estimated soil shear strength parameters.
The 2D sliding surface corresponding to every single transect will be further employed in order to provide a three-dimensional (3D) model of the sliding surface. It is performed in a Geographic Information System (GIS) by interpolation of the surface elevation and the depth of each 2D sliding curve. The entire volume of sliding mass will be finally calculated by making the difference between the surface elevation and 3D model of the sliding surface.

IPTA and SBAS time series analysis results
To illustrate the practical applicability of IPTA, 28 Sentinel-1A SLC images spanning between October 14, 2014, and October 27, 2016, are applied. All SLCs are co-registered and resampled to the master geometry. Two different lists consisting of ~ 240,000 points are generated as initial PT candidates based on thresholds of MSR (1.45) and spectral coherence (0.32).
The differential interferometric phase for PT candidates is used in a 2D regression analysis. The maximum amounts of the deformation rate and height corrections required for the 2D regression are set based on the preliminary knowledge of the Kahroud landslide. The points with high phase standard deviations, which are estimated from the regression analysis, are eliminated from the initial list resulting in ~ 108,000 PTs remained. The nuisance terms such as atmospheric signals and phase noise estimated from the regression model residuals are removed from the interferometric phase in order to improve the unwrapping results in the next iterations. To prevent the nonlinear deformation component from leaking into the atmosphere estimation, a temporal filter of small size (~ 35 days) is employed. The 2D regression analysis and filtering out the points is repeated until convergence of the algorithm. The total number of the stable points remained in the last iteration is ~ 62,000 (Fig. 6). The number of PTs located in the landslide area is dramatically decreased after the first iteration. This is probably due to the complex temporal behavior of the Kahroud landslide. The applied 2D regression model assumes a simple linear time dependence for deformation rates. Large model residual at a point, however, is an indication of the deviation of its temporal deformation behavior from the pre-defined linear model. The PT candidates with large model residuals which are mainly located in the upper part of the landslide are eliminated from the list even if they are stable points. This would be the main drawback of the IPTA.
The maximum deformation rate detected in the landslide area is 0.06 m/yr which mainly occurs in the lower part of the deformation zone. This area is going far away from the satellite. Surprisingly, some parts of the landslide zone which is mostly located in the middle area is going toward the satellite with the maximum deformation rate of 0.07 m/yr. Assuming the slope of the hillside, this seems impossible at first glance. However, a horizontal motion of soil mass toward the west direction makes the LOS deformation behave as an uplift signal. This issue will be more explained while presenting the SBAS results. In addition to such interesting behavior, some phase jumps between the PTs are clearly observed which might be due to the remaining unwrapping error.
Due to the lack of geodetic measurements, evaluation of the IPTA results using the ground truth information is not possible. Instead, thanks to ability to process very high coherent small baseline interferograms in the study area, the conventional SBAS time series analysis could be applied for cross-comparison. In case of consistency between the results obtained from both methods, the high performance of the IPTA in landslide measurements would be confirmed.
To implement the conventional SBAS time series analysis, 80 interferograms characterized by small spatial and temporal baselines are processed. The deformation rate estimated from the time series analysis results are shown in Fig. 7. The mean coherence map was used to filter out the low-quality pixels using a threshold value of 0.7. Accordingly, two distinguished deformation zones with different behaviors are clearly observed in the landslide area. The middle part is moving toward the satellite, while the lower one demonstrates the negative values referring to the ground surface motion away from the satellite. Due to low spatial density of PTs and induced unwrapping errors, this behavior could not be totally recognized by IPTA except at a couple of locations. The LOS deformation rate ranges from −0.09 m/yr to 0.04 m/yr. In order to make the comparison between the IPTA and SBAS results more easily, the same colorbar is used for both Figs. 6 and 7.
For better interpretation of such a complex behavior, the deformation map is superimposed on a 3D model of the ground surface as shown in Fig. 8a. In this figure, the look direction to the 3D model of the Kahroud landslide is from the southeast. The soil mass generally moves along the direction of the maximum slope. According to the aspect map of the study area (Figs. 3b and 8c), the hillside in which the Kahroud landslide is occurring is generally a southeast-facing slope. However, this has also a westward component in the upper and middle parts of the landslide. This claim was proved by precise investigation of all of the unwrapped interferograms used for generation of the velocity map. In all of the interferograms shown as the supplementary data, this counterintuitive signal is clearly observed. Moreover, the interferograms generated from a different track (descending) proved the same westward component as will be further explained. Therefore, the Kahroud landslide contains two different mechanisms which are roughly shown in Fig. 8(b). The soil mass in the upper part (blue color) moves toward the southwest. Consequently, the landslide includes a horizontal component in the west direction. The combination of the west and vertical components of the landslide makes the ground surface move toward the satellite along the LOS direction. The lower part, however, has an insignificant westward component making the soil mass move toward the southeast, i.e., the highest slope direction. Complexities in soil formation and sedimentation process as well as historical geological processes, often cause such interesting behaviors, i.e., landslides with a complex mode of deformation. For further analysis, a Sentinel-2 imagery of the study area which is shown in Fig. 8d is employed. According to this image, an area separating two different landslide mechanisms is located in the boundary of two different geological units which are illustrated in different colors in the Sentinel-2 color-composite. While most landslides happen to occur along the apparent slope of the ground (as the main cause is the gravity effect), more complicated cases of sliding along other directions can take place. It can be explained from a geotechnical engineering standpoint as well as geological features of the subsurface layers when, say, the sliding surface is not necessarily parallel to the apparent ground slope and inclined at other directions. One other reason may be attributed to different materials in the sliding mass belonging to two different geological units when the landform is constituted. More or less closely related phenomenon to this case is the complex two-wedge sliding (Cornforth 2005;Holtz et al. 2023) where the mass is divided into two parts one pushing the other. In three dimensions, the upper wedge may slide to the sides rather than pushing the lower part downward. Other reasons like the complex uneven geological stress systems in the upper part of the crust can also cause complicated landslides. In Fig. 8e and f, we presented simple schematic 2D cases where the sliding surface, due to complex geological feature, is essentially in an opposite direction to the apparent (steeper) Steeper ground (primary slope) Fig. 8 a Elevation contours overlaid on the landslide rate map; b 3D view of the Kahroud landslide on which the deformation map superimposed. The arrows roughly indicate the direction of soil mass movement; c 3D view of the landslide aspect. The legend of the aspect map is given in Fig. 3(b); d 3D view of the Sentinel-2A image of the landslide area with the color-composite of R: B11, G: B8 and B: B2. The Kahroud landslide is depicted as a black polygon. In all three 3D view, the look direction is from the southeast of the area. e Schematic view of a possible sliding surface in the opposite side of a steep ground due to subsurface geological features f Two different geological units with different properties with the one having a mild inclination, more prone to sliding ground slope. These are two extreme possible cases which may occur in three dimensions. In general, more investigations can be carried out through field observations to further understand the causes and mechanism. Figure 9 compares the IPTA and SBAS time series deformation at a selection of PTs the locations of which are shown in Fig. 6. The root-mean-square errors (RMSEs) calculated from both results at points P1 to P5 are 0.005 m, 0.005 m, 0.025 m, 0.025 m and 0.04 m, respectively. According to the deformation time series, a 3-month temporal gap (from August 22, 2015, to November 26, 2015 exists in the available dataset. Before the temporal gap, two time series analysis results matched significantly in all points except for P5. However, an increasing discrepancy is observed after the temporal gap. One probable reason is that the deformation time series of IPTA are corrupted by unwrapping error after the temporal gap due to inaccurate estimation of the ambiguities. The performance of IPTA as other persistent scatterer interferometry methods is highly affected by poor temporal sampling of data. This will sometimes result in under-estimation of the deformation rate as observed in P1, P2, P3 and P5. The presence of the 3-month gap leads to over-estimation of the velocity in P4. More reliable conclusions about the under-or over-estimation, however, should be drawn based on ground truth measurements. Other reason which may be taken into account is the incomplete removal of the nuisance terms such as atmospheric signal as well as the nonlinear deformation component which leads to the remained fluctuations and hence the misfit between the phase observations and the regression model.
It should be noted that in spite of disagreement between two results after the temporal gap, the behaviors of nonlinear deformation are rather similar.
One other important limitation of IPTA is an assumption of a linear simple model for the deformation evolution. Any deviation from the linear behavior will result in considerable unwrapping error in addition to removal of PTs whose temporal behavior is not linear. This effect is observed in the deformation time series of P4 whose temporal evolution according to the SBAS results is not linear. Fig. 9 Comparison between the IPTA and SBAS deformation time series at a couple of PTs whose locations are shown in Fig. 7. a P1, b P2, c P3, d P4 and e P5 1 3

Slope stability analysis results
The last part of this section is devoted to the presentation of the results of the stability analysis which is conducted in a deformation period of interest. The main input to the geotechnical modeling is a 3D deformation field which can be extracted by a combination of two different radar geometries. Two coherent ascending and descending interferograms spanning a desired time period of interest are used. Figure 10 shows the LOS interferograms from two different orbits. A common time interval between ascending and descending interferograms is employed, and the deformation of both interferograms is linearly interpolated to the common interval.
By decomposition of the LOS deformation, the east-west and vertical components are estimated which is shown in Fig. 11. As already mentioned, there is a westward horizontal component located in the middle of the sliding zone with the maximum deformation value of 0.048 m. This component makes the corresponding pixels move along the LOS direction toward the satellite. The main part of the landslide area, however, contains an eastward In the next step, a series of back calculations have been made by geotechnical analyses to present a possible mechanism of a landslide by estimating the depth of the sliding surface, its shape and geo-mechanical properties. Such a back calculation is traditionally made in the absence of available subsurface exploration data and geotechnical soil profile (Cornforth 2005). To do so, ten parallel transects with grid spacing of 90 m covering the entire landslide zone (1300*450 m 2 ) are generated as shown in Fig. 12. The east-west deformation components at each grid point are projected along the corresponding transect using the azimuth angle of its alignment (~ 36°).
Every transect with the horizontal and vertical amounts of deformation at each grid point is employed in the limit equilibrium method of analysis using some initial estimates of soil shear strength parameters including the soil friction angle as well as some cohesion intercept. From a geotechnical engineering standpoint, the cohesion intercept cannot be fully mobilized at the same time as the friction angle. However, for some reasons, we assumed some slight amount of residual cohesion strength together with the frictional strength, to be more consistent with physical observations and to cover the requirements of a complete Mohr-Coulomb failure criterion (also accounting for small but nonzero tensile strength). These initial estimates are determined based on the preliminary knowledge of Fig. 12 Locations of transects superimposed on the DEM; the Kahroud landslide is depicted in a red polygon the area and a general view of the subsoil materials consisting mainly of a combination of coarse and some amount of fine materials. In addition, we formally assume that the friction angle is close to its residual value as the sliding has already occurred.
An approximation of 2D sliding surface is initially estimated along each transect. The more precise values of the soil shear strength parameters are then found in a trial-and-error fashion so as the minimum of the F.S. in all analyses for a particular transect approaches 1.0 ( Table 2). The consistency between the estimated parameters at all transects demonstrates the high reliability of the inverse analysis. These values are then used to find the final 2D form of the sliding surface. The results so obtained, for transects 2 and 7, are shown in Fig. 13. The depth of the sliding surface at each grid is shown in Fig. 13(b) and (d). As shown in these figures and similarly in other transects, the sliding surface depth is more significant in the lower part of the slope with higher deformation rate, which is not unexpected.
The depths of the sliding surfaces estimated at every grid point are further interpolated to generate a 3D sliding surface. Figure 14(a) shows the final sliding surface with contours indicating the sliding surface height. Accordingly, the area is separated into two main parts: (1) the upper part in which the sliding depth is close to zero and is presented as "shallow sliding." In this area, the sliding surface is located on the ground surface; (2) the lower part with considerable deep sliding surface and we call it a "deep sliding." Figure 14(b) shows the sliding depth which is calculated as a difference between DEM and the sliding surface interpolation result.
One of the main applications of the estimated surface is to identify the volume of the sliding soil mass. Using the ground surface topography and the 3D model of the sliding surface, the volume of the soil mass involved by the Kahroud landslide is estimated to be ~ 11,500,000 m 3 which is mostly associate with the lower part of the area. It should be noticed, however, that this amount of soil mass corresponds to the deformation measured within 24 days.
One point to be taken into account is that the stability analysis performed in this study gives a rough estimation of the sliding surface. As stated earlier, other sliding surfaces may form with a more complex nature and geometry. Although the InSAR is able to provide high-quality measurements of the surface deformation, the method used for landslide modeling and sliding surface estimation requires precise knowledge about the  soil properties which are not available in this study. Such a precise information may be achieved through the geotechnical site investigations, though, in case of a fairly uniform soil profile in depth, a circular sliding surface often suffices and shows a reasonable agreement with real observations. Therefore, in spite of the possibly low accuracy of the preliminary results so obtained, it is expected that the results present a valuable overview of the landslide which can be employed for practical applications.

Conclusions
In this study, the performance of IPTA as a persistent scatterer interferometry method was investigated over Kahroud landslide with complex deformation behavior. The criterion of our judgment was simply the results obtained from the conventional SBAS time series analysis. The results obtained from both methods were considerably matched at some points before the temporal gap existing in the available dataset. In some other points, however, the IPTA results seemed to be highly corrupted by unwrapping error which is an indication of large residuals of the regression model. Although the phase measurements were iteratively compensated for the nuisance terms such as atmospheric and nonlinear deformation components, the remained fluctuations make such misfit between the measurements and the 2D regression model as well as the incorrect ambiguities estimation. Moreover, due to deviation of the deformation behavior of PTs from a simple linear model employed by IPTA, a large number of PT candidates were removed after the first iteration which makes retrieval of the sliding pattern rather impossible. The deformation rate identified by conventional SBAS time series analysis revealed two distinct sliding mechanisms in the study area. The reason for such a complex behavior may be sought in the soil formation and sedimentation process. The soil shear strength parameters and 2D sliding surfaces at a couple of transects were estimated by slope stability analysis. The results were then integrated in order to form a 3D sliding surface including shallow and deep surfaces in the upper and lower part of the hillside, respectively. This was finally followed by the calculation of the volume of the soil mass involved by the landslide. The estimated amount of the sliding soil mass can be a cause of concern if the slope fails to maintain its stability and sudden movement ensues in the downhill direction. The results obtained from InSAR measurements as well as the slope stability analysis will be a key element for further geotechnical investigations as complementary studies over the Kahroud landslide.