Practical Estimation of Landslide Kinematics Using PSI Data

: Kinematics is a key component of a landslide hazard because landslides moving at similar rates can affect structures or collapse differently depending on their mechanisms. While a complete deﬁnition of landslide kinematics requires integrating surface and subsurface site investigation data, its practical estimate is usually based on 2D proﬁles of surface slope displacements. These can be now measured accurately using Persistent Scatterer InSAR (PSI), which exploits open access satellite imagery. Although 2D proﬁles of kinematic quantities are easy to retrieve, the efﬁcacy of possible descriptors and extraction strategies has not been systematically compared, especially for complex landslides. Large, slow rock slope deformations, characterized by low displacement rates (<50 mm/year) and spatial and temporal heterogeneities, are an excellent testing ground to explore the best approaches to exploit PSI data from Sentinel-1 for kinematic characterization. For three case studies, we extract proﬁles of different kinematic quantities using different strategies and evaluate them against ﬁeld data and simpliﬁed numerical modelling. We suggest that C-band PSI data allow for an effective appraisal of complex landslide kinematics, provided that the interpretation is (a) based on decomposed velocity vector descriptors, (b) extracted along critical proﬁles using interpolation techniques respectful of landslide heterogeneity, and (c) constrained by suitable model-based templates and ﬁeld data.


Introduction
The assessment of kinematics is a crucial step for the complete definition of the style of activity of deep-seated landslides [1], which includes their displacement rate, mechanism and degree of internal segmentation associated with the presence of nested bodies or actively deforming sectors [2]. In fact, landslides moving at a similar rate can cause different types and amounts of damage to elements at risk (e.g., settlements and linear infrastructures) and are often characterized by different collapse and runout potential [2]. Therefore, a robust quantification of the global geometry and spatial heterogeneity, which reflects the presence of differentially deforming slope sectors, is required to gain a better understanding of the mechanisms of landslide deformation and failure from a risk analysis and mitigation perspective.
A complete characterization of the geometry and mechanisms of deep-seated landslides requires integrating different sources of information, including borehole, geophysical, geotechnical monitoring and remote sensing data, integrated in 4D framework (i.e., 3D geometry and activity) [3][4][5][6]. Nevertheless, in preliminary appraisal studies of several landslide sites, or during initial stages of landslide investigation, subsurface site investigation data and monitoring networks are often unavailable, and landslide kinematics must be inferred by combining information provided by field evidence and surface displacement measurements provided by remote sensing technologies [2].
The last decade has seen considerable advances in the application of remote sensing techniques to the study of landslides and especially satellite monitoring system has gained The main characteristics of slow RSD (i.e., size, low rate, kinematic complexity) make them suitable to test the efficacy of a kinematics assessment approach based on 2D profiles.
In addition, this approach is often the only feasible one, since these phenomena are so deep that they undergo boreholes and geophysical investigations only in special conditions, e.g., when interacting with extremely valuable assets as hydroelectric power facilities [20,30].
In the past, slow rock slope deformation study was limited by the low spatial coverage and resolution of monitoring techniques, either unable to detect very small displacement rates [31] or to capture displacement patterns over space and time. PSI techniques, including PS-InSAR™ [8] and SqueeSAR™ [32], proved their capability to measure small ground deformations with millimetric precision, making them suitable for regional-scale landslide mapping and inventory studies [7,11,[33][34][35]. However, despite the improvements and refinements introduced by these methods, several limitations [34] still affect the application of DInSAR to landslide investigation.
A major limit of the technique is the inability of the satellite sensor to record the real 3D components of ground displacement, catching only the satellite line-of-sight LOS projection of any possible 3D ground deformation. Therefore, when the true deformation vector differs from the LOS, the sensitivity decreases and the interpretation of InSAR deformation measurements becomes challenging [10]. The application to the study of slow RSD is then further complicated because of low signal to noise ratio, presence of vegetation that can reduce the number of available PS and ambiguities related to the complexity and heterogeneity of landslide mechanisms [9]. The characterization of slow heterogeneous movements of large RSD has thus been commonly limited to recognize their degree of activity (mean velocity), and their potential to infer their kinematics is not obvious.
In this work, using PSI data derived from C-band Sentinel-1 A/B we evaluate the suitability of different descriptors and profile extraction approaches to infer landslide kinematics and outline spatial variations in movement rates in order to provide a best practice approach that can be applied also to other kind of phenomena.

Case Studies
We selected as test sites three slow RSD (Figure 2a) representative of different deformation pattern, displacement rates and potential geohazard impact. These are Mt. Solena (Val Fraele Figure 2b), Corna Rossa (Valfurva, [29,36], Figure 2c) and Farinaccio (Val Grosina [11,13], Figure 2d), all actively deforming at different rates between 5 and 13 mm/year with local peaks reaching values of 20 mm/year. They are located in the north-eastern sector of Valtellina (Figure 2a, North Italy) in the Austroalpine domain. Mt Farinaccio and Corna Rossa slopes are made of polydeformed metamorphic rocks, mainly paragneiss, micaschists and phyllites belonging to the Grosina-Tonale and Campo nappes, respectively. Mt Solena is carved in Mesozoic limestone and dolomite of the Ortles Nappe [37,38].
All the three RSD affect high energy relief slopes (>1000 m) with a mean slope inclination of 30 • affecting areas ranging between 1 and 10 km 2 . Their onset can be linked to the post-LGM destabilization of the valley flanks by means of progressive failure processes [23]. They are characterized by evident superficial morphostructures (e.g., Mt. Farinaccio), are part of ground based monitoring networks and affect critical scenarios (e.g., Mt. Solena impending over Cancano Lake and dam) and provide good examples of complex and segmented phenomena undergoing differential style of deformation (e.g., Corna Rossa) with nested sectors possibly evolving towards catastrophic collapse.
For each of them we performed a geomorphological mapping by means of detailed stereoscopic photo-interpretation of aerial imagery (Regione Lombardia TEM1   The main goal of this semi-detailed mapping [1] is to maximize the information at the slope scale thus providing a consistent dataset to be used as support for a kinematic and long-term activity analysis and give clues on the deformation and damage degree of each phenomenon. The integration between morphostructural information and PSI data is then necessary for a complete assessment of the style of activity of the landslides and provides a reading key to correctly interpret the ongoing displacement.

Workflow and Analyses
Starting from the state of the art [10,11,39] we suggest a methodological workflow that combines PSI data, morphostructural field evidence and simplified 2D finite element (FEM) models to retrieve a reliable description of the kinematics of the phenomena.
We perform the following analyses ( Figure 3): 1-PSI data post-processing to extract quantities suitable to describe landslide kinematics; 2-spatial analysis of (point-like) PSI data to select profile traces representative of landside complexity (swath profiles) and distribute the information over landslide area (interpolation); 3-extraction of 2D profiles and assessment of the most suitable approach; 4-interpretation of 2D profiles using non-specific templates derived from simplified 2D FEM models and comparison to field evidence of selected case studies. Steps followed for the identification of the best kinematic descriptors from PSI data.

PSI Datasets Post-Processing and Kinematic Descriptors
For our analysis, we used commercially available SqueeSAR™ datasets (TRE Altamira Table 1) derived from Sentinel1 A/B dataset acquired between 2015 and 2017 in both ascending and descending mode ( Figure 4). The simple examination of LOS velocity (V LOS ) is the most straightforward way to investigate the style of activity of a landslide, but for a correct interpretation, it is necessary to take into account the LOS parameters and the slope topography (slope, aspect, etc.) to estimate how much of the true 3D displacement vector can be observed [11,40].
Some authors tried to overcome this issue by approximating the movement to a surface-parallel displacement, as in the case of ice flows analysis [41], projecting the LOS velocity along slope (V slope [2,40,42]). This facilitates the interpretation of VLOS data and maximizes the data availability, but since it assumes a global translational sliding, it hampers any unconstrained interpretation of the landslide kinematics [11,43]. This is especially true for phenomena such as complex landslides, in which the internal displacement pattern can vary and differ from a simple slope-parallel movement.
The integration of InSAR displacement data from ascending and descending satellite orbits can help increasing the sensitivity for displacements close to the blind plane spanned by the two LOS vectors. However, as the azimuth displacement components cannot be retrieved by the sensor, the ascending and descending radar LOS directions are simplified as both belonging to the East-West plane (Figure 4c).
We thus combined the two datasets through a 2DInSAR decomposition approach [10,44,45] and extracted the displacement vector components (Figure 4c,f).
First, since shallow movement associated with PSI point-like data can be induced by the instability of slope deposits (e.g., scree, glacial, or periglacial features) rather than by deep-seated deformations [11,43], we removed points inside superficial deposits based on available geomorphological maps [1] refined by original mapping on Google Earth™ imagery and recent orthophotos.
Moreover, as different acquisition geometries see different objects and each InSAR data stack has dissimilar sets of measurement points [46], we divided the area into regular square cells (size: 25 m) and for each cell assigned to the cell centroid (Pseudo PS) the average LOS velocity of PS and DS from the same acquisition geometry.
Grid size was selected after parametrically testing different cell dimensions ( Figure 5). Very small grid cell size 10 m (Figure 5a), below the spatial resolution of the sensor, in this case Sentinel-1, about 20 m, results in a reduced probability of finding Pseudo-PS (i.e., cell centroids combining information from both ascending and descending data). At the same time wide cells (50-100 m, Figure 5c,d) may extend over sectors characterized by different degree of activity and kinematics (i.e., over nested landslides boundaries) mixing up different signals and resulting in a smooth slope response. We thus selected an intermediate grid size of 25 m (Figure 5b) since it provides a fair discretization of the slope and it almost coincides with the mean distance computed between the points of the ascending and descending datasets (~24.58 m).
For each gridded point (pseudo-PS) we extracted the vertical (V v ) and horizontal (V e ) components and the 2D total displacement vectors (V T ) [10]: where V a and V d are the ascending and descending LOS velocities (mm/year), and θ a and θ d are the incidence LOS angles for the considered satellite platform in the two acquisitions geometries. Depending on the failure mechanism, the superficial slope movement has different components along the landslide body and consequently V v , V e and τ change too ( Figure 6a). Close to the main headscarp, the displacement vectors have a downward movement [11] and the total vector plunges at high angle into the slope (Figure 6b). In the middle and lower part, the horizontal component tends to become dominant and T vector usually becomes parallel to the slope or points upward in response of the toe uplifting ( Figure 6b). Therefore, the displacement distributions may be used as a tool for interpreting the different geometry of the sliding surface for landslides of different typology and to identify active structures on the slope.
A first assessment of the local (i.e., cell-scale) slope kinematics can be inferred by observing the difference between τ and the local slope dip (α) in each square cell, namely ∆ (Figure 6a,b): Values of ∆ close to zero indicate slope-parallel sliding (Figure 6b,c), negative values indicate daylighting or bulging movements while, where ∆ is positive, the movement is mainly downward dipping in the slope (Figure 6b,d).

Data Spatialization
To prepare the extraction of 2D kinematic profiles, we used three different spatial analysis methods to highlight the main issues possibly arising from each different approach.
As first approach, through an original Matlab™ tool we generated longitudinal swath profiles, corresponding to segments perpendicular to the down-slope direction in which the statistics of the displacement rate and pseudo-PS geometric information are calculated. Swath profiles are spatial analysis tools commonly used in geomorphometry [11,47,48], but also already applied to the extraction of landslide velocity values along 2D slope profile traces. Frattini et al. (2018) already used a similar method considering swaths extending up to the lateral boundary of the RSD. However, lateral variations due to changes in displacement vector direction with respect to the LOS, presence of secondary landslides or active structures, the geometry of the failure surface and the physical mechanical characteristics of the material [11] can significantly affect the summary interpretation of global kinematics along the considered slope profile.
To overcome this limitation and maximize the accuracy of the data, we implemented our Matlab™ tool to independently set the orientation and geometrical parameters of the swath (width, sampling step size) in a handy and interactive way. The profile trace can be arbitrarily drawn inside the landslide polygon, it is then subdivided in stripes perpendicular to the trace direction with variable length and width and inside each of them the mean value of PS data is calculated (Figure 6a,b).
As second spatial analysis method, using the ArcGIS "Point Statistics" tool, we extracted for each PS and PseudoPS the average of the values falling within a specified neighborhood (5 cells). The final output is a raster where the value of each cell is function of the surrounding ones around that location (Figure 7c,d). This is the most restrained spatial interpolator as it considers only the displacement values averaged on the PS itself and its closest neighbors. However, if points are too scattered, the raster map results discontinuous and the derived profile is strongly influenced by outliers or isolated values, thus compromising a correct interpretation.
As the third method, we tested a Natural-neighbor interpolation (Figure 6e,f). This latter finds the closest subset of input samples to a query point and applies weights to them based on proportionate areas to interpolate values [49]. This method preserves input data values and produces a continuous surface except at the sample points. We computed a 20 × 20 m grid and interpolated a surface bounded between PS and PseudoPS locations. No extrapolation was used to approximate values outside the convex hull.

Profiles Extraction
Because of internal segmentation, LOS velocity values are not usually constant through the slope and also the kinematic interpretation may vary from sector to sector.
In the case of complex landslides such as slow RSD, swath approach can thus be used as first tool to explore the spatial segmentation and to identify the most suitable profile trace to unravel the deep kinematics.
Variable dimensions (Figure 8a,b) highlight either the local movement of the sector close to the profile trace (narrow swath) or a mean displacement trend of a broader area (wide swath).
Considering PS LOS velocity information: keeping constant the longitudinal length of the swaths (100 m), a narrow swath size (50 m) gives a more precise response than a wider one (2000 m) as it emphasizes local changes in velocity with sharp peaks (Figure 8b), possibly corresponding to active morphostructures. A wide swath on the contrary averages a larger number of points that may belong to sectors with different kinematics and, as result, returns a smoother trend (Figure 8b). If the landslide has a complex behavior and internal segmentation, the use of wide swaths can be misleading for the analysis of the landslide activity and kinematics, whereas if the landslide is homogeneous and without strong strain partitioning, a wide swath [11] may be used to interpret the general deformation pattern. From this perspective, a comparison of profiles resulting from different swath dimension can be used as first tool to assess the activity heterogeneity and select the most representative profile traces, respectful of the spatial complexity. Swath width should be calibrated based on field and morphostructural observations. Once a suitable profile trace has been selected, PSI data need to be interpolated over the landslide area in order to extract continuous profiles. Using swath profiles, a mean value is extracted inside each strip and its value can be very different from the one obtained through neighborhood statistics or natural neighbor interpolation in which weighted values are computed starting from the given interpolation points. Neighborhood statistics (point statistics) results strongly influenced by points distribution, emphasizing local response of isolated points and sharpening the values changes along slope in a more jagged profile than the one resulting from Natural Neighbor interpolation (Figure 8c,d).
In the following analyses we adopted natural neighbor interpolation to extract velocity and kinematic profiles and, to select the most reliable kinematic descriptor, we compared the effectiveness of possible geometric indicators (V v , V e , V LOS , τ, ∆) considering both synthetic 2D finite elements (2DFEM) models and real case studies with known kinematics.

DFEM Interpretation Templates
To get nonspecific reference templates representative of typical landslide kinematics, we performed simplified, 2D Finite-Element numerical simulations using the software RS2 (Rocscience Inc., Toronto, ON, Canada). Although based on a continuum small-strain formulation, the adopted code is able to account for deformation and failure mechanisms in both continuous and discontinuous rock masses [50,51].
For the simulations, we considered simplified slope geometries with constant slope gradient (30 • ) and characterized by imposed failure surfaces with different shape (translational, rotational, and compound) introduced as a Goodman joint element (pseudo-joints in continuum-based modeling) to constrain landslide kinematics. Slope height was set to 1200 m and failure surfaces traced at depth comprised between 200 and 400 m to simulate large RSD. These models do not cover the wide range of possible geometrical and mechani-cal conditions but are simply meant to extract the distribution of displacement components on simple failure surfaces.
The model domain was discretized into six noded triangular finite elements. Boundary conditions were assigned in terms of displacements (i.e., fixed bottom and side displacements), and a gravitational stress field was initializedWe attributed to the sliding mass values of strength and deformability parameters representative of common rock types (Table 2) as gneiss, schists, granitoids, plus idealized "very stiff" and "weak" rocks in order to extract a generic displacement signature for each kinematics. The non-deformed stable slope was constrained imposing an elastic behavior and high strength parameters. We considered homogeneous materials characterized by an elasto-plastic behavior according to a Mohr-Coulomb failure criterion and we ran the simulations using the SSR (Shear Strength Reduction) technique, which allows to evaluate the Strength Reduction Factor associated with computed stress-displacement fields and failure mechanisms [52]. We then considered for each model the best SRF stage displaying the most evident kinematic deformation style and the critical stability conditions. For each model, vertical (V v ), horizontal (V e ) and 2D total displacement inclination (τ) values were extracted along slope and plotted in normalized distance-displacement graphs ( Figure 9).
A constant decrease in vertical values (Figure 9) from the top to the toe of the slope is typical of rotational kinematics as the sliding surface is steep in the upper slope sector and then becomes progressively parallel or gently dipping into the slope. A perfect rotational kinematics may also present daylight τ values at the toe, corresponding to bulging induced by rock mass push.
On the contrary, translational and compound mechanisms are characterized by almost constant vertical values ( Figure 9) that tend to stabilize in a plateau that becomes more evident for rigid rocks. Localized higher vertical components correspond to the headscarp sector as clearly shown in the compound mechanism (Figure 9) where the most of the deformation is accommodated in an active wedge at the top of the slope connected through antithetic structures to a sliding sector (Figure 1c). It must be noticed that the vertical component can be strongly influenced by local conditions or by the rheological influence of the model as well as by differential internal deformation due to the mesh properties or geometric setting. Similarly, the horizontal component V e (Figure 9) is less meaningful because it is strongly biased by the geometry of the model and does not provide clear signatures of different kinematic styles.
A flat pattern, corresponding to a constant dip angle, is typical of translational landslides (Figure 9) where the sliding movement is manly parallel to the slope, while for rotational kinematics it presents a highly decreasing angle going from the main headscarp (top) to the toe. Compound mechanism (Figure 9) is a combination of the previous two, as it shows highly plunging vectors in the headscarp sector, followed by translational movement along slope, mirrored by almost constant τ values. This suggests that, as τ, also ∆ parameter can be considered a valid indicator of slope kinematics since it is directly linked to τ value and the local slope inclination. (Equation (5)) and gives a more accurate insight on the local deformation pattern.

Interpretation Using Geomorphological Mapping
We compared the 2D kinematics profiles extracted from 1D V LOS distribution with geometrical information of V v and τ resulting from the 2DInSAR analysis. We did not consider the horizontal displacement rate V e since it represents the horizontal movement on a 2D E-W vertical plane and its values can be highly biased on slope with unfavorable orientation ( Figure 10). On the contrary the vertical velocity represents the upward/downward displacement rate and gives more significant information on the landslide gravitational movement.
The profiles extracted for the Mt. Solena slope (Figure 10a-c) reflect a drop in V v values (Figure 10g) close to the headscarp, that going downslope tend to stabilize around a steady value of about −3.5 mm/y, suggesting that the entire mass is uniformly moving and there are no active structures that induce important vertical movements. Similarly, τ profile (Figure 10g) confirms this observation with only few local fluctuations in the top slope sector and at the toe. Instead, LOS velocity profile (Figure 10g) sharply decreases in the upper slope portion to then rise towards less negative values towards the toe. The simple interpretation of 1D LOS profile may thus result misleading since the change in LOS velocity is not directly reflected by a change in kinematic style and a rise in V LOS may be improperly interpreted as signature of a rotational movement, while only corresponding to a faster sliding sector. Mt. Farinaccio (Figure 10d-f), shows very similar V v and LOS curves because the slope has an unfavorable orientation, almost N-S, with sliding direction towards S and E-W component almost null. As consequence, since in the 2DInSAR approach the N-S component is set to zero and the V e is very small, LOS velocity and the vertical component tend to coincide. τ profile can be more easily interpreted since it shows dipping displacement vectors in the upper sector (Figure 10h) of the slope, corresponding to the main headscarp, which then decrease downslope as the movement becomes less steep, suggesting sliding on a surface that progressively becomes more parallel to the slope.
These considerations are also supported by the results of 2D FEM models that point out as τ (i.e., the inclination of the 2D displacement vector) is the most stable parameter to interpret the landslide kinematics since it describes geometric variations in the mechanics of the phenomenon more than bare velocity changes ( Figure 11).
As consequence, the related ∆ value, that corresponds to the difference between τ and the local slope angle (Equation (5)) can be exploited to investigate local surface variations (uplift, subsidence, etc.). In Monte Solena (Figure 11a), ∆ values have an almost constant value and the 2D displacement vector remains almost parallel to the slope (Figure 11b), as in the case of a modelled simple planar surface in which the type of movement is strongly limited (Figure 11c).
Changes in the profile trend correspond to mapped areas of debris accumulation and nested shallower phenomena.
A different scenario can be depicted for Mt. Farinaccio (Figure 11d). ∆ shows highly dipping angles in the upper sector (Figure 11e) of the slope corresponding to the main headscarp and then decreases downslope, but always remaining steeply plunging into the slope. This trend is typical of rotational phenomena (Figure 11f) and local fluctuations are due to the presence of active morpho-structures highlighted in the mapping.
Corna Rossa (Figure 11g) is a more complex phenomenon [11,29,36] with a strong internal segmentation as anticipated in Figure 7 by swath profiles of increasing lateral length. This is due to a strong structural control that induce the onset of contiguous sectors actively deforming with different mechanisms expressed by distinctive morphostructural associations of scarps and counterscarps. The westernmost sector is characterized by an association of scarps in the upper sector and counterscarps and the relative ∆ profile shows a double pattern: first it decreases in the headscarp sector and then it stabilizes on a slightly decreasing trend (Figure 11h) from 2000 m.a.s.l. to the valley floor. Comparing this pattern with the 2DFEM template, we can describe the kinematic as compound with a rotational uppermost wedge followed by a mainly sliding part (Figure 11i).

Discussion
A robust kinematic characterization is fundamental to understand landslide mechanisms, interpret their controls, and predict their potential evolution to plan appropriate risk mitigation strategies. Despite having the same displacement rate, different landslides may in fact behave in different ways both in an evolutive and risk perspective according to their deformation style, involved rock mass volumes, interaction with at risk elements and collapse potential [53,54]. To this aim, PSI data have been widely applied to extract an activity and kinematic analysis [9][10][11], but an in-depth investigation of the most suitable geometric descriptor has never been conducted.
In this study, we exploited the use of commercial [8,33] PSI products to investigate their effectiveness in unravelling the kinematics of slow RSD that are among the most complex landslides and can be thus considered the most complete test cases on which validate the analyses [55].
We propose a practical workflow to: (a) identify the most suitable, unbiased descriptors of landslide kinematics; (b) maximize the potential of sparse PSI data to obtain continuous 2D profiles; (c) obtain constraints on the kinematic interpretation of 2D profiles.
First, to retrieve the real kinematic behavior it is necessary to refine PS datasets to remove those points lying in the superficial debris cover or periglacial forms that can have differential movements and introduce false signals in the analysis. Then, it must be taken into account that the near polar orbit of the satellite causes the impossibility to obtain readings in the NS direction, which in some case, such as Mt. Farinaccio RSD, is the principal direction of sliding, essentially impeding the analysis of horizontal velocity from PSInSAR™ measurements. However, exploiting these data through a 2DInSAR approach, thus combining the ascending and descending datasets in the same time span and with a high spatial density, it is possible to increase the sensitivity of the measure and extract geometrical information associated with the 2D displacement vector.
Using appropriate spatial analysis approaches, point-like information can be then visualized and investigated through along slope profiles arbitrarily traced. In particular, using swaths with variable dimensions it is possible to assess a first degree of the landslide heterogeneity and select the most representative profile trace. Afterward, the use of a Natura Neighbor interpolation allows to obtain a spatial coverage of the information over the landslide area, yet preserving the original value of the PS or PseudoPS point, and to extract continuous profiles.
Their interpretation was then supported by nonspecific templates from 2DFEM models and real case investigations, that proved the effectiveness of profiles in describing the landslide kinematics. The integration of these two approaches first gives an insight on the most suitable kinematic descriptors, that are found in the inclination of the 2D displacement vector (τ) and the associated ∆ value, corrected for the local slope angle. Furthermore, it allows to unravel both the general kinematic signature of the analyzed landslide, controlled by the rock rheology and the geometry of the sliding surface, and the presence of local spatial heterogeneities, induced by nested shallower phenomena, active morphostructures and highly damaged zones recognized from the mapping.
This approach can be thus applied to analyze even complex landslides with heterogeneous strain distribution. A peculiar example is provided by Corna Rossa (Figures 11g and  12) characterized by a strong structural control, mirrored by distinctive morphostructural associations and slope segmentation [29,36]. The NW slope sector is affected by several orders of scarps, while the SE area is dissected by steep scarps and counterscarps arranged in a graben system [29,36].
The transition between these two sectors is marked by the abrupt closing of the NW sector main scarp in correspondence of the graben system and by persistent scarps which are oriented towards NNW-SSE. Different features are ascribable to different deformation mechanisms and suggest a transition from a mainly sliding sector to a "spreading" one, characterized by dominant extension accommodated by symmetric and asymmetric graben structures [29,36]. This change in kinematic behavior, inferred from morpho-structural observations, is confirmed by τ and ∆ profiles (Figure 12b,c). From NW to SE we can recognize a transition from a compound sliding (Figure 12d,e) with 2D displacement vectors oriented almost parallel to the slope (Figure 12e), to a dominant rotation-rototranslation (Figure 12f,g) and a more complex kinematic (Figure 12h,i) strongly influenced by the presence of deep scarps and counterscarps arranged in a graben system (Figure 12i). τ and ∆ profiles well reflect the deep deformation style because, even if the LOS velocity is underestimated or biased, the combination of two different acquisition geometries provides a depictive inclination angle of the 2D displacement on the E-W vertical plane. The more the slope is favorably oriented, and the sliding occurs in the E-W direction, the more τ corresponds to the true inclination angle and the kinematic interpretation becomes straightforward. As we depart from the LOS plane, we can generally expect to underestimate the 2D displacement vector inclination, resulting in a less effective detection of rotational movements and the integration with morphostructural evidence becomes fundamental for a correct interpretation of the deformation style.
In general, when the cross-section is not aligned E-W and when the movement vectors are oblique, the E-W movement vectors are reprojected along the cross-section as if the maximum movement occurs along the steepest slope direction. However, if a landslide experiences oblique movements, when the vectors are reprojected along the cross-section, their horizontal component will be much more underestimated with respect to the vertical one [56]. As consequence, the interpretation of 1D V LOS profiles is not always straightforward in the assessment of kinematics because they can present peaks and fluctuations linked to heterogeneous velocity and isolated high values more than true kinematic transitions. Negative or positive peaks can outline active structures or nested phenomena that, despite having different displacement rates, keep the same deformation style (e.g., fast or slow sliding sectors).
By extracting 2D profiles of different descriptors, in principle suitable for kinematic characterization, we showed that a velocity analysis based on 1D LOS velocity values (V LOS ) is not suitable to unambiguously represent landslide kinematics, especially when with complex phenomena such as slow RSD. These are characterized by sectors with different activity, kinematics and heterogeneous strain fields, and the single use of V LOS values results less robust and partially ineffective in describing the response of each slope sector. At the same time, the kinematics also influences the percentage of movement that can be sensed along the LOS thus resulting in different representative LOS velocities that only capture part of the total displacement vector thus returning a partial velocity information. Instead, the combination of data from different geometries captured with different LOS increases the sensitivity to displacement and reduces the complexity related to interpretation of InSAR data [10].
A 2DInSAR analysis should therefore be preferred, using the derived τ angle in combination with digital elevation model to identify areas undergoing displacement into (∆ > 0) or out of slope (∆ < 0) [57]. A limiting side of this approach is that the resulting 2D displacement is "forced" to lay in the EW plane, affecting the estimated horizontal (EW) and vertical components of the 3D displacement vector. To correctly quantify the components of the 3D vector, additional information must be provided along other LOS directions such as additional satellite tracks, ad hoc UAVSAR acquisitions, GPS data, offset pixel tracking. Summing up, our analysis proved that even at local scale, commercial PSInSAR™ data combined with field observation and validated by means of 2DFEM models, can play an important role in interpreting the kinematics of landslides, but to be able to give exact estimates on the true inclination and magnitude of the 3D displacement vector, a combination of in situ instrumentations, detailed mapping of landforms and geological structures is needed [10].

Conclusions
Recognizing the kinematics of a landslide is a major task to characterize its evolution and predict its potential impact on infrastructures and elements at risk. However, the definition of a proper kinematic style becomes particularly difficult for complex slow RSD, and the study of this is hampered by the limited amount of available data and the low displacement rate (mm to cm per year).
In this study, we exploited commercial PSI data to identify the most suitable kinematic descriptor of slow RSD and we provided a general overview of the limitations and advantages of this technique in the study of this complex class of landslides.
We highlight how a single 1D LOS analysis is not sufficient to describe complex landslide kinematics and a multi geometry 2DInSAR approach should be preferred as it provides a powerful tool for visualizing local vector geometry. We compare different geometrical descriptors using data spatialization approaches and profile extraction techniques respectful of the landslide heterogeneity and we integrate their interpretation with morphostructural mapping and 2DFEM results, retrieving a final consistent kinematic characterization.
Author Contributions: Conceptualization, methodology, investigation and validation: C.C. and F.A.; formal analysis and writing-original draft preparation, C.C.; writing-review and editing, F.A.; supervision, funding acquisition, F.A. All authors have read and agreed to the published version of the manuscript.