Multi-temporal landslide activity investigation by spaceborne SAR interferometry: Polish Carpathians case study

The main goal of this research is the activity state verification of existing landslide inventory maps using Persistent Scatterer Interferometry (PSI). The study was conducted in Małopolskie municipality, a rural setting with a sparse urbanization 10 in Polish Flysch Carpathians. PSI have been applied using Synthetic Aperture Radar (SAR) data from ALOS PALSAR, and Sentinel 1A/B from different acquisition geometry (ascending and descending orbit) to increase PS coverage and overcome geometric effects due to layover and shadowing. The Line-Of-Sight PSI measurements were projected to the steepest slope, which allows to homogenize the results from diverse acquisition modes and to compare displacement velocities with different slope orientations. Additionally, landslide intensity (motion rate) and expected damages maps were generated and verified 15 during filed investigations. High correlation between PSI results and in-situ damage observations has been confirmed. Activity state and landslide-related expected damage map have been confirmed for 43 out of a total of 50 landslides investigated in the field. The short temporal baseline provided by Sentinel satellite 1A/B data allows increasing of the PS density significantly. The study substantiates usefulness of SAR based landslide activity monitoring for land use and land development, even in rural areas. 20

is the PSI matrix approach with diverse SAR sensors, where specific thresholding of the landslide velocity, acquired from specific PSI processing, are performed. A detailed description of the PSI-based matrix approach is presented in Cigna et al. (2013). Since 2014, when the ESA Copernicus Sentinel 1 satellite was launched, the availability of SAR data has changed for almost entire world. The first articles of using Sentinel-l SAR data for landslide monitoring have been presented by Monserrat et al. (2016) and Barra et al. (2016). Recently, Kalia (2018) applied 66 Sentinel acquisitions in a descending orbit mode for 70 classification of landslide activity on a regional scale using PSI at the Moselle Valley in Germany.
In reference to the abovementioned papers, the main objectives of the present study are: • Updating the pre-existing landslide inventory map in the area of Rożnów Lake in Poland from PSI results derived based on ALOS-PALSAR (2007, Sentinel 1A (2014-2016) and Sentinel 1A/B (2017) data; • Evaluating the effect of SAR geometry delivered from ascending and descending orbits from ALOS PALSAR and 75 Sentinel 1 and the sensitivity to measure deformation over the study area;

•
Evaluating the difference in landslide activity updated from three diverse data stacks, namely: L-band (ALOS), Cband with one satellite (Sentinel 1A, with a revisit interval of 12 days) and C-band with two satellites (Sentinel 1A and 1B with a revisit time of 6 days), respectively; • Creation of landslide intensity and expected damage maps based on PSI results. 80 Verification of the objectives and achieved results was done by field investigations. Landslide activity states for 43 out of a total of 50 landslides were investigated. Some examples of landslide activity evaluation inspected during the field work are presented.
This study is probably the first investigation of this type in the Carpathians region. Particularly, the study area is located in the catastrophic damage due to abundant rainfall and flooding in 1997 and 2010. This region is rural and, thanks to the breathtaking landscapes, the population has grown rapidly in recent decades.
Moreover, it is also attractive to tourists. Apart from scatterer residential areas, many compact tourist buildings can be found 100 in deforested areas. These areas are characterized by high elevations and drainage divides, which are characteristically affected by landslides. In 2010, catastrophic landslide activity damaged approximately 150 buildings and many roads, transmission lines, crops and orchards (Gorczyca et al., 2013). Based on Gorczyca et al. (2013), further landslide activity within the study area is likely, thus, it is important to monitor all landslide activity in this region.
The study area covers the surrounding hills of Rożnów Lake (Fig. 1). This figure shows the landslide distribution within the 105 study area as well their predefined activity states. The first preliminary results of PSI application in the region of Rożnów Lake have been provided by Wojciechowski et al. (2008), Perski et al. (2010) and Perski at al. (2011). Using ERS 1 images, the authors showed the potential of this technique for landslide movement detection in rural and challenging terrain for PSI analysis.

Geological and hydrological settings of the study area
From the geological point of view, the study area is located in the Outer Carpathians of the Magura Nappe, Silesian Nappe and Grybów Unit. These geological units are built from flysch rocks, such as clay, sand, gravel, sandstone and shale, of the Quaternary and Tertiary periods. The study area is located on the borderline between Beskid Wyspowy Mountain and 115 Carpathian Foothills. Beskid Wyspowy comprises low mountains and middle foothills with slopes ranging from 10o to 35o and elevations from 300 to 340 m. The Magura Nappe is highly susceptible to sliding, especially with respect to hieroglyphicbeds containing variegated shales (Burtan et al., 1991;Wójcik et al., 2015). Detailed visualization of diverse geological units is presented in Fig. 2 and the explanation can be found in the supplementary material. From the geological point of view, the study area is located in the Outer Carpathians of the Magura Nappe, Silesian Nappe and Grybów Unit. These geological units are built from flysch rocks, such as clay, sand, gravel, sandstone and shale, of the Quaternary and Tertiary periods. The study area is located on the borderline between Beskid Wyspowy Mountain and Carpathian Foothills. Beskid Wyspowy comprises low mountains and middle foothills with slopes ranging from 10o to 35o 125 and elevations from 300 to 340 m. The Magura Nappe is highly susceptible to sliding, especially with respect to hieroglyphic-beds containing variegated shales (Burtan et al., 1991;Wójcik et al., 2015). Detailed visualization of diverse geological units is presented in Fig. 2 and the explanation can be found in the supplementary material.

Methods
The methodology flowchart is presented in Fig. 3. The generation of landslide activity maps integrates the PSI data, the PSI post-processed delivered products and landslide inventory map. These three consecutive stages are described in the following 165 subsections.   Table 1 presents five different data stacks, which have been processed using PSI. Due to the limitation of traditional DInSAR 170 techniques related to degradation of radar signals resulting from atmospheric effects and temporal and geometrical decorrelation, the PSI approach has found widespread application in landslide studies (Perski et al., 2011;Bianchini et al., 2012). The PSI approach is based on exploration of stable natural or man-made reflectors which present stable backscattering signal over time. Actually, the PSI approach has been previously applied in monitoring Carpathian landslides (Perski et al., 2010(Perski et al., , 2011. However, due to the high temporal decorrelation, resulting from vegetative cover and short wavelength (X-band), 175 it was only partially successful (Perski et al., 2009;Perski et al., 2011). This is mostly related to low PS density due to the TerraSAR-X data application. Therefore, exploitation of C-band (Sentinel-1) and L-band (ALOS PALSAR) data can bring more advantages, especially in rural areas (Lu et al., 2018 In this work, we exploited the PSI approach introduced by Feretti et al. (2001). This approach utilized more than 20 images to separate diverse interferometric components which correspond to deformation components, atmospheric error or topographic error due to the difference between digital terrain modelling (DTM) used for topography estimation and real surface backscattering elevation. The PSI approach is based on the following steps: (1) interferogram generation with respect to a common master image, (2) PS candidate selection based on amplitude dispersion index, (3) first estimation of atmospheric 185 phase screen (atmospheric influence) and topographical and displacement components and (4) second estimation of interferometric components and final PS points selection. For precise description of the algorithms applied within each step, see Ferretti et al. (2000Ferretti et al. ( , 2001.

Radar data and PSI processing
The present approach was performed for ALOS and Sentinel data stacks independently. Some necessary pre-processing steps were applied for raw SAR data, namely, focusing, image cropping and compensating for zero Doppler centroid. SARscape® software was used for SAR data processing (Sahraoui et al., 2006). Numerous interferograms were created from single look complex (SLC) data in original resolution without applying any pixel averaging. The precise orbital files from the European Space Agency (ESA) were used and SRTM 3 sec digital elevation models3 sec with resolution of 90 m were applied for topographic phase removal.

PS post-processing phase 195
After PSI processing, all results for the five diverse data sets have been post-processed in order to retrieve the most adequate displacement information (see module 2 in Fig. 3). Subsequent post-analysis descriptions are provided in the following subsections.

PS suitability analysis
PSI measurements provide deformation information in the LOS direction (from the satellite to the ground as determined by 200 the incident and heading angles). Thus, it is impossible to retrieve 3D displacement from InSAR directly. Horizontal and vertical components of movement (assuming no N-S horizontal motion exists within the study area) can be retrieved by combining measurements from ascending and descending orbits. However, in mountainous areas, where PS coverage is low due to the geometrical conditions and distortions, velocity decomposition can be problematic. Therefore, conversion of LOS deformation into the most probable direction (direction of maximum slope), by assuming a pure translational movement 205 mechanism, is commonly used (Bianchini et al., 2012). In order to assess the sensitivity of InSAR LOS measurements in a specific point, we applied R index. This index represents the geometrical visibility of the specific area with respect to the SAR system used. R index takes into account morphology of the study area and the acquisition geometry of the particular SAR system and is calculated as ( Notti et al. (2010) ): where α is the angle of the satellite track from North, ϑ is the incidence angle, S is the slope and A stands for the aspect. The   R index close to zero shows low SAR sensitivity to measured displacements for specific locations. Based on testes for different test sites, Notti et al. (2014) have found out that distinguishing of four classes is the optimal solution. Therefore, R index has been divided into four classes, wherein we used the same interval length for all classes except the first. The first class (R index 220 <0) represents the location where detection of the PS points is difficult due to geometrical distortions. The Second class (R index 0 to +0.33) indicates the location where the slope geometry is not convenient for SAR measurements. The third group https://doi.org/10.5194/nhess-2020-112 Preprint. Discussion started: 17 April 2020 c Author(s) 2020. CC BY 4.0 License.
(R index 0.33 to 0.66) presents PS on slopes with acceptable geometry for SAR displacement measurements. The fourth group (R index 0.66 up to 1) represents high sensitivity of LOS deformation monitoring. This means that, in this location, the slope direction is almost perpendicular to the satellite LOS direction and SAR sensors can perfectly measure the deformation over 225 this location. Summarizing, in further studies, the PS points with R index higher than 0.33 were considered.

PS velocity projection along the steepest slope
Deformation retrieved from PSI is 1D measurement in the direction to the satellite. LOS deformations ( ) were projected along the steepest slope according to the equation: with β as the angle between the steepest slope and the LOS direction.
Despite the great advantage of the motion represented in the slope direction, this projection has some limitations. First, when β=90o, Vslope goes into infinity. Here we followed Herrera et al. (2013) and selected an absolute maximum value of β=72o, which is equivalent to cosβ=0.3, so that Vslope cannot be higher than 3.33 times than that of VLOS. In order to remove any Vslope exaggeration, we considered only PS points for which cosβ >0.3. 235 Moreover, we discarded PS points which show Vslope>0 because positives values represent uphill movements and it is not representative for small landslide movements, even though positive values exist within landslides, especially in the toe area.

Velocity thresholding for activity state estimation -PSI based matrix approach
For the activity and intensity assessment, some representative values have been retrieved from PSI analysis. These values are fixed values which are used for: (1) distinction between moving and non-moving landslides and (2) discrimination of extremely 240 slow from very slow moving landslides. It strictly depends on the study site considering the deformation processes and typology (Cigna et al., 2013). Commonly, the average of LOS velocity estimates is applied as representative of velocity (Bianchini et al., 2012). However, different thresholds are applied to assess the landslide activity. These values can vary with respect to SAR data wavelength and projection of the velocity. For VLOS, some authors (Righini et al., 2012;Herrera et al., 2013) utilized 2 mm/yr as the velocity threshold for landslide activity assessment for C-band data and 5 mm/yr for Vslope 245 (Cigna et al., 2013, Bianchini et al., 2013. These changes are mostly correlated with different velocity distribution patterns. For the LOS velocity, distribution is almost normal (Gaussian), while for SLOPE is second negatively skewed as a result of the PS reduction (Bianchini et al., 2013). Therefore, for activity state estimation, we applied 5 mm/yr as the Vslope threshold.
Based on the pre-existing landslide inventory map and PSI post-processed results, we applied the PSI-based matrix approach (Fig. 5). According to the PSI matrix (Fig. 5)

Landslide intensity estimation
Based on the representative values of deformation, landslide intensity was assessed relying on the Cruden and Varnes (1996) intensity scale. This scale is based on PSI velocities consisting of three categories: negligible, extremely slow and very slow ( Fig. 6). In reference to Cruden and Varnes (1996), landslides with sufficient information are divided into negligible (mean velocity 5 mm/yr), extremely slow (mean velocity between 5 and 16 mm/yr) and very slow (mean velocity between 16 mm/year 260 and 1.6 m/yr). It is also worth to mention, that some researchers use slightly lower thresholding (e.g. Bianchini et al., 2012) with the argumentation that PSI underestimates movements in comparison to the real movements. historical PSI estimated movements, respectively.

Landslide activity state and intensity map generation
The landslide activity and intensity maps integrate the post-processing derived products (R and Vslope) with the existing landslide inventory. For the assessment of landslide activity, the previously described PSI matrix-based approach was applied. 270 Results of the activity assessment are presented in Fig. 7. However, the activity state has been presented only for landslides where sufficient PS points have been found. At least four PS points within a landslide body were set up as the threshold. A number of landslides in the study area are complex landslides. Such landslides are divided into parts and represented in SOPO databases as separate objects. In this case, the mentioned threshold is related to each landslide object. After the post-processing phase, 3898, 5260 and 10,798 PS points were found within the landslide boundaries for ALOS PALSAR, Sentinel 1A and Sentinel 1A/B data stacks, respectively. This allowed us to estimate the activity state and intensity 285 scale for 128, 130 and 205 landslides using ALOS, Sentinel 1A and Sentinel 1A/B data, respectively. Fig. 9 presents, in the form of pie chart some statistics to landslide activity state assessment depicted in Fig. 7. The left column of pie charts in Fig. 9 corresponds to the maps shown in Fig. 7 (a), (c), (e) respectively. From the total number of landslides and landslides objects (according to SOPO database), 23%, 20% and 38% have been updated by PSI approach for ALOS, Sentinel-1A and Sentinel-1A/B data respectively. Percentages of particular activity states within updated landslides are 290 illustrated in the right column of pie charts in Fig. 9. These pie charts correspond to the maps in Fig 7 (b), (d), (f). Note, that the last maps show also activity state updated based on PSI-matrix. Therefore there are also landslides, which activity state has been updated based on historical or pre-existing data (SOPO database) if insufficient number of PS were detected on the landslide object (compare also Fig. 5).

Possible hazard assessment
Based on a literature review, a downstream investigation was performed and additional thresholds were set up in order to 300 assess possible hazards related to buildings and infrastructure located in landslide areas. For this purpose, we applied the method proposed by Mansour et al. (2011), i.e., the threshold of 10 to 100 mm/yr as a minimum landslide velocity which can cause moderate damage to infrastructure and buildings. Velocity rates higher than 100 mm/yr can cause major damage to infrastructure and buildings. Landslide with velocity below 10mm we classified as landslide with minor expected damages.
This thresholding has been adopted as an additional criterion in order to support environmental planning and management 305 strategies to areas which can be characterised by high landslide hazard and, consequently, should be addressed to potential damages protection. In Fig. 10, possible damages caused by mass movements in the study area are presented for three diverse PSI processing results.

Field validation
The confidence degree of landslide activity maps can be validated throughout their comparison with external independent sources, such as damage inventories, in situ measurements, field checks, etc. This evaluation aims at assessing whether 315 measured displacements represent landslide dynamics. It is worthwhile to emphasise that the reality of the PSI-velocity estimation, itself, is not assessed due to the lack of external field measurements, but whether this measurement is related to landslide activity or not. For this reason field verification was performed.
Due to the abundant number of landslide located within the study area (506), only 50 landslides, which are expected to produce moderate damage, were investigated in the field. Activity states have been confirmed for 43 landslides, which indicates around 320 86% success rate. However, even if no evidence was found in the field, it does not mean that the landslide is not active.
Landslide damage, especially in agricultural areas and to buildings, are immediately repaired and, therefore, not always visible to people living there. Some examples of field verification and descriptions of investigated landslides are presented in the following subsection. The field work was performed on November, 2018 and was led by an experienced landslide expert.

Landslide "Just-Tęgoborze," SOPO ID 23374 325
The landslide "Just-Tęgoborze" is a group of rock and debris slides occurring below St. Just pass. This landslide is located on the southern and eastern slopes. The landslide was developed on clay of various genesis and from the marginal Magura Nappe formations. The upper part is founded on the outcrops of the Magura sandstones and the lower part on the slate sandstones of the sub-Magura units. It has been documented as active for the last 50 years. The activity state is regularly confirmed by road damage. Essential movements occurred in June 2010 as a result of abundant rainfall. In addition, many residential buildings 330 have been damaged. Landslides tend to develop and increase activity over a large area. Presented landslides are very difficult to stabilise due to the activity, considerable thickness of the colluvium and slate shales occurrence in the ground (Wójcik et al., 2011). This landslide was assessed as active (in 2010), dormant (2014-2016)  speckled slates overfilled with thin sandstones (Paleocene-Eocene). In the south and southeastern directions, these settlements are also flat on the pieces, which belong to the Dukla unit. In the southern part of the landslide, there are thick-walled Cergowa sandstones (from Oligocene) and in the eastern part sandstones and shales of the Lower Krosno layers (Perski and Wójcik, 2010a;Wojciechowski et al., 2012).
In addition to the discontinuities associated with the slides, the geological structure of the landslide is marked by faults. The 350 Zbyszyce landslide has been active for at least 60 years. The communal road crossing the landslide is constantly destroyed in the middle and lower parts. This landslide was assessed as active (in 2010

Landslide "Lipie-Jelna" SOPO ID 73194
The Landslide "Lipie-Jelna" landslide is located between the towns of Lipie and Jelna. The landslide is a rock and debris slide. 360 The landslide material consists of clays, sandstones and shales from the Quaternary period. It is an old landslide with main scarps located more or less parallel to the extent of the landslides, ending with a 3 m landslide crown. The landslide area covers the whole area of the slope. When reactivated, this landslide covered about 80% of the landslide area. It tends to develop uphill and its activity depends on weather conditions. Many building cracks, as well as damage to roads and cultivated fields exist within the area of the landslide (Wójcik and Krawczyk, 2010). This landslide has been assessed as active based on all PSI 365 results with moderate damage caused by landslide activity. In Fig. 13, the extent of the Lipie-Jelna landslide is presented together with photographs taken during the field verification.