Assessing Urban Landslide Dynamics through Multi-Temporal InSAR Techniques and Slope Numerical Modeling

: Landslides threaten more than before the urbanized areas and are a worldwide growing problem for the already affected communities and the local authorities committed to landslide risk management and mitigation. For this reason, it is essential to analyze landslide dynamics and environmental conditioning factors. Various techniques and instruments exist for landslide investigation and monitoring. Out of these, Multi-temporal Synthetic Aperture Radar Interferometry (MT-InSAR) techniques have been widely used in the last decades. Their capabilities are enhanced by the availability of the active Sentinel-1 mission, whose 6-day revisiting time enables near real-time monitoring of landslides. Interferometric results, coupled with ground measurements or other approaches such as numerical modeling, signiﬁcantly improve the knowledge of the investigated surface processes. In this work, we processed the C-band SAR images of the available European Space Agency (ESA) satellite missions, using MT-InSAR methods to identify the surface deformations related to landslides affecting the Ias , i Municipality (Eastern Romania). The results (i.e., velocity maps) point out the most active landslides with velocities of up to 20 mm/year measured along the satellite Line of Sight (LOS). Following, we focused on the most problematic landslide that affects the T , ic ă u neighborhood and is well-known for its signiﬁcant implications that it had. To better understand its behavior and the sensitivity of the displacements to the environmental factors (i.e., rainfall), we carried out 2D numerical modeling using a ﬁnite difference code. The simulated displacement ﬁeld is consistent with the InSAR displacements and reveals the most active sectors of the landslide and insights about its mechanism.


Introduction
The human impact, through topographic and hydrologic modifications, is a triggering and preparatory factor for landslides. At the same time, as landslides evolve, they will control human activities, and a feedback loop will occur. Landslides limit and interfere with the infrastructure development, utility networks or may obstruct the expansion of urbanized areas. In the worst-case scenario, they destroy and inflict damages to buildings, infrastructure, and city heritage [1] and generate substantial economic losses [2][3][4][5][6][7]. Many sensing techniques, tools, and instruments are developing to detect, monitor, and analyze landslide processes. The Multi-temporal Interferometric Synthetic Aperture Radar (MT-InSAR) techniques are one of them, and more and more, they are used to investigate landslide hazards. The reliability and consistency of their results make them a great source of information for supporting decision-making policies for landslide risk and management.
Differential SAR Interferometry (DInSAR) techniques use radio's microwave properties to estimate the rate of ground displacement with millimeter accuracy. The measurement of relative surface displacement between two SAR acquisitions, calculated as the distance difference between a satellite antenna and ground objects, is recorded as the phase shifting of the backscattered echo [8][9][10]. The exploitation of the signal's phase property over a specific area results in the generation of an interferogram [8,11]. The interferogram represents a matrix of numerical values corresponding to phase variations, usable for the generation of Digital Elevation Models (DEMs) [12][13][14] or the detection of ground surface changes [10,15,16]. The working principle of DInSAR to quantify the displacements relies on the removal of topographic influence out of the original interferogram by using a DEM. More comprehensively, the backscattered signal is the sum of many contributors whose effects have to be considered and removed to isolate the displacement fringes. Besides the displacement component and topography, the interferometric phase is also subjected to the atmosphere's influence, the geometric and temporal decorrelations introduced by possible orbital errors, and the phase noise [11,[17][18][19][20]. In the last two decades, different sophisticated algorithms that consider more than three SAR acquisitions to calculate and remove such influences were developed. Although they have different strategies to estimate the displacement component, these multi-temporal interferometry techniques [21] such as Permanent Scatterer Interferometry (PSI) [22,23] and Small BAseline Subset (SBAS) [24] can overcome the DInSAR limitations and provide millimeter accuracy measurements over a long period. Hence, their potential to investigate landslides rapidly increased, not only for landslides mapping and monitoring their behavior and dynamics but also for detecting early failure indicators [18,21,.
The use of numerical simulations for landslide investigation aims to appropriately determine the stability conditions of the slope and its potential failure mechanisms [48]. Identifying the factors that destabilize the slope and lead to the triggering of a possible event can be used to assess and mitigate the hazard. Accordingly, effective policies and strategies may be proposed and implemented to stabilize the displaced material [49][50][51]. Numerical modeling of slopes makes use of mathematical equations to solve the mechanical response of the unstable mass. Regardless of the employed technique, i.e., Finite Elements (FE), Discrete Elements (DE), or Finite Difference Method (FDM), these approaches focus on deformation analysis as well as the safety factor of the slope [48].
Many studies integrate MT-InSAR techniques and numerical modeling due to their high effectiveness in monitoring and investigating the state of landslides dynamics [46,[52][53][54]. Using both approaches is highly beneficial as they allow us to evaluate and monitor the landslide activity over large areas and, additionally, to perform fast investigations of specific cases where fast-paced evolving surface deformations are detected. The outcome of MT-InSAR techniques consists of deformation velocity maps and displacement time series. Each measurement quantifies, with millimeter accuracy, the change of the distance between sensor and target in the satellite Line Of Sight (LOS). Thus, the MT-InSAR outputs provide information about the surface velocity and the extent of the process. At the same time, the simulations can evaluate multiple failure scenarios to explain the material's physical behavior of the observed displacements. Their capabilities are proven to work not only in the study of landslides but also in other fields such as dam monitoring [55] and mining activities [56,57]. Even though these two complementary techniques are employed more and more in the study of earth surface deformations, they are primarily used in areas supported by lots of in situ data [53,55,58,59]. In contrast, their use in regions lacking ground measurements is still challenging and not yet fully exploited.
This study firstly uses the MT-InSAR techniques to identify the critical areas within Ias , i Municipality showing landslide-related deformations. Further on, we focus on the T , icău landslide, one of the active landslides that affect the residential neighborhoods. The challenging part of the investigation is related to the lack of in-situ data required to model Remote Sens. 2021, 13, 3862 3 of 26 the behavior of landslide material. To understand its failure mechanism and dynamics, we set up a two-dimensional Finite Difference section along the slope and run several scenarios which potentially describe the observed deformations and pattern. Lastly, we analyzed the simulations and the MT-InSAR results, the velocity maps, the displacement time series, and the magnitude of displacements to point out the landslide body's active and most dangerous parts and support the proposed sliding mechanism and the geomorphological landslide type.

Study Area
The Municipality of Ias , i (91.5 km 2 ) is the most important city of the north-eastern part of Romania ( Figure 1). This area is constantly facing landslide-related deformations being not only one of the landslide hot-spots of the nation [60], but it is recognized at the European scale as well [61]. In the northern part of the city, the terrain morphology is characterized by gentle hillslopes and cuesta landforms with elevations of 200-220 m asl. In contrast, a hilly morphology defines the southern side with steep escarpments and higher altitudes of up to 404 m asl. scenarios which potentially describe the observed deformations and pattern. Lastly, we analyzed the simulations and the MT-InSAR results, the velocity maps, the displacement time series, and the magnitude of displacements to point out the landslide body's active and most dangerous parts and support the proposed sliding mechanism and the geomorphological landslide type.

Study Area
The Municipality of Iași (91.5 km 2 ) is the most important city of the north-eastern part of Romania ( Figure 1). This area is constantly facing landslide-related deformations being not only one of the landslide hot-spots of the nation [60], but it is recognized at the European scale as well [61]. In the northern part of the city, the terrain morphology is characterized by gentle hillslopes and cuesta landforms with elevations of 200-220 m asl. In contrast, a hilly morphology defines the southern side with steep escarpments and higher altitudes of up to 404 m asl. Figure 1. Location of Iași Municipality (Northeastern Romania) and its terrain morphology. The landslide inventory is from [62], and the numbers in parentheses represent the landslide type based on the classification of [63]. (high-resolution image is available at https://doi.org/10.6084/m9.figshare.14915103.v2).
The morphology of the Iași area developed on a monoclinic geological structure [64][65][66] under the migration and incision of Bahlui River and its tributaries towards the south, Figure 1. Location of Ias , i Municipality (Northeastern Romania) and its terrain morphology. The landslide inventory is from [62], and the numbers in parentheses represent the landslide type based on the classification of [63]. ( The morphology of the Ias , i area developed on a monoclinic geological structure [64][65][66] under the migration and incision of Bahlui River and its tributaries towards the south, which generated asymmetric hills and structural valleys [67][68][69]. The present-day morphostructure, particularly the lithology, coupled with the small depths of the phreatic level and springs [70], serves as key-conditioning factors for landslide reactivations. The lithology consists of claystones and mudstones interbedded with sandy layers deposited during the Middle Miocene [65,71,72]. In the southernmost part of the area, where the elevations are higher than 320 m asl, a caprock of calcarenites, oolitic calcarenites, and quartz arenites interbedded by sands and claystones outcrops [65,[71][72][73]. Overlaying these deposits is a 40 m thick layer of sand. The above described geological formations are covered on ridges and gentle hillslopes by quaternary deposits consisting of fluvial gravels and clays with varying thicknesses from 2 to 5 m, and loess that measures up to 25 m in some cases [70,74,75]. The dry climate is specific for the area and has mean annual rainfall values of 560 mm/year [76]. Rainy periods were recorded ( Figure 2) in 1920-1950 and 1960-2000, a situation which might indicate a 30-40 years humid cycle separated by shorter dry periods [77]. These rainy periods caused the increase of the groundwater table, favoring the occurrence of landslides [77][78][79].
The slopes are shaped by landslides [80][81][82] which might reactivate after rainfall events, particularly during the spring season [62,77]. Many cases of landslide reactivations were reported by various authors [62,78,79,83,84] or by local authorities and mass media, as they repeatedly damaged the residential buildings, road network, and utility infrastructure.
which generated asymmetric hills and structural valleys [67][68][69]. The present-day morphostructure, particularly the lithology, coupled with the small depths of the phreatic level and springs [70], serves as key-conditioning factors for landslide reactivations. The lithology consists of claystones and mudstones interbedded with sandy layers deposited during the Middle Miocene [65,71,72]. In the southernmost part of the area, where the elevations are higher than 320 m asl, a caprock of calcarenites, oolitic calcarenites, and quartz arenites interbedded by sands and claystones outcrops [65,[71][72][73]. Overlaying these deposits is a 40 m thick layer of sand. The above described geological formations are covered on ridges and gentle hillslopes by quaternary deposits consisting of fluvial gravels and clays with varying thicknesses from 2 to 5 m, and loess that measures up to 25 m in some cases [70,74,75].
The dry climate is specific for the area and has mean annual rainfall values of 560 mm/year [76]. Rainy periods were recorded ( Figure 2) in 1920-1950 and 1960-2000, a situation which might indicate a 30-40 years humid cycle separated by shorter dry periods [77]. These rainy periods caused the increase of the groundwater table, favoring the occurrence of landslides [77][78][79].
The slopes are shaped by landslides [80][81][82] which might reactivate after rainfall events, particularly during the spring season [62,77]. Many cases of landslide reactivations were reported by various authors [62,78,79,83,84] or by local authorities and mass media, as they repeatedly damaged the residential buildings, road network, and utility infrastructure.
The slope investigated in detail ( Figure 3) was affected by a relevant landslide reactivation in the spring of 1942. More than 350 houses were destroyed by the event, along with the road network. The event created the main scarp of 6-7 m in height and about 800 m in length [78]. Downslope, the material reached the Cacaina River floodplain and generated typical morphology with mounds and micro depressions of 1-2 m of height difference, out of which spring waters emerged more than one month after the event. The triggering factor was the prolonged rainfall, which started in the previous autumn, coupled with the snowmelt.  [85], and 1961-2018 data from ECA&D [86,87]). The red line indicates the polynomial loess smooth trend with a 0.5 span. At the bottom are represented the landslide events in the Municipality of Iaşi as reported by [77,79].   [85], and 1961-2018 data from ECA&D [86,87]). The red line indicates the polynomial loess smooth trend with a 0.5 span. At the bottom are represented the landslide events in the Municipality of Iaşi as reported by [77,79]. (high-resolution image available at https://doi.org/10.6084/m9.figshare.14915148.v1).
The slope investigated in detail ( Figure 3) was affected by a relevant landslide reactivation in the spring of 1942. More than 350 houses were destroyed by the event, along with the road network. The event created the main scarp of 6-7 m in height and about 800 m in length [78]. Downslope, the material reached the Cacaina River floodplain and generated typical morphology with mounds and micro depressions of 1-2 m of height difference, out of which spring waters emerged more than one month after the event. The triggering factor was the prolonged rainfall, which started in the previous autumn, coupled with the snowmelt.
After the event, several remedial works, such as drainage systems and passive mechanic infrastructures, were carried out to drain the area and stabilize the sliding material. However, the passive stabilization infrastructure and drainage systems designed to prevent further displacements do not work anymore, probably due to degradation, and slope deformations are currently active. After the event, several remedial works, such as drainage systems and passive mechanic infrastructures, were carried out to drain the area and stabilize the sliding material. However, the passive stabilization infrastructure and drainage systems designed to prevent further displacements do not work anymore, probably due to degradation, and slope deformations are currently active. At the moment, the area is characterized by hummocky-like landforms and is affected by very slow-moving slides. Pieces of evidence of surface displacements were identified during in-situ surveys and manifest as wall cracks, wall bulging, cracks in the buildings' outer walls, damages to passive stabilization infrastructures, tilted trees, leaned electric poles, road cracks, and bumps ( Figure 4).   At the moment, the area is characterized by hummocky-like landforms and is affected by very slow-moving slides. Pieces of evidence of surface displacements were identified during in-situ surveys and manifest as wall cracks, wall bulging, cracks in the buildings' outer walls, damages to passive stabilization infrastructures, tilted trees, leaned electric poles, road cracks, and bumps ( Figure 4).

SAR Data and MT-InSAR Methods
The landslide deformations that affect the Municipality of Ias , i were investigated for more than 25 years, from 1992 to 2018, using C-band SAR data acquired by different sensors (Table 1). We processed the available ERS-1, ERS-2, ENVISAT, and Sentinel-1 SAR images provided by the European Space Agency (ESA). Both satellites acquisition geometries, ascending and descending, were used for our analysis to cross-validate our outputs and increase the spatial extent of information retrieval. Additionally, we used a 5 × 5 m highresolution LiDAR DSM, provided by Prut-Bârlad Water Administration, to remove the topographic phase and geocode the results. The data stacks were processed using the Interferometry Stacking tool of Sarscape v5.4. This tool includes Permanent Scatterers Interferometry (PSI) [17,22,23] and Small BAseline Subset (SBAS) [24,88] processing techniques. The use of one method or the other depends on the number of SAR images available and the ground conditions [21,[89][90][91]. These advanced techniques employ different strategies to provide high accuracy results and have their pros and cons. Their development overcomes the limitations of DInSAR, such as the temporal and spatial decorrelations, the phase unwrapping errors, and the artifacts due to the atmospheric component [18,20,21].
PSI is the first stacking interferometry algorithm developed to overcome some of DInSAR shortcomings [17,22,23]. It considers the entire dataset of SAR acquisitions to create interferograms between a so-called master image and all the other slaves with no restrictions in terms of temporal and spatial baselines (Figure 5a). After the connection graph generation, described in the previous sentence, the PS workflow implemented in Sarscape consists of the following steps: (i) interferometric process, which executes the following sequence automatically: coregistration, interferogram generation, flattening using the high-resolution DSM, and PSs selection by considering the Amplitude Dispersion Index (ADI); (ii) first inversion during which the topographic residuals and the displacement velocity are estimated by using a linear velocity model and removed from the complex interferograms; (iii) the second inversion uses the linear model products calculated previously to estimate the atmospheric components by using a low-pass filter (1000 m) and a high-pass filter (365 days), and fit the final displacement velocity model; (iv) geocoding of the final PS products. Remote Sens. 2021, 13, x FOR PEER REVIEW 8 of 26 The permanent scatterers (PS) are targets on the ground that are not affected by temporal decorrelation and maintain their signal phase-coherent throughout the analyzed period. Because of this, in Sarscape, the approach is limited to a certain number of sufficiently high coherent targets and requires a large number of images. This characteristic might be a constraint of technique for non-urban areas and areas dominated by vegetation. However, it is reliable in urbanized areas as the buildings maintain a coherent backscattering signal in time. In our case, we considered only the PS measurements with coherence greater than 0.75 and discarded those lower than this threshold. Thus, this approach is suitable for our study case, especially with Sentinel-1A/B datasets characterized by a short repeating cycle between acquisitions.
The SBAS implementation in Sarscape exploits the phase difference of spatially distributed scatterers to measure the ground deformation. Compared to the PSI approach, this technique is feasible to be used with a smaller stack of images as it pairs all the SAR images defined by user-specified temporal and spatial baselines ( Figure 5b) to reduce the decorrelations and phase noise. The multilooking of data and the small subsets connection increase the spatial information achievable by analyzing a smaller stack of SAR images. However, as an effect of multilooking, this approach is not suitable to detect local deformation. In our case, we applied the SBAS technique to ERS 1/2 and ENVISAT datasets as the number of available images was limited.
The workflow for the SBAS module in Sarscape is: (i) connection graph generation, explained earlier; (ii) interferometric step, which generates a stack of interferograms, filters them using the Goldstein adaptive filter, and then the phase is unwrapped using the Delaunay Minimum Cost Flow method; (iii) refinement and re-flattening, which estimates and removes the topographic and constant phases as well as the possible phase ramps from the unwrapped stack; (iv) first inversion, the first estimate of the displacement rate and residual heights are calculated and used to re-flatten the interferograms. A second unwrapping is also performed; (v) the second inversion removes the atmospheric phase component to clean and calculate the final displacement velocity and time series; (vi) geocoding of the final products.
Since the LOS displacements are 1D, we used the method of [29] to project the velocity along the slope for the Sentinel-1 PS. To resolve the components of the velocity, we also used [29] approach to compute the E-W and vertical displacements, using the information from both orbits also for Sentinel-1 PS (points in 10-meter vicinity were selected). The permanent scatterers (PS) are targets on the ground that are not affected by temporal decorrelation and maintain their signal phase-coherent throughout the analyzed period. Because of this, in Sarscape, the approach is limited to a certain number of sufficiently high coherent targets and requires a large number of images. This characteristic might be a constraint of technique for non-urban areas and areas dominated by vegetation. However, it is reliable in urbanized areas as the buildings maintain a coherent backscattering signal in time. In our case, we considered only the PS measurements with coherence greater than 0.75 and discarded those lower than this threshold. Thus, this approach is suitable for our study case, especially with Sentinel-1A/B datasets characterized by a short repeating cycle between acquisitions.
The SBAS implementation in Sarscape exploits the phase difference of spatially distributed scatterers to measure the ground deformation. Compared to the PSI approach, this technique is feasible to be used with a smaller stack of images as it pairs all the SAR images defined by user-specified temporal and spatial baselines (Figure 5b) to reduce the decorrelations and phase noise. The multilooking of data and the small subsets connection increase the spatial information achievable by analyzing a smaller stack of SAR images. However, as an effect of multilooking, this approach is not suitable to detect local deformation. In our case, we applied the SBAS technique to ERS 1/2 and ENVISAT datasets as the number of available images was limited.
The workflow for the SBAS module in Sarscape is: (i) connection graph generation, explained earlier; (ii) interferometric step, which generates a stack of interferograms, filters them using the Goldstein adaptive filter, and then the phase is unwrapped using the Delaunay Minimum Cost Flow method; (iii) refinement and re-flattening, which estimates and removes the topographic and constant phases as well as the possible phase ramps from the unwrapped stack; (iv) first inversion, the first estimate of the displacement rate and residual heights are calculated and used to re-flatten the interferograms. A second unwrapping is also performed; (v) the second inversion removes the atmospheric phase component to clean and calculate the final displacement velocity and time series; (vi) geocoding of the final products.
Since the LOS displacements are 1D, we used the method of [29] to project the velocity along the slope for the Sentinel-1 PS. To resolve the components of the velocity, we also used [29] approach to compute the E-W and vertical displacements, using the information from both orbits also for Sentinel-1 PS (points in 10-m vicinity were selected).

Statistical Post-Processing of PS Measurements
To objectively evaluate the areas that show a consistent deformation trend, we used the non-parametric Mann-Kendall test [92] to identify the PS points that show a displacement time series trend (which might be different from a linear trend, respectively monotonic) and the slope of the linear trend which assesses the magnitude of the deformation for automatic filtering purposes. We employed the Man-Kendall test (MK) in Rstat through the mk.test function of the trend package [93]. Each point's resulting p-value is used to select those that have a significant trend (p-value smaller than 0.0000005). Further filtering is done using the linear regression slope, which should be higher/smaller than 0.15/−0.15 (depending on the satellite's orbit) to eliminate the points that show a weaker trend. For the selection of the clusters characterized by high deformation rates, the kernel density (density function) of the PS points was estimated using velocity as a covariate (rhohat function) in the spatstat R package [94,95]. This approach can point out areas that show deformation rates over larger spatial extents and remove single or a few deformation points. That might be related to land-use change (in our study, this situation is related to new buildings constructed during the monitoring period) or particular topographic cases (new buildings constructed over disturbed material fillings that do not reflect the hillslope's state).
To investigate the relationships between displacement and rainfall, we used the breakpoint and trend calculation R package greenbrown [96,97]. In our case, the breakpoints and trends were computed based on the quantile regression to the median, using the Trend function from the greenbrown R package. For the detection of the breakpoints, the approach uses iteratively the ordinary-least squares moving sum (MOSUM) test (if the test indicates a significant structural change at p ≤ 0.05), the number and location of the breakpoints are estimated by minimizing the Bayesian Information Criterion (BIC), and by reducing the residual sum of squares of this regression respectively [96][97][98][99]. The displacement rate trends are correlated with the precipitation of ECA&D data [86,87] for the Iaşi meteorological station to potentially identify the behavior of representative targets over the T , icău landslide after rainfall events.

Ground-Based Data, Model Setup, and Numerical Modeling
The investigated area faces many shortcomings in terms of geotechnical data and insitu measurements. Currently, the landslide is not well assessed as there are no instruments installed over the site or monitoring surveys carried on. The available information consists of several exploratory borings drilled in 2006 and published in the Ph.D. thesis of [100].
Their description provides information about the lithological limits, the water table level measured during the drilling, the unit weight of materials, plasticity indices, and intrinsic strength parameters.
The slope stability analysis was done using the 2D Finite Difference code along the section crossing the landslide (Figure 2a). We subdivided the model into seven units of geo-materials ( Figure 2c): a succession of five layers in the upper part of the slope, the landslide debris unit in the middle and lower sectors, and the geological bedrock. To these units of geo-materials, we assigned the properties available in the literature.
The model consists of elements, nodes, and zones, to which specific physical and mechanical parameters are attributed and used to compute the internal stresses and strains within the slope. To perform the simulation of the behavior of the slope, the initial state of the model is subjected to various parametric changes related to its geometrical and mechanical characteristics, cyclic stresses, and potential anthropic activities or the increase in rainfall intensity which might lead to landslides triggering. The typically achieved solution concludes when the internal stresses reach the state of equilibrium, and no more strains develop, or when the collapse occurs.
Due to the many parameters required to perform numerical simulations, we completed our database by making some assumptions and empirical correlations ( Table 2). The massdensity of soil was easily derived through the formula: where ρ is mass density, γ is unit weight, and g is gravity. Further on, we assumed the tensile strength of material high enough for the upper layers and the bedrock to create a relatively stable medium, as suggested by the PSI results. For the landslide unit, which we mostly stressed in our simulations, we took the tension as one-fifth of the cohesion. Based on the literature review and the authors' knowledge, we empirically assigned the Poisson's ratio (see Table 2). Young's modulus was calculated as a relationship of oedometric modulus, available in the borings description, and Poisson's ratio: where E isYoung's modulus, M is oedometric modulus, and ν is Poisson's ratio. We employed the Finite Difference code in the FLAC (Fast Lagrangian Analysis of Continua) environment to carry out our modeling. We followed the general workflow starting with the grid generation, choice of the constitutive model and assignment of material properties, setting the boundary and initial conditions, sequential modeling, interpretation of results, and, if necessary, model refinement. We performed a static analysis using the simple elastoplastic Mohr-Coulomb model, which assumes that the maximum shear stress controls failure and that this failure shear stress depends on the normal stress. Hence, failure occurs when the applied shear stresses equal the available resistances: where τ is shear stress, c is effective cohesion, σ is the normal effective stress, and φ is the effective friction angle. Given the variability of the available indices and parameters observed in the log description of the wells, we performed various simulations with different values to reduce their uncertainty. As our most interest is related to the landslide unit, we considered a range of variables for the intrinsic parameters (Table 2) to identify corresponding changes in the dynamics and mechanics of the landslide. Specifically, we simulated the landslide mass behavior for different combinations of friction angle, cohesion, and tension. Moreover, we considered a possible worst-case scenario similar to the one that occurred during the 1942 reactivation. In this case, we raised the level of the water table close to the topographic surface.

MT-InSAR Outputs
The results of Sentinel-1 SAR data analyzed with the PSI approach consist of mean velocity maps of single ground targets and time series of displacements. The velocities obtained over the entire study area range from −24.3 to +16.9 mm/year in the ascending orbit ( Figure 6a) and from −24.2 to +12.9 mm/year in the descending one (Figure 6c). Positive values point out that the distance between target and sensor decreases over time, while the negative values suggest an increase of target-to-sensor distance as the target moves away from the satellite. The identified PSs in the descending orbit are 131,670, of which around 92.7% are stable points. In the ascending orbit were identified 125,022 PSs, of which 93.1% are stable. We consider as stable points those with a velocity between −2 and +2 mm/year in our analysis. This threshold was chosen based on the precision of MT-InSAR techniques [18,23,101] and on the specific observation on the PS results of Copou (stable plateaux) and T , icău area (instable hillslope) in [102].
The velocity of deformations recorded over the affected slopes ranges from ±7 mm/y to ±24 mm/y, measured along the satellite's LOS, depending on the orbit's geometry and the slope angle. North-eastern and eastern facing unstable slopes display positive changes in the descending orbit and negative values in the ascending orbit. Contrarily, south-western and western facing slopes show positive deformations in the ascending orbit or negative values in the descending one.
In the affected area of T , icău landslide, the mean velocities of the PSs have similar behavior with the above-described pattern, with high mean velocity values from ±7 mm/y up to ±24 mm/y (Figure 6b,d). The identified points are located mainly in the upper and the middle parts of the slope, as the lower part of the hillslope is covered by dense vegetation. Specifically, stable points are found at the top of the hill and above the landslide crown, while the unstable points are over the basal part of the landslide scarp and the depleted area.
The projected velocities along the slope show ( Figure 7) similar patterns with the 1D velocities. The up and E-W components also confirm the conclusions obtained from the ascendent and descendent 1D velocities. Interestingly, these components are not as noisy as 1D data, the landslided hillslope being distinguishable from the adjacent stable hilltop. In the same time, coverage is lost because some points from the ascendent and descendent orbits are not in close proximity.
In the case of ENVISAT and ERS-1/2 SAR data, analyzed with the SBAS technique, the results are relatively poor compared to Sentinel 1 outputs, especially for the ERS-1/2 dataset where the low coherence of the signal significantly reduced the spatial information. The mean velocity map of the ENVISAT ascending dataset (Figure 6e,f) shows values ranging from −4 to −9 mm/year over the investigated event, and it confirms the continuous slowmoving state of the landslide. Due to the much better results of the Sentinel 1 datasets, we will focus our analysis on these results.
The Sentinel 1 SAR data processing allows the detection of landslides' affected areas, updating and improving the available database, and monitoring the critical areas (Figure 6a,c). Most of the detected deformations overlay the landslide-affected areas within the landslide hazard zonation of Ias , i Municipality delineated by the Ias , i Municipality Administration (green polygons in Figure 8) and which cover 61.45% of the city surface. Seven out of the nine hazard zones delineated by the Ias , i Municipality Administration are affected by slope displacements. we measure negative velocities in the ascending geometry while the descending records are positive ( Figure 6).  The statistical post-processing of the SAR velocity pinpoints the present-day clusters of displacements related to active or dormant landslides mapped by [62,79] (Figure 8orange and magenta polygons). We considered the empirical thresholds of [63,103] for landslide velocity to classify the rate of movement of the detected scatterers as follows: >16 mm/year are the velocities of very slow landslides that require maintenance, while in contrast, velocity <16 mm/year is characteristic of extremely slow landslides. Hence, two clusters of moving points with a velocity higher than 16 mm/year arguably classify as very slow-moving landslides located on hazard zones 9 and 2 of Figure 8.  Specifically, these clusters of displacement that we consider hot-spots of landslide activity are the red polygons in Figure 8 over the Țicău Neighborhood, in the North-Eastern side of Copou Hill and the NE hillslope of Galata Neighborhood. As highlighted in Figure 8, some other localized scatterers with similar velocities of the same class are associated with the construction of recent residential buildings over existing landslides. Such areas include the Munteni Neighborhood on the Western Copou Hillslope between Viticultori and Dealul Zorilor streets-hot-spot c; Păcurari Neighborhood on the S Copou Hill between Șipoțel and Cazărmilor streets-hot-spot d; Bucium Neighborhood on the E Socola hillslope between Mihail Galino and Margareta Baciu streets-hot-spot e; W Păun hillslope between Ciprian Porumbescu and Păun streets-hot-spot f. Other cases are associated with site-specific activities occurring over the hillslopes (e.g., Păcurari Neighborhood on S Copou hillslope, an area with intense activity due to the new houses built between Iancu The availability of both orbits, ascending and descending, improves our results' reliability through cross-validation of the measurements obtained by the first geometry with the results of the second one and by the possibility to compute the components of the LOS displacement. Knowing that the T , icău landslide is sliding from southwest to northeast, we measure negative velocities in the ascending geometry while the descending records are positive ( Figure 6).
The statistical post-processing of the SAR velocity pinpoints the present-day clusters of displacements related to active or dormant landslides mapped by [62,79] (Figure 8-orange and magenta polygons). We considered the empirical thresholds of [63,103] for landslide velocity to classify the rate of movement of the detected scatterers as follows: >16 mm/year are the velocities of very slow landslides that require maintenance, while in contrast, velocity <16 mm/year is characteristic of extremely slow landslides. Hence, two clusters of moving points with a velocity higher than 16 mm/year arguably classify as very slowmoving landslides located on hazard zones 9 and 2 of Figure 8.
Specifically, these clusters of displacement that we consider hot-spots of landslide activity are the red polygons in Figure 8 over the T , icău Neighborhood, in the North-Eastern side of Copou Hill and the NE hillslope of Galata Neighborhood. As highlighted in Figure 8, some other localized scatterers with similar velocities of the same class are associated with the construction of recent residential buildings over existing landslides. Such areas include the Munteni Neighborhood on the Western Copou Hillslope between Viticultori and Dealul Zorilor streets-hot-spot c; Păcurari Neighborhood on the S Copou Hill between S , ipot , el and Cazărmilor streets-hot-spot d; Bucium Neighborhood on the E Socola hillslope between Mihail Galino and Margareta Baciu streets-hot-spot e; W Păun hillslope between Ciprian Porumbescu and Păun streets-hot-spot f. Other cases are associated with site-specific activities occurring over the hillslopes (e.g., Păcurari Neighborhood on S Copou hillslope, an area with intense activity due to the new houses built between Iancu Flondor, Nicolae Oblu and Bucovina streets-hot-spot g; Cas , in street-hot-spot h; Ciric Sports Base built on an excavation site upslope the landslide scarp-hot-spot i).
The other areas that show ground deformations are classified as extremely slowmoving landslides. As the lower boundary of this class, we considered the 2 mm/year threshold representative for the investigated area. In our analysis, areas with velocities lower than 2 mm/year, classified as relatively stable, are located on the plateaus of the hills (Figure 8-stable areas j and l) and floodplains (Figure 8-stable areas k, m, and n).
The use of the Multi-temporal InSAR techniques enables the possibility to discriminate between the landslide-induced deformations over large areas from the isolated cases unrelated to the process. The former case is associated with the mapped landslide ( Figure 8-the yellow polygons), while the latter relates to single or clusters of buildings located on slopes unaffected by landslides (red and blue squares not included in the yellow polygons in Figure 8). Other deformation hot-spots, usually on floodplains, are associated with residential or non-residential planning projects (blue polygons in Figure 8). We used the Mann-Kendall test to compute the trend for all the PS points in both ascending and descending orbit to discriminate between these hot-spots. Only the PS points with MK test p-value smaller than 0.0000005 (a trend is present) and a slope of the linear trend higher/smaller (depending on the orbit geometry) than 0.15/−0.15 were selected (trend present with a slope that indicates consistent deformation). Thus, were sorted out 1288 points in the ascending orbit and 951 points in the descending orbit (red and blue squares in Figure 8). The remaining points are classified as extremely slow-moving landslides delineated by the yellow boundary or isolated instabilities seen as single features in Figure 8.
The Bucium East Hazard Zone (hazard zone 6) is one of the city's growing neighborhoods currently characterized by the expansion of the residential buildings. It has one very slow-moving landslide hot-spot and three extremely slow-moving landslides located over the previously mapped landslides. The Bucium West Hazard Zone (hazard zone 5) has one very slow-moving landslide hot-spot and four extremely slow-moving landslide hot-spots, similar to Bucium East Hazard Zone. The Galata Hazard Zone (hazard zone 2) has three hot-spots of very slow-moving landslides and six hot-spots of extremely slow-moving landslides. Most of them are on the NW hillslope of Galata Hill, where many residential buildings exist. The Copou West-Păcurari Hazard Zone (hazard zone 1) has four hot-spots of very slow-moving landslides and eleven hot-spots of extremely slow-moving landslides, known as recent landslides, and inhabited by the neighborhood's community. The Copou East Hazard Zone (hazard zone 9) has the largest hot-spot identified within a slow-moving landslide and four hot-spots of extremely slow-moving landslide hot-spots. The Aviat , iei-Airport-Moara de Vânt Hazard Zone (hazard zone 7) has one slow-moving deformation hot-spot and three hot-spots of extremely slow-moving landslides, the Central-Tătăras , i Hazard Zone (hazard zone 8) has only one hot-spot associated with extremely slow-moving landslides, and the other two hazard zones, namely the Cetăt , uia and Manta Ros , ie Hazard Zones (hazard zones 3 and 4) do not present landslide-related deformation hot-spots.
Based on the PS time series displacements, we can argue the sensitivity of slope dynamics to rainfall. The evolution pattern of the PS measurements over T , icău landslide is defined by trends associated with accelerating or decelerating periods (Figure 9). The displacements recorded during the investigated period confirm that T , icău landslide is active with rates of displacements that exhibit different trends separated by generalized breakpoints for the considered PS points. that groundwater response to rainfall is susceptible to a time delay [104], our interpretation of rainfall's influence on displacement through groundwater recharge is the most probable. The complexity of the topography (roads, embankments, built surfaces) in the study area could explain the multiple breakpoints of the displacement trend during 2015 and 2016. That explains why some of the surface scatterers have a different sliding response in time. The displacement behavior related to the distribution of the daily and the cumulative rainfall amount shows that prolonged rainfalls of medium intensity and two to three weeks long are much more significant for displacements initiation and acceleration compared to heavy rainfalls. Even though the short period of heavy rainfalls, of about one to two days, produces a more significant total amount of precipitation, it mainly occurs during drought periods, and due to the high intensity, the runoff of the water is more significant than the infiltration. In this context, PS displacements can be monitored together with the rainfall time series to search for increasing displacement trends that might indicate a possible reactivation event. The reactivation of 1942 is one such example that took place during the spring season. In that case, the landslide reactivation occurred due to the high accumulated amounts of precipitation. The rainy period started in the previous The analysis of precipitation series indicates for the years 2014, 2016, and 2017 rainy springs and summers (both daily and as 30 days cumulative rainfall), while 2015 is drier during all seasons. In this context, the displacement trend decelerates from 2014 to 2015 and accelerates again during 2016. The location of the trend breakpoints (the rug plots from Figure 9) is consistent across all the PS for the spring 2015 and autumn 2016 trend changes. While there is no groundwater monitoring in the study area, and considering that groundwater response to rainfall is susceptible to a time delay [104], our interpretation of rainfall's influence on displacement through groundwater recharge is the most probable. The complexity of the topography (roads, embankments, built surfaces) in the study area could explain the multiple breakpoints of the displacement trend during 2015 and 2016. That explains why some of the surface scatterers have a different sliding response in time.
The displacement behavior related to the distribution of the daily and the cumulative rainfall amount shows that prolonged rainfalls of medium intensity and two to three weeks long are much more significant for displacements initiation and acceleration compared to heavy rainfalls. Even though the short period of heavy rainfalls, of about one to two days, produces a more significant total amount of precipitation, it mainly occurs during drought periods, and due to the high intensity, the runoff of the water is more significant than the infiltration. In this context, PS displacements can be monitored together with the rainfall time series to search for increasing displacement trends that might indicate a possible reactivation event. The reactivation of 1942 is one such example that took place during the spring season. In that case, the landslide reactivation occurred due to the high accumulated amounts of precipitation. The rainy period started in the previous year's autumn and continued in the spring of the year after coupled with the snow melting. These factors led to the soil's complete saturation due to the water level increase and finally to the landslide triggering.

Numerical Simulations
We carried out a parametric study to verify the possible changes of landslide dynamics due to variations in the material properties and optimize our model. In this way, we reduce the uncertainty of the input data and minimize our assumptions to set up a more realistic model. In the simulations, we stressed the landslide unit (Layer 6-Landslide debris) because the displacements expected to occur are within this layer as they have been detected on the surface through MT-InSAR techniques and field evidence as well. We also took into account the water level measured during the execution of the boreholes even though it is susceptive to temporal depth variations. Thus, we performed several trial-and-error simulations to assess the response of the landslided material to the variation of parameter values. Although the properties are only slightly changed, the landslide behavior varies significantly.
In terms of displacement magnitude (Figure 10), the model behaves consistently with the increment of material properties values. Specifically, by increasing the values of friction angle (φ) and the effective cohesion (c ) of material, the total displacements recorded by the model are decreasing. The simulation considered representative for the investigated landslide (φ = 20 • , c = 5 kPa, t = 1 kPa) records displacements of a maximum of 2.5 cm and points out the middle sector of the slope as the most prone to deformations, compared to the upper and lower ones. These results significantly increase our understanding of the landslide in-depth behavior and sliding mechanism and can be used as prior information for future stabilization and development works.
Analyzing the displacement magnitude field, the middle sector of the slope results as the most active part of the landslide regardless of the scenario. The displacement magnitude history of some points distributed along the landslide body ( Figure 11) supports the described behavior also. Moreover, the behavior of in-depth history points indicates higher displacements as we get closer to the surface.
To determine the landslide typology and mechanism, we used the maximum shear strain increment, which is a good indicator of the failure mechanism [105]. The identified shearing zones are associated with the sliding surface, the main and secondary scarps along the landslide body ( Figure 12). The secondary scarps increase in magnitude after they are subjected to the water level increase. In fact, in conditions of high-water level, similar to the 1942 reactivation, both magnitude of displacements and maximum shear strain increment significantly increased compared to the scenario where the water table was considered the one measured in the boreholes ( Figure 13). In this"worst" scenario, the total displacements are much higher than 4.5 cm, reaching a maximum of 9 cm in the middle sector. Based on the maximum shear strain increment, the sliding surface and the secondary scarps aggravate, and the main scarp is susceptible to reactivation. Although it is not possible to directly compare the velocities measured through interferometry techniques and the displacements obtained in the numerical simulations, we observe that the results of both methods point out the median part of the landslide as the most active sector.
The advantage of numerical modeling consists of the continuous displacements field that we obtain, while in the case of PS measurements, we get discontinuous data but are highly accurate. On the other hand, performing numerical modeling requires as many insitu measurements as possible. The lack of such data limits the time-dependent modeling to observe the evolution of displacements, comparable to the temporal displacement series of MT-InSAR results. Because of this, we performed our simulations using the Mohr-Coulomb elastoplastic model. The results give us information about the displacements' behavior, their magnitude, and the sliding mechanism.

MT-InSAR Outputs and Numerical Simulations
The investigated Țicău landslide behaves as a very slow-moving translational slide from SW to NE direction. The typology of the landslide is concluded based on the simulations we carried out. The maximum shear strain increment reveals the failure surface Although it is not possible to directly compare the velocities measured through interferometry techniques and the displacements obtained in the numerical simulations, we observe that the results of both methods point out the median part of the landslide as the most active sector.
The advantage of numerical modeling consists of the continuous displacements field that we obtain, while in the case of PS measurements, we get discontinuous data but are highly accurate. On the other hand, performing numerical modeling requires as many insitu measurements as possible. The lack of such data limits the time-dependent modeling to observe the evolution of displacements, comparable to the temporal displacement series of MT-InSAR results. Because of this, we performed our simulations using the Mohr-Coulomb elastoplastic model. The results give us information about the displacements' behavior, their magnitude, and the sliding mechanism.

MT-InSAR Outputs and Numerical Simulations
The investigated T , icău landslide behaves as a very slow-moving translational slide from SW to NE direction. The typology of the landslide is concluded based on the simulations we carried out. The maximum shear strain increment reveals the failure surface and the presence of several minor scarps that accentuate as the water level increases. The mean annual velocities, depicted through the processing of SAR data, measure from 10 to 20 mm/year over the landslide body. Even though the measurements are in the Line of Sight of the satellite, it is a good indicator of the deformations that occur in the area, especially the ascending orbit case. Due to the satellite's flight path and the sliding direction of the landslide, the measured displacement values are close to the actual displacements that occur over the landslide.
Simulations of slope displacements reveal the landslide mechanism and allow the investigation of multiple scenarios to assess the process's dynamics properly. The pattern of displacements is consistent with the simulations we performed, intending to reduce the uncertainties related to the data we used.
The results of both used techniques are in agreement and reveal that the middle sector of the slope is the most active part of the landslide. The analysis of the displacement time series indicates that acceleration periods take place during long-term rainfalls. The simulations also enhance this aspect, as they clearly show an increase of total displacements of the landslide and a much more prominent sliding surface.
The numerical modeling results regarding the displacement pattern of the landslide argue the landslide typology as a shallow translational body. This outcome is contrasting to what previous studies suggest, a deep rotational landslide [100]. The borehole information on which these previous conclusions are based (i.e., granulometry, unit weight) is deficient and unconvincing. They are evaluated only from an uncomplete geotechnical point of view without considering the stratigraphic and geomorphologic data (bedrock, surficial deposits, and landslide sliding surface).
A better investigation of the subsurface is necessary to improve the model from a geophysical and geotechnical point of view and in-situ measurements by installing necessary instrumentations. These will allow the use of an advanced model that considers the materials' viscosity characteristics and thus proceed to a time-dependent creep analysis.

Benefits and Limitations of the Analysis
To carry out our investigation, we exploited both the employed techniques at best by integrating their results to overcome the flaws of one method with the strengths of the other. Regarding the MT-InSAR analysis, the main problem encountered during the image processing is related to dense vegetation and loss of coherence for the analyzed period. When coupled with the lack of images, this leads to substantial loss of spatial information, as in the ERS and ENVISAT data stacks. The basal sector of the T , icău landslide is covered by forest, and we can notice that it was not possible to detect deformations over that area, not even with Sentinel-1. However, we overcame this problem using numerical simulations, which indicate that the slope's basal sector is also moving even though it is not as active as the middle part.
Due to the lack of reliable in-situ data for our numerical modeling, we had to make some assumptions worth mentioning to complete our material properties database. We considered the landslide debris as a singular, homogenous entity. In reality, the remolded material most likely has lithological and implicitly geotechnical variations, both in-depth and laterally. To reduce the uncertainty of made assumptions, we performed trial-anderror simulations of various combinations of the intrinsic properties, showing that these assumptions are relevant and justified.
The groundwater level we employed in our modeling is the measured one during drilling, which might be at its lowest depths, considering that the boreholes were executed during the summer. We should keep in mind that it has seasonal fluctuations, manifesting as accelerations and decelerations of deformations, as suggested by the MT-InSAR measurements.
At the same time, in our analysis, we did not consider the anthropic implications such as the existence of the residential buildings that exert an additional load over the unstable body nor the already existing stabilization infrastructure, such as drainage system, retaining walls, and gabions at the base of the main scarp. About this aspect, it is necessary to get a detailed investigation about the type of building foundations, the year of construction, to evaluate the state of degradation, the height of buildings, the construction material.
However, the numerical modeling results are supported by MT-InSAR results and field information, indicating that the displacements obtained are consistent and validate our model. Further geophysical prospects will improve the geological model and will help to generate an improved numerical model that could better explain the spatial and temporal variations of the displacements patterns.
Nonetheless, our results prove that the use of MT-InSAR techniques and SAR data coupled with numerical simulations has the potential to deliver critical information about landslide processes in areas where the lack of in-situ data is the major shortcoming, not only in the well-investigated ones. The MT-InSAR provides data about the spatial extent and the displacements behavior in time, while the simulations give a better understanding of the landslide type and its mechanism and in-depth behavior. Their outcome is beneficial not only for the landslide monitoring and landslide hazard assessment but also as prior information for future planning and policy implementation.

Conclusions
The application of Multi-temporal SAR Interferometry techniques (MT-InSAR) and numerical modeling of slope failures in urban areas is a must in the current context of urban territorial expansion to successfully identify and monitor possible slope deformations that might threaten the infrastructure and people. For the case of Ias , i Municipality, the benefits of coupling MT-InSAR results (i.e., surface velocity maps) and slope numerical simulations allowed us: (1) to identify the potentially dangerous areas affected by slope instabilities; (2) to properly delimitate the extent of the active sector tampering the T , icău neighborhood's integrity; (3) to monitor the landslide activity and its behavior related to the rainfall amount and water level changes; (4) to argue the landslide typology based on its velocity and failure mechanism.
The integration of InSAR and simulation results coupled with field surveys highly facilitates the characterization and improves our understanding of the landslide dynamics and evolution. Thus, we can say that T , icău landslide behaves as a very slow-moving landslide with a translational mechanism sliding from SW towards NE. The displacement pattern points out the middle sector of the slope as the most prone to deformations that are accelerating after extended periods of rainfall. The mean annual velocities over the landslide body vary from 10 mm/year to over 20 mm/year in the middle sector of the slope underlined by the material simulations. These outputs should be deemed as prospective and used as prior information to support the implementation of stabilization policies, as the local authorities recently approved funding for such strategies. Nonetheless, our results enhance the idea of integrating the results of MT-InSAR techniques and numerical modeling as it has the potential to be used in areas lacking ground measurements. Their capabilities to acquire and generate information about slope instabilities are valuable to assess the landslide hazard.
Future works include continuing the landslide monitoring using MT-InSAR techniques and, especially, improving the numerical model. To this end, we will consider the buildings and infrastructure that load the sliding material. There are plans to carry out shallow seismic surveys and electrical resistivity tomography (ERT) to increase the quality of in-situ data. The access to better data and knowledge will allow, in the long run, to model the entire slope and use an advanced viscous model to assess the time-dependent behavior of the landslide.