Shallow Landslides Physically Based Susceptibility Assessment Improvement Using InSAR. Case Study: Carpathian and Subcarpathian Prahova Valley, Romania

: A multi-temporal satellite radar interferometry technique is used for deriving the actual surface displacement patterns in a slope environment in Romania, in order to validate and improve a landslide susceptibility map. The probability the occurrence of future events is established using a deterministic approach based on a classical one-dimension inﬁnite slope stability model. The most important geotechnical parameters for slope failure in the proposed study area are cohesion, unit weight and friction angle, and the triggering factor is a rapid rise in groundwater table under wetting conditions. Employing a susceptibility analysis using the physically based model under completely saturated conditions proved to be the most suitable scenario for identifying unstable areas. The kinematic characteristics are assessed by the Small BAseline Subsets (SBAS) interferometry technique applied to C-band synthetic aperture radar (SAR) Sentinel-1 imagery. The analysis was carried out mainly for inhabited areas which present a better backscatter return. The validation revealed that more than 22% of the active landslides identiﬁed by InSAR were predicted as unstable areas by the inﬁnite slope model. We propose a reﬁnement of the susceptibility map using the InSAR results for unravelling the danger of the worst-case scenario.


Introduction
Landslides are among the most significant natural events that cause human and material losses worldwide. Conventionally, landslides are defined as the movement of a rock or soil mass along a slope [1]. In Romania, landslides are the most frequently occurring slope processes in the Cretaceous and Paleogene flysch sector of the Eastern Carpathians, and in the Subcarpathian region of Neogene molasse deposits e.g., [2][3][4][5][6]. The Carpathian and Subcarpathian Bend area (the Curvature Carpathians) is especially susceptible to landslides, amid intense local seismic activity coupled with neotectonics, both of which favor instability e.g., [7,8]. Furthermore, pluviometric regime specificities, i.e., alternating wet and dry cycles, play a role in setting off shallow planar and translational slides [6,7,9]. According to Surdeanu [9], over the past 130 years, four pluviometric excess periods (1912-1913; 1939-1942; 1970-1972 and 2004-2005) induced most of slope failures in Romania's flysch and molasse areas. Most rainfall-triggered landslides are shallow and directly linked to subsurface hydrological processes. Given the broad countrywide extent of the phenomena, in the 1980s, Romania joined the global cluster of countries where critical prediction studies are conducted on shallow rainfall-induced landslides. Most predictions of rainfall-induced landslides are presently based on landslide susceptibility maps and on the correlation between landslide occurrence and rainfall thresholds [10,11]. Landslide susceptibility maps, validated by direct field observations, are crucial elements of prevention studies and risk management applications for landslide-prone areas. The traditional technique we employ is the one-dimension infinite slope stability model, the oldest and most physically (geotechnically) simple model, based on the limit equilibrium theory, which is used to ascertain a given slope's safety factor [11,[14][15][16][17][18][19][20][21]. The safety factor (SF) is computed as the dimensionless ratio between the main forces that act on a slope: shear stress, which provokes motion along a potential failure surface, and shear The traditional technique we employ is the one-dimension infinite slope stability model, the oldest and most physically (geotechnically) simple model, based on the limit equilibrium theory, which is used to ascertain a given slope's safety factor [11,[14][15][16][17][18][19][20][21]. The safety factor (SF) is computed as the dimensionless ratio between the main forces that act on a slope: shear stress, which provokes motion along a potential failure surface, and shear strength (described by the Mohr-Coulomb approach). In the assessment of the infinite slope model, it is assumed that the slip plane and water table are parallel to the slope's topographic surface [21,22]. The infinite slope model is suitable for shallow sliding with planar failure surfaces, where the sliding depth is considerably lower than the slope's length [23]. Shallow land-Remote Sens. 2021, 13, 2385 3 of 22 slides are common translational slips and are characterized by a thin soil cover (usually no more than a few meters [24]). These slope movements dominate landslide processes in hilly areas, as evinced by numerous specialized studies e.g., [25]. The Carpathian and Subcarpathian Prahova Valley, the proposed study area, is a mountainous and hillslope environment that is highly fitting for testing the infinite slope model's robustness, as well as for determining its limitations (Figure 1).
In order to enhance the assessment of landslide susceptibility maps across extensive areas, the infinite slope model can be integrated in Geographic Information System environments (GIS), as this is a pixel-based calculation that considers each raster cell separately e.g., [26][27][28][29][30][31][32][33][34][35]. However, a variable degree of uncertainty is linked to the complex and inherent variability of geological features when applying physical models over large areas. The validation process is thus a necessary phase of determining the extent to which the result reflects an accurate image of the real settings according to available field data. The present study tests the landslide susceptibility map against field survey landslide inventory database entries alongside interferometric synthetic aperture radar (InSAR) deformation maps, in a statistically oriented manner.
InSAR is a satellite monitoring method that enables the first-ever capability to effortlessly derive surface displacement over extensive areas at a low cost [36]. The main condition for detection is for the displacement to affect radar coherent targets, such as man-made features, open rocks or land [37,38]. An additional noteworthy benefit of this technique consists of enabling users to derive displacement rates of high vertical accuracy (sub-millimeter) [39], with a relatively high temporal resolution of up to 6 days in case of Sentinel-1 A/B satellites above Europe and global availability of free satellite data [40]. InSAR monitoring is also able to deliver a different outlook over a given study area, either by identifying new, undiscovered developing landslides, or by validating others that had been identified previously [37]. What sets InSAR deformation maps apart from susceptibility maps generated using conventional methods is the kinematic characterization of landslides, which is useful for validating and bringing pre-existing landslide inventories up to date. In addition, dynamic characteristics can be used as an input for mechanical models describing the dynamics of landslide processes.
In the last two decades, several multi-temporal interferometric methods have been successfully applied for updating and validating landslide inventories by combining traditional geomorphological information with observations of ground displacement analysis, from which the most recurring were the Persistent Scatterers interferometry [41][42][43][44][45][46] and the Small BAseline Subset interferometry [47][48][49][50][51][52]. Due to its enhanced capabilities, the InSAR multi-temporal technique has been used in a multitude of studies for landslide detection. Several satellite data are normally used for landslides studies as the satellite imagery acquired in C-band, such as historical satellites ERS and Envisat [46,[52][53][54][55], or the more recent data from satellites Radarsat and Sentinel-1 [45,[56][57][58], X-band satellites-TerraSAR-X [44,59] and L-band satellites-ALOS PALSAR [48,54]. These studies use satellite interferometry as the main tool for assessing slope instability, and usually compare these data with geodetic measurements, field recognition or optical imagery for validation. Righini et al. [46] use deformation maps obtained with satellite interferometry to update a pre-existing landslide inventory in Biferno Basin, Italy. Bardi et al. [59] use InSAR techniques to detect slope instability and combine the InSAR results with independent statistical analysis and activity models to characterize the activity state of the earth surface in specific regions.
In Romania, the use of satellite interferometry for studying landslides is recent; currently the only study is by Necula [60], which focuses on the detection of active landslides in a habited area from Ias , i county using multi-temporal techniques applied for Sentinel-1 imagery.
In this study, we propose applying the SBAS interferometric technique for detecting slope instability. Additionally, this technique is also used as a complementary tool to quantify its performance in terms of predicting landslide-prone areas compared to the Remote Sens. 2021, 13, 2385 4 of 22 infinite slope model. We propose a comparison between the susceptibility maps derived using the classical one-dimension infinite slope stability model and the InSAR deformation map. Because of InSAR characteristics, the study area focuses on inhabited communities which present coherent radar targets. As data sources, we selected the C-band Sentinel-1A and B European Space Agency (ESA) satellites. These satellites were considered suitable for the current study due to their free availability and data continuity since 2014, as well as their enhanced temporal resolution. For data processing we use the Berardino [47] SBAS algorithm due to its ability to enhance coherence, reduced noise of interferograms, and to identify nonlinear movements.

Tectonic and Geological Settings of the Study Area
The study area consists of the Carpathian and Subcarpathian Prahova Valley, which separates the Eastern and Southern Carpathians over 56 km (Figure 1). The Prahova Valley is among the most densely populated sectors of Romania's Carpathian and Subcarpathian region, as well as being the location of an important trans-Carpathian roadway connecting Balkan countries and Western Europe [8,35]. Over the past century, a wide range of infrastructure projects have been implemented along the Prahova Valley, including roads, railways, bridges, watercourse regulation and bank protection measures, as well as oil and in-channel gravel extraction activities, all of which changed the local natural environment to varying degrees [61].
The Carpathian stretch of Prahova river consists of a series of narrow reaches spaced out by minor tectonic-erosive depressions. The mountain valley slopes feature four fluvial terraces. Diluvial slope formations primarily consist of sandstone fragments and limestone, cemented by a clay-silt matrix (the flysch facies) that determines slope instability notably in areas with no abundant alpine vegetation. A 2 m diluvium layer prevails, below which a highly fractured and altered flysch facies is present, down to a depth of approximately 5 m. Flysch deposits are folded in serried synclines and anticlines. Deep-seated landslides and rock falls spread on the hanging wall of reverse faults. Most of these slides are of Pleistocene age, stabilized under forested conditions [9,[62][63][64][65].
One of the defining elements of the present study area is the large Subcarpathian erosive basin spread between the Comarnic-Breaza and Câmpina settlements in which Prahova river carved four terraces [8,14]. The erosional river valley is lowered in the northern part, in the Barremian-Aptian flysch zone. This northern part is represented by a marl-calcareous flysch with thin micro-breccia layers, known as the Comarnic formation, flanked on its southern edge by the Senonian red marls (Gura Beliei formation) and a Paleogene marl-sandstone flysch formation (the Sotrile flysch). Downstream, the confluence of the Belia and Prahova rivers outcrops the post-tectonic bedcover [8,63]. From the locality of Breaza to Cornu, the Breaza syncline extends in a west-to-east direction (Figure 1), has developed symmetrical flanks and is axially faulted. The structure's southern sector is elevated along the Breaza fault [8,14], where is represented by lower-Miocene deposits, such as the Cornu strata and the Brebu conglomerates. In the southern part of the Cornu locality, the molassic deposits are present with the Miocene-Pliocene-Lower Pleistocene strata, and they form the Miocene-Pleistocene zone ( Figure 2).
In general, shallow-to-medium deep landslides are quite frequent in the Subcarpathian area. The active processes, affecting over 60% of scarp surfaces, are triggered by intense or prolonged rainfall periods [8,14]. Rainfall-induced landslide reactivation also affects bodies and scarps of older landslides, and the slopes where the land is used for intensive grazing. In this study, we use a detailed landslide inventory of over 400 landslides (Figure 1) for model evaluation along with geologic and morphologic information. The landslide inventory resulted from repeated research campaigns carried out by Prof. Armas , in the study area over a 5 years period. A detailed description of the environmental conditions, landslide types, forcings (precipitation, groundwater, earthquakes, river erosion, human impact), motion and evolution pattern can be found in previous work of the author [8,61,[66][67][68]. From the landslide database, only shallow, slow-moving slides are selected for this study. landslides ( Figure 1) for model evaluation along with geologic and morphologic information. The landslide inventory resulted from repeated research campaigns carried out by Prof. Armaș in the study area over a 5 years period. A detailed description of the environmental conditions, landslide types, forcings (precipitation, groundwater, earthquakes, river erosion, human impact), motion and evolution pattern can be found in previous work of the author [8,61,[66][67][68]. From the landslide database, only shallow, slow-moving slides are selected for this study.

Spatial Database and Input Parameters
The deterministic approach necessitates a high volume of data, e.g., a high-resolution digital elevation model (DEM) and geomorphological specificities (slope angles, sine and cosine), existing thematic maps, a series of detailed field surveys, geomorphological

Spatial Database and Input Parameters
The deterministic approach necessitates a high volume of data, e.g., a high-resolution digital elevation model (DEM) and geomorphological specificities (slope angles, sine and cosine), existing thematic maps, a series of detailed field surveys, geomorphological investigations and mapping procedures, and numerous laboratory tests. Our DEM is obtained by the automatic interpolation of 5 and 10 m contour lines derived from existing 1:25,000 and 1:5000 topographic maps, enhanced with data from a 1 m resolution digital elevation model achieved by airborne LiDAR, specifically performed in 2012 along the valley bottom [35]. The geotechnical input parameters such as cohesion and friction angle of slope materials are obtained from laboratory analysis. The data we use are part of a larger 2014 geotechnical study for the Comarnic-Brasov motorway, and consists of geotechnical drillings, vertical electrical surveys (VES), geological mapping and geotechnical laboratory analysis, together with a detailed geological mapping completed in the autumn of 2017. In the Remote Sens. 2021, 13, 2385 6 of 22 2014 geotechnical study, drilling was performed with wireline systems, with freshwater fluid and bio-degradable additives. The drilling diameter is 146 mm with a 100 mm core diameter, and both disturbed and undisturbed samples were collected. In the case of disturbed samples, they are collected directly from the core box. The undisturbed samples ore obtained with Shelby tube from the drilling hole. The undisturbed samples are 87 mm in diameter and are 500 mm long. The following laboratory tests are performed by a private company in a grade I (ISO9001 9001) geotechnical laboratory: • Granulometry analysis; • Natural moisture analysis (w); • Plasticity analysis-Atterberg limits (liquid limit W L , plastic limit W P , consistency index I C , plasticity index I P ); • Structure index analysis (moist volumetric weight γ, dry volumetric weight γ d , voids index e, porosity n); • Monoaxial compression resistance (compression resistance Rc).
Although these laboratory parameters might represent point singularities and unique characteristics, they can be extrapolated as a general description of each lithological unit identified in the study area (Table 1).
Another important set of data are obtained from a detailed vertical electrical sounding (VES) performed by a private company (Search Corporation SRL, Bucharest, Romaniahttp://www.searchltd.ro/, accessed on 1 May 2021). VES provide the electrical resistivity (conductivity) of rocks using the flow of an induced electric current inside the earth. VES is performed at 46 stations along the Prahova Valley, and provides valuable information regarding the resistivity distribution, which can be interpreted in terms of depth-location and thickness of several main geological layers ( Figure 2). Therefore, we obtain a detailed spatial image of various subsurface lithologic units, and also estimate the depths of bedrock surfaces.
Our geological mapping is performed using data from walk-over surveys, geological site investigations in outcrops and anthropic openings and cuttings, excavations, and several existing water wells. Finally, we merge all the results with those from 1/200,000 and 1/50,000 geological maps and propose the stratigraphic succession shown in Figure 2 and Table 1.
Field-mapping, with stereoscopic interpretation of aerial photographs, reveals the location of unhooks and normal and reverse faults. We obtain a detailed complex tectonic pattern consisting of 29 well identified unhooks, 342 normal faults and 144 reverse faults ( Figure 2).

InSAR Data
The Prahova study area is covered by Sentinel-1 images acquired in two frames, on descending and ascending orbits ( Table 2). In the study area, slow-developing landslides are expected, therefore stacks of images with a temporal resolution of 1 month were considered in order to reduce the computing time. Two stacks of 123 images acquired between October 2014 and October 2018, out of which are 62 on ascending orbit and 61 are on descending orbit, were processed separately in order to derive surface motion from different geometries. Table 1. Litho-stratigraphic classes and geotechnical parameters used in the study area.

Stage/Age Series/Epoch Classes Cohesion [Pa] Friction Angle [ • ] Unit Weight [N/m 3 ]
Undifferentiated on the geologic map in Figure 2 Pleistocene  Figure 3 presents the extent of the area covered by the images acquired on each orbit. The temporal resolution of the Sentinel-1A satellites was initially of 12 days, but increased later to 6 days after the launch of Sentinel-1B satellite. The extent of the area selected for InSAR processing is conditioned by the limitations of the remote sensing technique and by the results obtained from the deterministic approach. Thus, the Carpathian area of the Prahova Valley is eliminated mainly due to a lack of coherent targets for the radar signal. In addition, a large difference between the Doppler centroids of the SAR images acquired in ascending orbit limit the common coverage area by all the images in the stack to the Subcarpathian reach of the valley.

The Infinite Slope Model
In order to compute the safety factor (SF) of slopes and to derive the landslide susceptibility map for a selected study area susceptible to translational landslides, a onedimensional deterministic stability model (infinite slope model) is applied in a raster-based GIS environment (ILWIS) [22]. Figure 4 presents an outlook on the way in which slope specific stability is approached. The workflow is separated into three phases: the preparation (preliminary) phase, the analytical phase and the validation. The infinite slope model is seated on the limit-equilibrium approach and on the principle that failure occurs when disturbing forces exceed the capacity to resist of the earth materials [69]. The process is enhanced by the effects of water. A detailed description of the geotechnical, hydraulic, and geometric principles for the SF calculation is provided by Mergili et al. [70]. In this study, we calculate the SF using the same expression in Equation (1) as in [35], written in effective stresses without pore water pressure [22,71]: The infinite slope model is seated on the limit-equilibrium approach and on the principle that failure occurs when disturbing forces exceed the capacity to resist of the earth materials [69]. The process is enhanced by the effects of water. A detailed description of the geotechnical, hydraulic, and geometric principles for the SF calculation is provided by Mergili et al. [70]. In this study, we calculate the SF using the same expression in Equation (1) as in [35], written in effective stresses without pore water pressure [22,71]: where: c' = effective cohesion (Pa = N/m 2 ), γ = unit weight of soil (N/m 3 ), z = depth of failure surface, m = the ratio of the groundwater level to the soil depth (dimensionless), γw = unit weight of water (N/m 3 ), β = slope gradient ( • ) and φ = effective angle of shearing resistance ( • ). Every parameter (except "m") in Equation (1) represents a space variable, and we infer the SF for each pixel. The saturation depth "m" is the only time-dependent variable of the FS equation [35]. This ratio ranges from 0 to 1; 1 when the soil is completely water saturated, and 0 when the soil is dry. In the present study, we solely tested the worst-case scenario of fully saturated soil, as partial saturation is more challenging to approach geotechnically.
Despite the fact the study area is within the Vrancea seismic region, which is Europe's most seismically active zone [72,73], in our model no external forces, such as seismic loading, are taken into account. Our approximation is valid because active, shallow landslides along the Prahova Valley are rainfall-triggered events, as identified in previous studies [3,7,8,35,[66][67][68].
The landslide susceptibility maps are classified into three stability classes. For SF less than 1, the slope is in an unstable condition, and when SF is over 1.5, the slope is in a stable condition. For intermediate values of 1 or less then 1.5, the slope is considered to be in a meta-stable condition, and minor undermining factors can easily lead to instability [35].

InSAR Based Method
The SBAS technique [47] consists of an approach that reduces temporal and spatial decorrelation by deriving small baseline interferograms. Instead of using discrete coherent targets such as in the case of classic PS methods [41,42], SBAS exploits coherence over larger areas using Delaunay triangulation [74]. Because SBAS exploits coherent areas instead of discrete coherent targets, the technique is more useful for rural and natural areas.
The SBAS algorithm is applied using the SARscape software, provided by SARMAP. The steps of applying the SBAS algorithm involve pairing images for the generation of interferograms, and deriving residual height using the NASA Shuttle Radar Topographic Mission (SRTM) DEM with 1 Arc-Second Global, displacement velocity and displacement time series after the removal of atmospheric effects. Paired images are called masters and slaves. During the interferometric process, the images are co-registered in the geometry of the masters, taking into account temporal and spatial baseline-imposed restrictions. For the ascending stack, from a number of 62 images, we generate 807 interferograms that respect the constraints of a maximum temporal baseline of 365 days and a spatial baseline of 1200 m. After imposing the same constraints on the 61 images from the descending stacks, a number of 865 interferograms are also obtained. In the next step, the interferograms are filtered in order to remove the atmospheric, topographic and orbital influences on the phase variations. After that, each interferogram is inspected manually. A number of 3 low quality images from the ascending stack and 1 image from the descending stack and their corresponding interferograms are removed from the processing chain because of being affected by too large atmospheric artefacts and orbit errors that could not be corrected.
The final displacement time series for different coherent targets in the study area are obtained after applying a second atmospheric filtering for generating the final and cleaned displacement time series.

Deterministic Approach
In this study, the one-dimension infinite slope stability model is employed to assess a shallow landslide susceptibility map across a mountainous and hilly region of Romania. For determining the safety level for valley slopes, an overall factor of safety of 1.0 is applied as a criterion for stability. Using slope angles, individual maps for the sine, cosine and the square of the cosine are calculated. Geotechnical data of slope materials are obtained from geological surveys, VES, boreholes, laboratory tests, and also from the literature. In our slope stability model, we use a series of thematic maps computed specifically for slope materials: soil thickness, cohesion, unit weight, and friction angle, and the slope failure susceptibility map using the deterministic approach under completely saturated conditions ( Figure 5).
Soil thickness refers to the thickness of the overburden and is computed from the depth to bedrock (Figure 5a). We consider the depth of the contact between the overburden and the bedrock as the slip plane of failure along a slope. The contact is established according to several data sources such as drillings, VES and specific literature e.g., [8,35,[66][67][68]. Another source of data we use is the landslide inventory in which depths of the slip plane are reported. In the Subcarpathians, depths of the slip plane are between of −1 and −4 m, and in the Carpathians, a shallow depth of contact is measured, between −0.1 and −1 m. The difference is due to the geological particularities of the two areas. Hills are in general characterized by soft rocks with alternations of permeable and impermeable layers, while mountain slopes are covered by thick periglaciar debris and diluvial soil. Therefore, in the hilly area, there is a higher variation in the depth of the contact. In the mountains, the minimum soil depth is located on the Bucegi Mountains' steep slope towards the Prahova river.
If we compare the Carpathian and the Subcarpathian areas in terms of geotechnical indices that directly influence landslides, we observe that in the Carpathians, where hard rocks (schist, sandstone, hard marls and conglomerates) are visible, the cohesion and the unit weight values are high, but the friction angles values are low. Hence, this area is prone only to rock sliding. On the other hand, in the Subcarpathians, characterized primarily by predominantly soft rocks (clays, marls, weathered conglomerate and sandstone, sands and silts), the cohesion and the unit weight values are lower, and the friction angle values are higher, favoring rainfall-induced landslides. When soil depths increase, slope angle, internal friction angle and wetting conditions gain importance [35,75,76]. We find that soil cohesion is more variable than the friction angle, as reported in previous studies [77]. Figure 5e displays the results of the landslide susceptibility under completely saturated conditions. The spatial distribution map reveals that this method is relevant only for hilly regions in our study area, which have more homogenous lithology, alternating permeable/impermeable soft rocks and larger deforested areas [8]. The Prahova river has a N-S oriented course along an important fault that cuts transversal the Miocene-Pliocene anticline and syncline system [35] (Figure 2). The slopes of the upper terraces are transformed into an erosion glacis, continued with a large colluvial-proluvial glacis that is developed on the terrace II tread. Terrace II is the largest and best-preserved fluvial terrace in the Subcarpathian sector. The slopes of terrace II are affected by active landslides that heavily erode the lowest terrace tread, which appears to act as a set of stability shoulders between the sliding areas [8]. This geomorphological pattern is well represented on the landslide susceptibility map, which classifies these slopes as unstable or in a critical state. In this study, we show that this pattern is confirmed also by the InSAR results.  If we compare the Carpathian and the Subcarpathian areas in terms of geotechnical indices that directly influence landslides, we observe that in the Carpathians, where hard rocks (schist, sandstone, hard marls and conglomerates) are visible, the cohesion and the unit weight values are high, but the friction angles values are low. Hence, this area is prone only to rock sliding. On the other hand, in the Subcarpathians, characterized primarily by predominantly soft rocks (clays, marls, weathered conglomerate and sandstone, sands and silts), the cohesion and the unit weight values are lower, and the friction angle values are higher, favoring rainfall-induced landslides. When soil depths increase, slope angle, internal friction angle and wetting conditions gain importance [35,75,76]. We find that soil cohesion is more variable than the friction angle, as reported in previous studies [77]. Figure 5e displays the results of the landslide susceptibility under completely saturated conditions. The spatial distribution map reveals that this method is relevant only for hilly regions in our study area, which have more homogenous lithology, alternating permeable/impermeable soft rocks and larger deforested areas [8]. The Prahova river has a N-S oriented course along an important fault that cuts transversal the Miocene-Pliocene anticline and syncline system [35] (Figure 2). The slopes of the upper terraces are transformed into an erosion glacis, continued with a large colluvial-proluvial glacis that is developed on the terrace II tread. Terrace II is the largest and best-preserved fluvial terrace

InSAR Results
After processing the datasets acquired from ascending and descending orbits for the Prahova Valley, two deformation maps are achieved. An average number of 430,000 scatterers are obtained from each of the datasets within an area of 120 km 2 and a spatial density of 850 point/km 2 . The realized maps show an average velocity of -0.5 mm/yr with a standard deviation of ± 0.4 mm/yr in the satellite's line of sight. According to the literature, the stability threshold depends on the radar sensor used and the standard deviation of the PS dataset obtained, which is around ± 0.5 mm/yr. The most commonly defined value interval defined as stable for Sentinel-1-derived deformation maps is either between ± 1.5 mm/yr or ± 2 mm/yr [59,78]. In Figure 6, stable areas with a low movement amplitude (± 2 mm/yr) corresponding to stability are depicted in green, while areas affected by instability are emphasized in red if the target moved away from the sensor and blue if the target moves towards it. Both deformation maps show consistent movement patterns for areas with coherence values of at least 0.70, which indicates the high reliability of the results [42,47]. Some small artefacts can be observed in the deformation maps for the ascending and descending orbits ( Figure 6), with differences occurring in the type of pattern shown at the same location (either movement away or towards the sensor). There are unstable areas that are clearly depicted in both deformation maps and areas with a limited degree of correspondence from one map to another. The hilly area in the north of the valley shows an inversion of the positive and negative movement patterns of points found at the same location in both maps. This is caused by the dependency between the satellite's line of sight geometry in relation to the terrain slope, which influences the direction of the displacement depicted by the satellite. From a descending orbit, the sensor perceives the targets on the left slope as approaching it, while the targets on the right slope are distancing from the sensor. The influence of the satellite's view geometry on movement patterns derived for hilly surfaces have been discussed also by other authors in the past, who recommended using both orbits and decomposing the results for a better depiction of the actual direction of the observed displacement [50]. However, the fact that negative and positive movement trends are found at the same location in both deformation maps is an indication that these areas show consistent movement patterns.
Because deformation maps ( Figure 6) appear initially to exhibit sub-vertical displacements instead of translational movement as would be expected in the case of shallow landslides, the resulted displacement vectors obtained from ascending and descending orbits are decomposed on horizontal and vertical directions using points located in corresponding areas of 100 m 2 from both images and the known viewing angle for each location ( Figure 7A,B). From the representations of horizontal versus vertical linear velocities, we conclude that the predominant displacement in the study area is horizontal, surpassing the vertical velocity in areas affected by instability. In this case, the horizontal movement can be associated with a displacement down the slope, implying the existence of a translational displacement, characteristic to shallow landslides. Remote Sens. 2021, 13 Figure 7A, negative values depict movement towards the west, while positive values show movement towards the east. In Figure 7B, negative values correspond to vertical subsidence while positive values represent vertical uplift.

The Prediction-Rate Curve
To confirm the validity of results, we use the prediction-rate curve [79]. The prediction-rate curve involves computing the conformity between the landslide susceptibility results and the empirical observed shallow landslides, as the inventoried landslide distribution in each susceptibility class [79][80][81].
The multi-temporal landslide inventory map we employ to validate the deterministic approach is obtained from a combination of field surveys and aerial imagery. For model evaluation, we focus on shallow active and latent landslides, considering an inventory of over 400 sites (Figure 1). In our validation approach, due to the mobility of shallow landslides especially during rainy periods, we use the entire landslide affected area.
The prediction-rate curve results are shown in Figure 8, using cumulative statistics.
The prediction-rate curve shows that almost all landslides along the Subcarpathian terrace slopes are located within the most hazardous 18% area of the susceptibility map, where the InSAR processing was also employed.  Figure 7A, negative values depict movement towards the west, while positive values show movement towards the east. In Figure 7B, negative values correspond to vertical subsidence while positive values represent vertical uplift.

The Prediction-Rate Curve
To confirm the validity of results, we use the prediction-rate curve [79]. The predictionrate curve involves computing the conformity between the landslide susceptibility results and the empirical observed shallow landslides, as the inventoried landslide distribution in each susceptibility class [79][80][81].
The multi-temporal landslide inventory map we employ to validate the deterministic approach is obtained from a combination of field surveys and aerial imagery. For model evaluation, we focus on shallow active and latent landslides, considering an inventory of over 400 sites (Figure 1). In our validation approach, due to the mobility of shallow landslides especially during rainy periods, we use the entire landslide affected area.
The prediction-rate curve results are shown in Figure 8, using cumulative statistics.
The prediction-rate curve shows that almost all landslides along the Subcarpathian terrace slopes are located within the most hazardous 18% area of the susceptibility map, where the InSAR processing was also employed.
The comparison of mapped landslides and modelling results reclassified into two different classes (unstable/stable) is presented in a contingency table (Table 3). True positive refers to landslides located in areas classified as unstable and critical, and false positives outline areas predicted as susceptible but without landslides. True negative is when areas predicted as safe have no landslides, and false negative delineate the case when landslides occur in areas predicted as safe. Several quality measures are provided by linking together correct and incorrect classified positives (i.e., unstable areas) and negatives (i.e., stable areas) ( Table 3). The comparison of mapped landslides and modelling results reclassified into two different classes (unstable/stable) is presented in a contingency table (Table 3). True positive refers to landslides located in areas classified as unstable and critical, and false positives outline areas predicted as susceptible but without landslides. True negative is when areas predicted as safe have no landslides, and false negative delineate the case when landslides occur in areas predicted as safe. Several quality measures are provided by linking together correct and incorrect classified positives (i.e., unstable areas) and negatives (i.e., stable areas) ( Table 3). Table 3. Prediction errors (% of correctly and incorrectly predicted pixels) and accuracy statistics.

Odds ratio (a + d)/(b + c)
Ratio between correctly and incorrectly classified pixels 42.47

Positive predictive power a/(a + b)
The proportion of true positives in the total of positive predictions 0.7

Negative predictive power d/(c + d)
The proportion of true negatives in the total of negative predictions 0.99 Although the deterministic model is data demanding and requires a high degree of simplification of input values, the accuracy statistics derived from the contingency matrix show that the model is able to correctly classify especially safety areas (93% being correctly classified) and non-observations of the sliding phenomenon (0.3%). The result can be also related to the fact that the susceptibility map is a prediction for the hazard of future landslides and not the hazard itself.

Odds ratio (a + d)/(b + c)
Ratio between correctly and incorrectly classified pixels 42.47

Positive predictive power a/(a + b)
The proportion of true positives in the total of positive predictions 0.7

Negative predictive power d/(c + d)
The proportion of true negatives in the total of negative predictions 0.99 Although the deterministic model is data demanding and requires a high degree of simplification of input values, the accuracy statistics derived from the contingency matrix show that the model is able to correctly classify especially safety areas (93% being correctly classified) and non-observations of the sliding phenomenon (0.3%). The result can be also related to the fact that the susceptibility map is a prediction for the hazard of future landslides and not the hazard itself.
Considering the area in Figure 7 as stable should be subjected to future scrutiny, as there might be errors generated by the limitations of the method we used. The large apparent stable area in the reach of the Prahova Valley is the result of old deep-seated landslides, downfalls and debris, reinforced by vegetation, providing a non-realistic slope stability. The future motorway plans along the Prahova river corridor include massive deforestation on the left side of the range, which can result in a plethora of processes with the potential to create a catastrophic impact on the environment and the human habitat.

Comparison between Susceptibility Map and InSAR Detected Landslides
In order to assess the model performance in terms of quality and predictive power of the deterministic approach, we implement a GIS-based analysis (Figure 9). We test the level of correspondence between unstable areas detected by the InSAR technique and the landslide susceptibility map.

Discussions and Conclusions
In this study, we combine classic and modern methodologies and their cross-validation in order to obtain a reliable landslide susceptibility map. The relative probability for spatial occurrence of mass movements is derived from the computation of the SF based on the geotechnical properties of slope materials. The most important geotechnical predictors for slope failure are cohesion, unit weight and friction angle. The trigger specific for the Subcarpathian study area is the rise in groundwater table, as a result of water infiltration during heavy rainfall episodes. Therefore, our landslide susceptibility map is InSAR technology depends on radar coherent back scatterers, which are mainly found in built-up environments. However, in our case study, the landslide-susceptible areas are located on slopes, with a limited number of buildings and infrastructure. Additionally, active shallow landslides are characterized by the presence of vegetation, which is a main cause for signal decorrelation. Therefore, it was taken into account that the obtained InSAR deformation maps depict instability mainly in built-up areas, where the phenomenon is related to subsidence and uplift processes, which are in most cases caused by human impact, such as groundwater exploitation and construction works.
Points found on horizontal surfaces displaying velocities higher than +2 mm/yr in the sensor LOS (line of sight) were assumed to correspond to subsidence or uplift on a vertical direction and attributed to other tectonically driven processes. In order to eliminate the instability that is unrelated to landslides, the InSAR points found on horizontal and quasi-horizontal surfaces are discarded, starting from the point that landslides develop on slopes. For this purpose, all InSAR points located on areas with a slope gradient lower than 3% were eliminated from the analysis. The calculation of the slope gradient in the study area was based on a 1 m spatial resolution Lidar DEM resampled to 10 m resolution to match the spatial distribution of the InSAR points. A geomorphological map was applied to discard all InSAR points located along the river meadow and terrace treats. This procedure drastically reduced the number of available InSAR points derived from each orbit, decreasing the density of InSAR points from approximately 850 points/km 2 to 120 points/km 2 in the study area.
The unstable InSAR points from images obtained from both orbits are selected by imposing a threshold of stability between −2 mm/yr and 2 mm/yr, respectively, considering the InSAR results accuracy of ± 0.4 mm/yr and the assumption that movement lower than 2 mm/yr can be associated with stability. Points that display velocities outside the stability interval are associated with potential slope failure. By applying this binary classification, two InSAR deformation maps at 20 m resolution with stable/unstable condition classes are generated. The stability condition classes derived from each orbit were then merged into a single map containing InSAR-derived stable/unstable conditions. The susceptibility map is also classified into stable/unstable condition classes by merging the unstable and quasi-unstable classes into one singular class and resampling the data to a 20 m resolution raster. The two binary maps are then spatially cross correlated to derive the accuracy of the susceptibility map. In order to predict the ongoing landslides as indicated by In-SAR and assess the correspondence between all defined classes, we present a confusion matrix in Table 4. The comparison between InSAR and the susceptibility map reveals that the susceptibility analysis is able to predict more than 22% of the active landslides identified by InSAR (Table 4). This result is also conditioned by the constraint of back scatterers in coherent areas. In our view, most of the 78% missed ongoing displacements assemble all kind of ground movements, not only landslides (for example human impacts such as ground water extraction). As expected, the highest accuracy is obtained for the estimation of stable areas, which the deterministic map detected with an accuracy of around 82%.

Discussions and Conclusions
In this study, we combine classic and modern methodologies and their cross-validation in order to obtain a reliable landslide susceptibility map. The relative probability for spatial occurrence of mass movements is derived from the computation of the SF based on the geotechnical properties of slope materials. The most important geotechnical predictors for slope failure are cohesion, unit weight and friction angle. The trigger specific for the Subcarpathian study area is the rise in groundwater table, as a result of water infiltration during heavy rainfall episodes. Therefore, our landslide susceptibility map is valid under slope saturation conditions. In the Carpathian area, the most important characteristic of the lithology is the hardness of the rock, hence the values of cohesion and unit weight are high, and the values of friction angles are low. For this reason, even if we increase wetting conditions, the overall effect in terms of slope failure is negligible. Moreover, the mountain area of the Prahova Valley, characterized by flysch formations, is still well-forested, and more susceptible to deep-seated landslides. Old landslides located under forest are present, and this might affect our ability to capture all dynamic patterns of slope processes.
On the other hand, in the Subcarpathians, the rocks are soft, with alternating permeable/impermeable layers that trigger slope failure both in saturated and unsaturated conditions. For these soft rocks, cohesion and unit weight values are low, but friction angle values are higher compared to the mountain area. Increases in the water content of soil also increase the friction angle, thus the saturated condition has a much higher influence in slope failure. Since the Subcarpathian Prahova Valley is much more inhabited than the Carpathian area, an additional underground water source from water sewage spill and wastewater might cause higher landslide susceptibility classes.
The study reveals that in the Subcarpathians, the Breaza syncline, which is filled with a succession of Miocene clays, marls and sandstones, represents the most susceptible area and records the highest density of shallow-to-medium deep landslides affecting the overburden and a small part of the bedrock [8,35]. Some variations in type, size and orientation of landslides are induced by local strata gradients and the presence of faults. The pattern for future landslides is controlled by two morphological alignments: the continuous slope of the second terrace and the glacis of the upper third and fourth terraces.
In the hilly areas, landslides develop in relation to the length of the slope. On the second terrace slope, larger landslides are embedded on the southern flank of the syncline, forming amphitheater-shaped landslide heads. On the northern wing of the syncline, where Brebu conglomerates outcrop, there are shallower, stabilized and drained landslides. Overland flows become prevalent here, because of impermeable rocks [8].
Results are validated using a statistical approach and InSAR. Although the infinite slope stability model applied in this study is influenced by uncertainties in measured geotechnical characteristics, the validation assessment and InSAR approach revealed a good performance in reproducing the observed landslide areas in the Subcarpathian reach of the Prahova Valley. In this hilly region, landslides are found in unstable and critical predicted areas, and the statistical validation of modelling results shows a good predictive capacity. For those areas with higher coherence for the radar signal, the SBAS-InSAR technique provided a good confirmation of slope dynamics, proving to be an adequate tool for landslide analysis.
The almost 23% of spatial overlay between the deterministic susceptibility map and the SBAS-InSAR ground deformation map is a result heavily conditioned by the presence of back scatterers. This result represents mostly areas characterized by regressive-developing scarps in the built-up areas or captures the presence of man-made features located on unstable slopes. The result of the validation for a potentially extreme scenario through InSAR is suited for revealing the slow ongoing movements in their initial developing state, with all the limitations imposed by the presence of signal coherent areas. Our SF map represents the probability of the occurrences for future landslides. The absence of movements detectable by InSAR for over 80% of the stable areas indicated as landslidesusceptible does not necessarily imply that these areas will not be affected by landslides in the future.
In order to improve our classic deterministic analysis, the susceptibility map is updated by integrating it with InSAR-detected instabilities. The areas defined as stable, quasi-unstable or unstable in the deterministic map are updated to a confirmed instability class if InSAR technique detects active movement ( Figure 10).
Since built-up areas show the best candidates for radar scatterers, the proposed methodology brings the most reliable results for landslide detection in populated environments, while capabilities of the current satellite SAR sensors are still limited in natural, highly vegetated areas. Several authors have proposed improved algorithms that increase coherence in natural areas (SqueeSAR) [82][83][84], but in the current study area, the presence of high vegetation causes a temporal decorrelation of the C-band signal regardless of the selected method. However, we show that the use of InSAR for monitoring slow-developing landslides affecting inhabited and constructed areas is meaningful and produces satisfactory results. Considering that the study area presents good reflectivity and coherence, the proposed technique offers a highly accurate estimation of ground displacement on sliding slopes.
The characteristics of each method makes them more suitable for different types of environments. While InSAR is able to confirm active movements in built-up environments, the deterministic method focuses on vegetated slopes that limit the applicability of the remote sensing technique. Therefore, refined maps of landslide susceptibility obtained by integrating these two complementary methods offer a more complete view of an area of interest, and enhance landslide susceptibility estimations in areas with existing infrastructure and buildings, or where future key infrastructure projects are planned.

Data Availability Statement:
The data generated and studied in this article can be obtained upon request to the principal author.