Surface Instability Mapping in Alpine Paraglacial Environments Using Sentinel-1 DInSAR Techniques

Current climate warming leads to widespread glacier shrinkage in high alpine terrains and associated changes in surface dynamics of deglacierized environments. In consequence, slope instabilities increasingly develop along retreating glaciers through debuttressing effects or degrading permafrost conditions. In the context of associated hazards to the local environment and infrastructure, a thorough analysis of slope instabilities is highly relevant. Affected regions are mostly inaccessible and cover large areas, therefore remote sensing techniques such as differential interferometric synthetic aperture radar (DInSAR) are valuable tools to monitor surface movements and assess their evolution. We apply standard and advanced DInSAR methods using Sentinel-1 SAR data from 2015 until late 2021 to map and classify slope instabilities in three glacierized regions in the European Alps. The final products include an inventory per region, with a total of 815 mapped slope instabilities, of which 38% move <3, 9% move 3–10, 42% move 10–30, and 11% move >30 cm/yr. An additional assessment of four landslides occurring along shrinking glaciers shows time series with recent accelerations in 2018/19. Validation of Sentinel-1 derived slope movement products is performed by comparison with shorter wavelength TerraSAR-X and optical Sentinel-2 derived data using offset tracking. Results clearly show the suitability of Sentinel-1 DInSAR methods to detect a range of slope movements in high alpine terrain, yet also highlight the limitations. We therefore recommend a combination of advanced Sentinel-1 DInSAR and Sentinel-2 offset tracking methods to develop a comprehensive inventory of alpine slope motion.

Digital Object Identifier 10.1109/JSTARS.2023.3287285 lakes continuously develop. The exposed, generally rough topography and steep lateral moraines in glacial and paraglacial environments provide a predisposing setting for slope instabilities of various kinds, including landslides and rockfalls. These often form through glacier debuttressing effects [7], [8], [9] or progressively degrading permafrost conditions and rock glaciers [5], [10]. The detection and analysis of slope instabilities developing in such environments is therefore increasingly important [3], [11], [12]. Glaciers are generally hardly accessible and in situ observations of surface movements are therefore limited. As a consequence, remote sensing methods are necessary to ensure a systematic analysis of slope instabilities over large areas. Effective methods to map and analyze surface movements from space include differential interferometric synthetic aperture radar (DInSAR) as well as optical and radar offset tracking methods. Multiple studies applying such techniques have been performed worldwide and for the European Alps in particular [9], [14], [15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25]. These techniques enable the detection of a wide range of slope movements that occur at scales of small slopes to entire mountain flanks, including rock and soil slides, debris flows and rock slope deformations, as well as rock glaciers and other permafrost-related movements [21], [24], [26], [27]. Involved displacement rates show a similarly broad range, with very slow landslides moving few millimeters per year to rapid movement of several meters per year observed on glaciers, rock glaciers and very fast-moving landslides (e.g., the Moosfluh landslide in Switzerland [9]). This large variability, together with the widespread occurrence of movements in all orientations, calls for the integrated use of remote sensing methods.
Since the operational phase of the European Space Agency's (ESA) Sentinel-1 SAR and Sentinel-2 optical satellites in 2014/15, DInSAR and offset tracking techniques are particularly efficient and suitable for the detection of ground motion. The free distribution of Sentinel-1 (S1) and Sentinel-2 (S2) images allows a continuous and systematic analysis of surface movement in alpine regions, with S1 providing images every 6 to 12 days and S2 every 5 days. SAR techniques, however, remain affected by certain limitations. These are related to the characteristics of the sensor, such as flight direction, wavelength, and resolution, as well as the imaged area, e.g., geometry and land cover [28]. A systematic overview of the possibilities and limitations of SAR techniques in high alpine regions and an assessment of their suitability for systematic screening purposes prior to potential slope failures is crucial when studying the effects of changing  climate conditions on alpine surface dynamics. Particularly in regions of dense infrastructure and numerous settlements, as in the steep topography of the European Alps, such an assessment is becoming increasingly important.
In the following, we present a systematic assessment of available SAR methods to detect surface movement in glacial and paraglacial environments. Specifically, we apply multiple standard and advanced DInSAR techniques using S1 data to develop inventories of slope instabilities in three selected study regions of the European Alps. Based on the results, we outline the possibilities and limitations of our methodology and the suitability of S1 data in high alpine terrain. We compare results to higher resolution TerraSAR-X and optical S2 offset tracking results for accuracy assessment and validation. The work is performed as part of the AlpGlacier project of the ESA (https://alpglacier.geo.uzh.ch/; last accessed 19/01/2023). We focus on paraglacial environments in the Mattertal region in Switzerland, the Mont-Blanc massif region in France and Switzerland, and the Ötztal Alps of Austria (see Fig. 1, Table I). All regions are characterized by steep mountain relief, occurrence of permafrost, abundant infrastructure, and surface displacements. Each site has experienced glacier retreat since the end of the LIA, leading to the potential formation of glacier lakes and unstable slopes. Glacier outlines exist for periods since the LIA [1], [29] and from 2015/16 [2]. Systematic S1 and S2 images at a high temporal resolution are available for all regions.

II. STUDY REGIONS
The Mattertal region in Switzerland is characterized by high and steep relief, introducing a high potential for mass movements. Permafrost is widespread at elevations higher than approx. 2 500 m a.s.l., particularly across North-facing slopes (Alpine Permafrost Index Map, APIM) [30]. Glacier lakes developing in the last century have been studied and inventoried by [31] and show a wide distribution. This correlates with a  II  PROCESSED DATA FOR ASCENDING AND DESCENDING S1 TRACKS COVERING THE THREE STUDY REGIONS profound glacier retreat since 1850 and increased glacier velocities [32], [33]. Numerous studies have been conducted on the region to assess surface movements [16], [26], [34]. An inventory of surface movement at regional scale for the canton of Valais in Switzerland shows a large range of movement types and sizes in the Matter valley (https://sitonline.vs.ch/dangers/ dangers_geologiques/de/; last accessed 19/01/2023). However, this survey is not exhaustive and focusses mainly on regions with infrastructure and settlements.
The Mont-Blanc massif in France, Switzerland, and Italy includes the highest mountain peak in Europe, the Mont Blanc at 4 810 m a.s.l., as well as some of the largest glaciers, e.g., Mer de Glace and Argentière. The region is characterized by a very high and steep relief with extensive snow and ice cover. Permafrost is widespread at elevations above approx. 2 400 m a.s.l. (APIM) [30] and glaciers show a pronounced retreat since the LIA [1], [32]. While glacier lakes have not been extensively mapped in this region, numerous studies focus on rockfall and surface instabilities related to changes in permafrost conditions and thinning of glaciers [35], [36], [37].
The study site in the Ötztal Alps in Austria is also marked by a high relief and potential for mass movements. Permafrost occurs above approx. 2 400 m a.s.l. (APIM) [30]. An inventory of glacier lakes above 1 700 m a.s.l. compiled by [38] shows a widespread distribution of lakes. Glacier outlines from 1998 [29] and 2015/16 [2] indicate a significant glacier retreat of 66% between the LIA and 2006 [39]. Regional studies have been conducted on glacier dynamics, the formation of glacier lakes as well as the development of slope movements and associated hazards [13], [40], [41].

A. Synthetic Aperture Radar (SAR)
ESA's S1A/B single look complex (SLC) interferometric wide (IW) swath TOPS mode data is used for the DInSAR analyses presented in this study. The three study regions are covered by the tracks and acquisition dates listed in Table II. Images are generally available from late 2014/early 2015 until present, with acquisition time intervals of 12 days until early 2017, 6 days until the end of 2021 and again 12 days until present due to the failure of Sentinel-1B. We therefore consider data until late 2021. Sentinel-1 operates at C-band with a wavelength of 5.6 cm, ensuring a spatial resolution of a few meters and the detection of surface movements of up to 84 cm/yr for 6-days repeat intervals, with an expected subcentimetric accuracy [28].
In addition, we consider data from higher resolution and shorter wavelength X-band TerraSAR-X and CosmoSkyMed sensors. Data from these sensors are only available for the Swiss study regions. They typically enable the detection of slow as well as fast movements, which cannot be accurately measured at C-band. Furthermore, their use allows us to validate the results obtained with S1. TerraSAR-X (TSX) has a revisit time of 11 days (minimum) and a wavelength of 3.1 cm with a very high spatial resolution, enabling the detection of small, slow surface movements up to 25 cm/yr. Cosmo-SkyMed (CSK) also has a wavelength of 3.1 cm and detects movements up to 70 cm/yr with a revisit time of 4 days. Both X-band sensors provide data preceding S1 (i.e., before 2014) and thus allow an analysis of past behavior of landslide dynamics.

B. Optical
In addition to SAR data, optical data from S2 are used for validation purposes in all three study sites. Optical sensors overcome some of the limitations associated with SAR sensors and are thus a valuable complementary data source. S2 operates with a 5-days revisit time and ensures a spatial resolution of 10 m, with expected sub-pixel accuracies of approx. 1 m/yr. Optical sensors can be considered for the analysis of very fast motion of large landslides.

C. Digital Elevation Models
Digital elevation models (DEMs) are used during processing and for geocoding of the SAR and optical data and results. For the Swiss Mattertal region, the available SwissAlti3D DEM from

IV. METHODOLOGY
DInSAR techniques are widely applied for the detection of surface displacement at millimeter to subcentimeter accuracies and at local to regional scales [7], [9], [17], [19], [21], [22], [24], [27]. To detect slow and fast slope movements of various sizes, multiple DInSAR techniques are necessary (see Table III). Standard DInSAR involves the generation of differential interferograms, which provide a good initial overview of surface displacement. To generate displacement maps and time series of velocity, advanced DInSAR methods are required. These include multi-looked temporal interferometry (MTI) techniques such as stacking and multi-baseline [27], [42], [43], [44], which show surface movements of several tens of centimeters per year. Additionally, persistent scatterer interferometry (PSI) is applied to highlight slow movements of few millimeters per year [45], [46], [47], [48]. Standard and advanced DInSAR methods are thus equally required to ensure the detection of a wide range of slope movements and assess their temporal and spatial evolution.

A. Standard DInSAR Processing
Standard DInSAR combines two SAR images acquired at different times and covering largely the same area with only a slight offset in orbit [49]. The difference in phase information of the two combined images contains a topographic and interferometric phase component. The topographic phase is removed through the use of an external DEM. The remaining interferometric phase reflects surface movement along the satellite's look direction, i.e., line-of-sight (LOS), and is displayed as colored fringe cycles in the resulting differential interferogram [7], [28], [49]. One fringe cycle represents movement on the order of a 2π ambiguity, relating to half the satellite wavelength, i.e., 2.8 cm for S1. Fringe signals in differential interferograms thus give a good first indication of surface movement. To retrieve an absolute deformation map, the modulo 2π ambiguity of the interferometric phase is resolved through phase unwrapping by means of a statistically based minimum cost flow approach [50], [51]. Particularly in areas of high relief, phase unwrapping often introduces errors and requires pixel averaging and coherence thresholds to yield reliable results [44], [51]. Additional phase contributions are given by noise (leading to decorrelation and coherence loss) and atmosphere [28]. Such effects are corrected and removed through adaptive and atmospheric filtering [51], [52]. The final unwrapped and filtered differential interferograms act as input for MTI stacking and multi-baseline methods. They are thus a requirement for advanced DInSAR techniques.
Standard DInSAR techniques achieve a relatively high accuracy of few millimeters to centimeters on detectable displacements. The accuracy is mainly limited by the atmospheric and phase noise contributions (see Table IV). As a typical atmospheric distortion, we consider a relatively high phase error of π/2, which corresponds to 0.7 cm for S1. In differential interferograms of 6 or 24 days this corresponds to 43 or 11 cm/yr, respectively. Phase noise is related to coherence, where the error attributed to noise for high-coherence areas is quite low (e.g., 1-2 mm) [15]. A total error of minimum 1 cm can thus be expected for displacement rates extracted from S1 differential interferograms. This relates to roughly λ/8, which for 6-days S1 differential interferograms with a spatial resolution of a few meters gives an accuracy of approx. 40 cm/yr [53], [54].
We generate differential interferograms according to two pairing strategies. The first strategy involves interferometric pairs based on the shortest possible temporal baseline, i.e., typically 6 or 12 days. These are expected to show fast movements of up to approx. 80 cm/yr. The second strategy involves pairs with 24-days temporal baselines to highlight slower movements. Few pairs with temporal baselines of one to two years are also generated as a basis for the multi-baseline technique (cf. Section IV-C). However, such long temporal baselines generally lead to strongly decorrelated results in mountainous regions.

B. Stacking
Stacking is an MTI technique that combines multiple unwrapped differential interferograms of selected temporal baselines over a specified time frame. The results show the average  [42], [43]. Through this integration, the signal-to-noise ratio increases and true surface deformation is better highlighted. Atmospheric disturbances and artefacts with random variability are removed [28], [44], [55], [56]. The stacked images are referred to a specified point of zero deformation to avoid phase bias and are processed according to a least-squares solution [42]. Time series are extracted from the final displacement map to show the spatio-temporal evolution of surface movements. Stacking results thus yield more information on slope instabilities than when only considering standard differential interferometry.
Stacking is suitable for seasonal and annual deformation analyses, where changes in average surface displacement can be highlighted from one time period to the next. The expected accuracy is derived as for standard DInSAR but is also dependent on the number of input interferograms [15]. When considering in-series input data, i.e., subsequent differential interferograms of the shortest temporal baseline, the estimated error increases with an increasing number of input images and is on the order of 1-2 cm/yr (see Table IV). We test stacking procedures for three different time intervals, considering differential interferograms of the shortest possible temporal baseline (6 or 12 days) and of 24 days as input. In a first approach, we stack all pairs during 2015-2021 to show the average displacement throughout this period. The second approach considers pairs from the summer seasons (July-October) of 2015-2021 combined to highlight the average summer displacement during these years. The third approach considers pairs from each summer separately to assess differences between each year and generate a time series of average summer velocity. The constraint on data from July to October ensures higher coherence due to avoidance of snow cover and snow melt images. We additionally apply a coherence threshold to mask out all areas with a coherence below 0.4.

C. Multi-baseline Interferometry
Multi-baseline, or small baseline subset interferometry, is equally an MTI technique and similar to stacking. Multi-baseline also involves the combination of multiple unwrapped differential interferograms to mitigate noise and highlight true surface movement in resulting displacement maps and time series. In addition, however, Multi-baseline enables the generation of progressive time series during individual summer periods, which is not achieved in stacking results. Multi-baseline LOS displacement maps show the progressive change of surface motion in a single reference approach, i.e., for time periods always considering the first input date with all subsequent dates (A-A, A-B, A-C, A-D, etc.). It thus complements stacking results by enabling such selective time series. The input interferograms, however, must be temporally connected and cannot have any gaps during the covered time period. Furthermore, this technique requires the use of small spatial baselines to reduce spatial decorrelation [28], [44]. These are given by the generally short spatial baselines of no more than 200 m of S1 data. The expected accuracy is on the order of 1-2 cm/yr as for stacking (see Table IV).
Multi-baseline analyses are suitable when studying seasonal deformation patterns but can also be applied to detect long-term displacement. We test two approaches, each covering a different time period and based on differential interferograms of the shortest temporal baseline. In a first seasonal approach, we consider each summer separately, i.e., July to October of 2015 to 2021, to derive seasonal time series. We select the same input interferograms as used during stacking to enable a direct comparison of the results. To obtain the average summer displacement during 2015-2021, we select and stack the maximum displacement maps of each summer. These typically cover the entire coherent summer period, i.e., from the first input date (beginning of July) to one of the last input dates (end-September to end-October). In a second multi-annual approach, we consider the summer seasons of all years combined, i.e., covering July to October from 2015 to 2021, to understand the longer term seasonal evolution. In order to maintain the required temporal connectivity, we additionally consider few differential interferograms of one to two years temporal baselines to link each summer to the next. Use of these is expected to lower the overall coherence but highlight slow movements on the order of few millimeters to centimeters per year. Areas of low coherence are again masked out, while remaining large scale atmospheric errors are spatially filtered.

D. Persistent Scatterer Interferometry
PSI enables the detection of slow surface movement of few millimeters per year at millimeter-scale accuracies. The technique is based on the selection of pixels termed persistent scatterers (PS) that exhibit a scattering behavior in multiple differential interferograms and remain coherent over time [45], [57]. PS are scarce in nonurban and vegetated/snow-covered areas but relatively abundant over alpine terrain in summer as well as over settlements and infrastructures. As such scatterers are often smaller than the resolution cell size, coherence can be preserved also when considering differential interferograms of large spatial baselines [44], [45], [58]. PSI can thus retrieve phase displacement with a minimal contribution of other phase components, in particular noise [59]. Additionally, the combination of many differential interferograms minimizes the effects of atmospheric disturbances and increases overall image coherence and accuracy [28], [45], [57]. Furthermore, local displacement can be recognized through single PS, allowing the detection of very small and detailed features within larger slope instabilities that may move at different rates. The PSI technique is thus required in addition to stacking and multi-baseline to detect slow and small movements.
The accuracy of PSI results has been widely studied and is best summarized by findings presented in paper [60]. They estimate a limit to the accuracy for C-band data of 1-1.8 mm/yr over a time period of 5 years for the average displacement rate and of 4.2-6.1 mm/yr for single deformation measurements (see Table IV). Standard deviations, on the other hand, are estimated as 0.4-0.5 mm/yr on displacement velocities.
Our PSI approach follows the interferometric point target analysis processing as developed by Gamma Remote Sensing AG. The technique is based on co-registered SLC data, which are oversampled to ensure better point coverage and quality for the S1 IW TOPS acquisition mode. Persistent scatterer candidates are selected based on their spectral diversity and temporal intensity variability. Various spectral and intensity thresholds are tested to find the optimum solution, with an overall coherence threshold of 0.75. The PS displacement is determined iteratively in two-dimensional regressions using multi-and single-patch approaches. The solutions are unwrapped and atmospherically corrected to yield a final point-based displacement map. We follow a single reference approach, which limits the detectable velocities to few centimeters per year. Again, summer images covering July to October of 2015-2021 are selected, to avoid snow and ice cover and thus preserve sufficient coherence for the point selection. This gives a high point density in mountainous regions despite leading to gaps during winter in the resulting time series. If winter data were included, as is the case for example in the European Ground Motion Service EGMS (https://egms.land.copernicus.eu/), very few scatterers would be detected in mountainous regions, resulting in barely any information on surface movements.

E. Optical and SAR Offset Tracking-Validation and Integration
Offset tracking, also termed feature tracking or digital image correlation, tracks displacement in multi-temporal imagery by comparing pixels or distinct surface features in the input images [9], [20]. The technique is typically based on data acquired with high-resolution optical or SAR satellite sensors. Surface deformation is measured and quantified in two directions, orthogonal to the LOS and near-horizontal. In SAR coordinates, this relates to movement in slant-range and azimuth directions. The use of orthorectified optical imagery retrieves movement in North-South and East-West directions [9], [18]. Offset tracking thus enables the detection of rapid, large, and heterogeneous slope movements in two dimensions and complements and validates the results obtained with the applied DInSAR methods [9], [14], [18], [59], [61], [62].
We use offset tracking based on optical S2 and SAR TSX images. Offset tracking with S2 can detect rapid displacement of up to meters per day, with possible displacement accuracies of 1/10th to 1 pixel. It is restricted to cloud-and snow-free images, while changes in vegetation cover and solar illumination can lead to decorrelation [20], [61]. Using the normalized crosscorrelation method with optimized parameters on S2 images [63], we generate displacement maps showing average annual horizontal movement from one summer to the next for 2016-2021 [64]. Offset tracking based on very high-resolution TSX data ensures a high ground resolution and sampling, with a subpixel accuracy of 1/10th pixel. We use temporal baselines of 11, 22, and 33 days, resulting in an accuracy on the order of a few meters per year [18]. Late summer image pairs with short perpendicular baselines are considered with a time interval of about six years from 2008-2013. To improve the quality of the results, stacking of multiple offset tracking data of similar time intervals is applied. Note that TSX data is only available for Switzerland and thus a validation using this data is not performed for the Mont-Blanc massif and Ötztal study regions.

A. Standard DInSAR Results for the Three Study Regions
Differential interferograms generated for the shortest possible temporal baselines, i.e., 6 or 12 days, generally show good coherence for the summer period from July to October in all three study regions. Examples for the Swiss site are shown in Fig. 2. Glaciers, snow, and vegetation cover are strongly decorrelated as expected. This is particularly the case in the Mont-Blanc massif area, where large extents of the glacierized areas have low coherence and thus yield less reliable results on slope instabilities compared to areas of higher coherence. The overall coherence markedly decreases in differential interferograms with 24-days temporal baselines. This indicates that surface movements in the study regions are too fast to be accurately detected with standard DInSAR during periods of 24 days.

B. Stacking
Stacking results of the first approach, using all available differential interferograms of the shortest temporal baselines from 2015-2021, have a generally low coherence. This is due to the inclusion of winter data. Results from the second approach, i.e., based on all summer images from July to October 2015-2021, show much higher coherence and spatial coverage. Considering  the stacking results from the third approach, i.e., each summer season individually, years with a smaller number of input differential interferograms generally show lower coherence compared to years with a larger stack. This is due to less averaging of random noise and variability. The number of input images per stacking approach, study site and orbit are listed in Table V. Overall, the results clearly highlight slope movements in all three study regions. LOS displacement rates range from few centimeters to up to half a meter per year as detected in combined summer 2015-2021 results (see Fig. 3).
The coherence notably decreases in stacking results based on 24-days differential interferograms. For further analyses and especially a comparison with multi-baseline and PSI results, stacking based on the shortest temporal baselines is thus more reliable.

C. Multi-baseline Products
Seasonal multi-baseline results largely match those from stacking, highlighting the same surface movements with very similar LOS displacement rates. Fig. 4 shows the average summer displacement during 2015-2021, i.e., the stacked maximum displacement maps from each summer. Due to strong decorrelation in the Mattertal region in the summers of 2018 and 2019, we have split each of these two summer seasons into two periods to preserve coherence and the required temporal connectivity. This break in the time series allows us to observe higher velocities at the beginning of summer and a slight decrease toward the end of summer (cf. Section V-E).
Multi-annual multi-baseline results include one-to two-years differential interferograms to link the summer seasons of 2015-2021. They highlight slow movements visible at a scale of +/-1  cm/yr, which were not detected in seasonal multi-baseline results. A comparison with PSI results should ideally show similar slow movements. The use of both multi-baseline approaches thus enables the detection of a wide range of slope movements.

D. Persistent Scatterer Interferometry
The PSI results show a good distribution of persistent scatterers in the Mattertal and Ötztal study regions. Fewer scatterers are obtained for the Mont-Blanc massif, due to the apparent lower overall coherence in the paraglacial areas (see Fig. 5). Slow movements on the order of 1-2 cm/yr are detected. These are mostly confirmed by signals in the multi-annual multi-baseline results.

E. Integration and Inventory
Based on S1 DInSAR results, we develop a non-exhaustive inventory of 815 slope instabilities in the three study sites and classify detected features according to movement type and Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.  velocity. This allows an interpretation of recent movement in the studied regions and an assessment of the suitability of each S1 technique. Slope instability mapping per study region is manually done in QGIS following a grid-wise approach. The movement type is determined based on a geomorphological analysis by studying surface features in orthophotos [21], [24]. This classification reveals that landslides, rockslides, rock slope deformations, and rock glaciers predominate in the studied regions (see Table VI). Approximately 20% of the detected features cannot be classified into a known movement type due to poor orthophoto image quality (cloud/snow cover, low resolution) or lack of characteristic surface features. Slope movements are additionally classified according to two velocity classes, provided by the Swiss Federal Office for the Environment (FOEN) guidelines "Protection against Mass Movement Hazards" [65] and the classification scheme of the International Permafrost Association "Rock glacier inventories and kinematics Action Group" [66]. We perform this classification by visually comparing S1 seasonal summer displacement maps and assigning a possible velocity range for each mapped surface feature. Note that this represents LOS velocity extrapolated per summer and thus an approximation of true surface motion. In some cases, wrapped differential interferograms also give an indication of the occurring surface movement. The majority of mapped slope displacements move 10-30 cm/yr (see Table VII) and are largely detected in stacking and multi-baseline results. In the Swiss and Austrian sites, approx. 12% and 15% of all detected features, respectively, move 30-100 cm/yr. This indicates that advanced multi-temporal DInSAR techniques using S1 data are highly suitable for the detection of relatively slow movements [67] across large regions. In the Mont Blanc region in France, on the other hand, half of all detected features move at rates of no more than 2 cm/yr and are mainly detected with PSI. This highlights that PSI is suitable for the detection of very slow  [67]. The results show that an integration of stacking, multi-baseline and PSI techniques enables the generation of a comprehensive inventory of a wide range of surface movements, with all methods complementing each other. Very rapid landslides, on the other hand, are not well detected with C-Band sensors such as S1 due to decorrelated signals at the available temporal baselines. Assigning a velocity range to such features requires the use of other sensors with shorter time intervals or longer wavelengths.

F. Four Example Landslides in the Three Study Regions
In the following, we show results of the temporal and spatial evolution of four landslides detected along retreating glaciers in the study regions. For each slide, we extract the average LOS velocity for each individual summer of 2015 to 2021 from the corresponding S1 stacking and multi-baseline displacement maps. To derive more representative values of the true displacement field, we correct the velocity under the assumption of largely slope-parallel movement and display the resulting values as time series (see Figs. 6, 8, 10, and 12). The correction is based on the slope orientation and inclination and generates scaling factors for any given location. To validate the S1 time series, we also include average velocities from optical S2 offset tracking results. For the Mattertal region, we additionally display TSX offset tracking and CSK PSI data.

1) Findel Slide, Mattertal:
The Findel slide near Zermatt is located on the northern flank of the Findel Glacier and covers an approximate area of 0.347 km 2 (see Fig. 6) The time series in Fig. 7 shows the temporal evolution for a point at 410787E, 5096059N (UTM 32T; black star in Fig. 6) based on S1 stacking and multi-baseline, S2 and TSX offset tracking and CSK PSI results. Velocity values are given as the summer average for S1, TSX and CSK and annual average for S2. A very good correlation between the S1 DInSAR methods is observed, showing a distinct acceleration in 2018/19, followed by a subsequent deceleration. This trend is confirmed with S2 and TSX offset tracking. CSK results provide an indication of movement prior to S1 data availability, highlighting an overall acceleration of the slide since 2015. Very high velocities of >1 m/yr in 2018/19 lead to decorrelation in S1 results and a break in the time series, resulting in two periods (July and October in 2018, August and September/October in 2019). This break shows faster movement at the beginning of summer and a possible decrease in velocity toward the end of summer, suggesting some interseasonal variability. TSX and S2 results confirm very high velocities of several meters per year and are therefore required in addition to S1 DInSAR to understand the temporal evolution of this slide. The estimated errors for S1 and CSK are very low with only few millimeters to centimeters [15], [60]. S2 errors are on the order of 1/10th of a pixel, i.e., approx. 1 m/yr [68], [69], [70]. TSX errors of 1/10th of a pixel may exceed few m/yr for the considered acquisition intervals [18].
2) Mer De Glace Slide, Mont-Blanc Massif: The Mer de Glace slide in the Mont-Blanc massif is located on the eastern flank of the Mer de Glace glacier and covers an approximate area of 0.497 km 2 (see Fig. 8). The time series in Fig. 9 shows the temporal evolution of a point at 340487E, 5083846N (UTM 32T; black star in Fig. 8) for S1 stacking and multi-baseline and S2 offset tracking results. A very good correlation between S1 results is noted, with a low error of only few millimeters to centimeters. No TSX or CSK data are available for this site. We observe a steady increase in velocity since 2015 with a strong acceleration in 2018 and subsequent deceleration. The overall velocity extracted with S1 is generally lower than for the Findel slide and does not exceed 60 cm/yr. S2 does not show a clear movement signal near the black star in Fig. 8. This is due to observed low velocities of +/-1 m/yr (see Fig. 9), which lie within the S2 error of 1/10th of a pixel. A strong signal, however, is noted to the North (white diamond in Fig. 9), corresponding to an area of no data in the S1 displacement maps. When studying the orthophoto, this area shows a lateral moraine with distinct downslope lineaments that cannot be observed further South. Such surface features indicate surficial movements that may occur due to relatively fast sliding and exceed velocities detectable with S1, resulting in decorrelation and no data in this area. This example strongly highlights the need for both S1 DInSAR as well as S2 offset tracking methods to detect slope motion.

3) Saleinaz Slide, Mont-Blanc Massif:
The Saleinaz slide is located on the southern flank of the Saleinaz Glacier in the Mont-Blanc massif and covers an approximate area of 0.056 km 2 Fig. 6. The Findel slide is visible in the Google Earth satellite orthophoto and all available SAR and optical datasets overlain on the SwissALTI3D 5 m DEM. A velocity time series is extracted at the point indicated by the black star (see Fig. 7), the 2015/16 glacier extent is shown in pale blue with black outlines [2]. The three color scales at the bottom indicate 1) a velocity range of +/-30 cm/yr for Cosmo-SkyMed PSI (CSK) and Sentinel-1 stacking (S1) results, 2) the 2π fringe cycles for Sentinel-1 interferograms (S1 Int) and 3) the velocity range of 0-8 m/yr for TerraSAR-X (TSX) and Sentinel-2 (S2) offset tracking results. CSK PSI shows little motion during 2008-2013. A subsequent acceleration since 2015 is detected with S1 stacking results (top row). An additional strong acceleration in 2018/19 leads to decorrelation in S1 interferograms and stacking and requires use of TSX offset tracking (middle rows). In 2018, S1 interferograms show coherent motion at the beginning of summer (July), but appear decorrelated in August due to very fast displacement. The movement decelerates again toward the end of summer (October), indicating an inter-seasonal change in velocity. A similar pattern is observed in 2019. Overall, the moving area clearly increases in velocity in 2018/19 followed by a deceleration in 2020/21, with the eastern part of the slide showing more prominent and faster movement than the western part. S2 offset tracking results confirm this trend (bottom row).   (Fig. 9), the 2015/16 glacier extent is shown in pale blue with black outlines [2]. The color scales at the bottom indicate a velocity range of +/-30 cm/yr for Sentinel-1 stacking (S1) results and a velocity range of 0-8 m/yr for Sentinel-2 offset tracking (S2) results. S1 stacking results clearly show that the Mer de Glace slide accelerates from 2016 onwards, followed by a deceleration in 2019/20 and renewed acceleration in 2021 (top and middle rows). S2 offset tracking results do not show a signal here (bottom row), due to low displacement rates of <1 m/yr that lie within the S2 error. S2 results rather indicate movement further North (white diamond), an area that corresponds to no data in S1 results, suggesting high displacement rates exceeding1 m/yr. (see Fig. 10). The time series in Fig. 11 shows the temporal evolution for a point at 349829E, 5092730N (UTM 32T; black star in Fig. 10) for S1 stacking and multi-baseline, TSX stacking and S2 offset tracking results. No CSK data are available for this region. S1 and TSX errors are again only few millimeters to centimeters, and S2 errors are approx. 1 meter. A good correlation is observed between the S1 results, apart from a slight offset in 2017. The slide appears to accelerate in 2018, with no valid data points in 2020 and 2021 due to very strong decorrelation and low coherence. A continuing acceleration is thus assumed. TSX provides some historic data and shows similar velocities in 2013 as in 2015. The markedly lower velocity in 2016 detected with S1 does not follow this trend; the data point may be false due to processing issues such as Fig. 10. The Saleinaz slide is visible in Google Earth (2015) and Bing Aerial (2020) orthophotos as well as SAR and optical results overlain on the Copernicus 10 m DEM. A velocity time series is extracted at the point indicated by the black star (see Fig. 11), the 2015/16 glacier extent is shown in pale blue with black outlines [2]. The color scales at the bottom indicate a velocity range of +/-30 cm/yr for TerraSAR-X (TSX) and Sentinel-1 stacking (S1) results and a velocity range of 0-8 m/yr for Sentinel-2 offset tracking (S2) results. TSX and S1 stacking results indicate a slight deceleration from 2013 to 2017 (top row). A distinct acceleration in 2018 is then noted in S1 stacking results (middle row). Later years appear decorrelated, suggesting high velocities. The S2 offset tracking results do not show movement here due to low velocities <1 m/yr that lie within the S2 error (bottom row). unwrapping errors, or else the slide experienced a pronounced deceleration in 2016. The overall velocity is similar to that of the Mer de Glace slide. Due to velocities below 1 m/yr and thus within the calculated S2 error, S2 results show no clear signal in Fig. 10.

4) Marzell Slide, Ötztal:
The Marzell slide in Ötztal is located on the northwestern flank of the glacier named Marzell Ferner and covers an approximate area of 0.166 km 2 (see Fig. 12). It appears to be part of a larger, more slowly moving complex with a total area of 0.264 km 2 . The time series in Fig. 12. The Marzell slide is clearly visible in the Google Earth orthophoto and available SAR and optical datasets overlain on the Copernicus 10 m DEM. A velocity time series is extracted at the point indicated by the black star (see Fig. 13), the 2015/16 glacier extent is shown in pale blue with black outlines [2]. The color scales at the bottom indicate a velocity range of +/-30 cm/yr for Sentinel-1 stacking (S1) results and a velocity range of 0-8 m/yr for Sentinel-2 offset tracking (S2) results. An increase in velocity in 2017-2018 can be observed in S1 stacking (top row), followed by a subsequent decrease (middle row). This trend is partially followed by S2 offset tracking results (bottom row), which show a more pronounced signal with faster movement of >1 m/yr further downslope (white diamond).   Fig. 12) for S1 stacking and multi-baseline and S2 offset tracking results. No TSX and CSK data are available for this site. The S1 errors are few millimeters to centimeters, S2 errors approx. 1 meter. We again observe a good correlation between S1 results, apart from data points in 2021. The slide appears to accelerate in 2017 and again in 2019. Maximum summer velocities remain below 50 cm/yr in S1 and 1 m/yr in S2. Due to these relatively low displacements, S2 results do not show a distinct signal at or surrounding the studied location (black star in Fig. 12). They do, however, show a more pronounced signal further downslope, particularly in 2016 and 2017 (white diamond in Fig. 12). This area corresponds to no data in S1, while the orthophoto shows fractured rocks as well as debris cover. Both indicate a possibly faster part of the observed movement with mainly superficial sliding.

VI. DISCUSSION
The results obtained with the applied methods clearly show that S1 DInSAR techniques are suitable to detect surface motion in glacial and paraglacial environments. They enable the detection and classification of >800 slope instabilities of various sizes and velocities in the three study regions. Standard DInSAR is a necessary initial method to generate differential interferograms, which are used for advanced MTI techniques. Both stacking and multi-baseline generate displacement maps and allow the extraction of velocity time series, thus yielding important information on the spatial and temporal evolution of slope instabilities. In particular, these two techniques based on S1 data enable the detection of a wide range of velocities of <10 to >30 cm/yr. As both techniques largely detect the same features and velocities, they clearly validate each other. However, each method has specific advantages. Stacking is not restricted to temporal connectivity and small spatial baselines and can be applied for seasonal and annual analyses. Multi-baseline products require temporal connectivity, yet enable an analysis of inter-seasonal displacement changes in progressive time series. The use of both methods is therefore recommended. Slower and faster movements are typically not detected with standard DInSAR, stacking and multi-baseline and require the use of alternative methods. We successfully apply PSI to detect movements of <3 cm/yr at millimeter-scale accuracy. Fast movements are best detected through use of S2 offset tracking methods. Our results therefore show that a combination of DInSAR and optical methods is necessary to attempt a comprehensive understanding of slope movement in high alpine terrain.
While S1 DInSAR techniques have demonstrated a high suitability for slope instability mapping in paraglacial environments, they are limited by different factors. The main shortcomings include the following: 1) the detection of displacement only in one dimension in the LOS; 2) a limitation to the detectable velocities due to phase aliasing effects; 3) restricted spatial and temporal coverage by reason of phase decorrelation and layover/shadowing. These limitations can be clearly observed in our obtained results, with their effect varying between the three study regions.
1) Due to the near-polar orbital flight path of most satellites, including S1, and displacement measurements only in the LOS direction, horizontal as well as movement in North-South direction is largely underestimated [57]. To derive a displacement field closer to reality, we apply a correction under the assumption of slope-parallel movement. Based on locally varying slope orientation and inclination, we calculate a correction factor to transfer LOS velocity to assumed slope-parallel movement. 2) Phase aliasing limits the maximum detectable velocity to a fourth of the wavelength λ, indicating a clear dependency on the sensor [28], [71]. For S1 differential interferograms, detectable velocities are limited to 84 cm/yr when considering 6-days repeat intervals. Movement exceeding this velocity will appear aliased in differential interferograms and cannot be classified. The existence of movements exceeding this threshold in the study sites therefore remains uncertain.
3) The spatial and temporal coverage is limited by sensor characteristics, such as spatial and temporal resolution, as well as geometric effects and surface cover of the imaged area. In mountainous regions, we generally achieve a spatial coverage on the order of 50% as a result of widespread areas of low coherence. The spatial coverage is further reduced by approximately 10% through geometric effects [28]. Rugged terrain and high topographic relief introduce layover and shadowing effects. We partly overcome this by using data from ascending and descending orbits, yet some areas of no data remain. Coherence loss occurs in all three study sites due to snow and ice cover in winter and vegetation cover at lower altitudes in summer. Therefore, we focus on available images from July to October, introducing significant data gaps during the winter months. The overall coherence varies between the study regions, despite similar topographic relief, which is also apparent in the number of detected slope movements. The Mont-Blanc massif holds the lowest number of detected features. This is because the study area is smaller compared to the Swiss and Austrian sites, but also due to low coherence throughout the region. The Mattertal site, despite being larger in area, also has a relatively low number of detected features. This is likely due to the much larger glacier coverage in the studied region. The Ötztal region in Austria, on the other hand, counts the highest number of detected features. This can be attributed to its larger area, but also overall higher coherence and less decorrelation in the S1 DInSAR results.
We overcome some of the limitations associated with S1 DIn-SAR techniques, in particular phase aliasing effects, by using alternative methods and sensors. Our analysis of four landslides occurring along retreating glaciers clearly demonstrates this. Considering results from SAR and optical sensors and methods, accelerating and decelerating phases of various velocities are determined. The observed trends indicate a relationship between increasing velocity and changing climatic conditions due to resulting effects on surface dynamics in high mountain areas [3], [6], [25], [35], [72]. This is particularly apparent for the Findel slide in the Mattertal region (see Figs. 6 and 7). We can observe a pronounced acceleration in velocity in 2013/14 relative to earlier years by comparing historic CSK and more recent S1, S2, and TSX data. This coincides with a significant retreat of the Findel Glacier in the last 20 years, suggesting increasing downslope displacement, e.g., through debuttressing effects. The velocity trend indicated by S1 data is also apparent in S2 and TSX and thus validated by these datasets. The velocity values, however, greatly differ. Maximum detectable velocities are limited to approx. 80 cm/yr in S1 data due to phase aliasing effects [71]. The consequence of this becomes apparent during the summers of 2018/19, when the slide significantly accelerated in velocity. This acceleration leads to coherence loss in S1 differential interferograms, so that the velocity can no longer be accurately detected. S2 and TSX offset tracking methods, however, are not limited by this effect and indicate much higher velocities of >7 m/yr. They are thus crucial to detect the full velocity range of the Findel slide. This example demonstrates that a combination of multiple SAR and optical sensors and methods is recommended to attempt a full comprehension of the velocity of surface instabilities with remote sensing techniques.
Results for the Mer de Glace, Saleinaz and Marzell slides equally highlight the need for a combined use of S1 DInSAR and S2 offset tracking methods. Particularly the Mer de Glace (see Figs. 8 and 9) and Marzell (see Figs. 12 and 13) slides appear to include areas that move at different velocities. For both slides, S1 DInSAR results clearly indicate a moving area of approx. 60 cm/yr, for which the corresponding time series show an acceleration in 2018/19. No strong decorrelation is observed, indicating that velocities remain within the S1 phase aliasing threshold. Correspondingly, S2 offset tracking results do not show a signal at these locations, as the velocity lies within the S2 error of 1 m/yr. For the Mer de Glace slide, however, an S2 signal is observed further North, with movement exceeding 7 m/yr (see Fig. 8). Such high velocities cannot be detected with S1 due to phase aliasing and result in areas of no data. A comparison with orthophotos shows that the S1 signal corresponds to a large rock outcrop area, covered by sparse grass vegetation and including few scree slopes. The presence of grass indicates that movements are slow enough for vegetation to grow. The S2 signal corresponds to a scree/debris surface with a very distinct upper scarp and downslope striations. Such surface features indicate rapid surficial sliding of the lateral moraine due to the glacier thinning. Whether these two movement compounds are part of a same large rock slope deformation or move with separate mechanisms cannot be determined from the available data and requires additional information from in situ field data.
For the Marzell slide, an S2 signal with velocities exceeding 3 m/yr is detected further downslope and closer to the retreating Marzell Ferner (see Fig. 12). Again, such high velocities cannot be detected with S1 DInSAR methods due to phase aliasing. Both sensors are thus necessary to identify multiple parts of the slope instability that move at different velocities. Here, a comparison with orthophotos shows a distinct upper scarp along the mountain ridge and a grassy, rocky surface, which likely defines the boundary of a larger, more slowly sliding complex that is detected with S1. Other scarps further downslope within this slide correspond to S1 signals of higher velocities, where loose blocks and debris dominate the surface. The S2 signal corresponds to a rock-and scree-covered slope with sliding and rockfall surface features that result in the observed high velocities. Based on this information, we interpret that the Marzell slide moves as one large sliding complex that includes areas of different displacement rates.

VII. CONCLUSION
The comparison of standard and advanced DInSAR methods as presented in this study successfully highlights the possibilities and limitations of detecting surface displacements in paraglacial environments of the European Alps. Based on S1 data, we mapped >800 slope instabilities in the three study regions. The resulting inventories show the wide variability in size and velocity of ground displacement observed in paraglacial environments and demonstrate the need for multiple remote sensing methods to analyze surface movement in such remote and inaccessible areas. Numerous movements of 3 to >30 cm/yr can be accurately detected with multi-looked temporal interferometry methods such as stacking and multi-baseline. We recommend the use of both techniques, as they validate each other and have specific advantages. To also detect very slow movements of <3 cm/yr, we recommend the use of S1 PSI or else high-resolution, short-wavelength sensors. The detection of fast movements exceeding >1 m/yr, on the other hand, requires the use of high-resolution TSX or optical S2 offset tracking. The effects of changing climate on surface dynamics in paraglacial environments can be observed through the analysis of four landslides that occur along retreating glaciers. The temporal and spatial evolution of these instabilities indicates a progressive increase in displacement, likely due to thinning glaciers, exposed rock and debuttressing effects. Such slope instabilities may include several, possibly independent parts that move at different rates and which can only be detected through the combined use of multiple sensors and methods. Overall, S1 and S2 data show a good correlation in terms of landforms that can be monitored and greatly enhance systematic observations and analyses of surface displacement in alpine regions. Where possible, we suggest the additional use of historic data to understand the temporal evolution of surface movement, particularly in the context of rapid glacier retreat. This is apparent for the Findel slide, where historic CSK data reveal an increase in displacement in the last 15 years.
Various further work is suggested to increase the range of detectable slope instabilities in terms of size and velocity. The use of longer wavelength satellite data such as L-Band ALOS-2 PALSAR-2 enables the detection of faster and larger motion and in vegetated areas. This can enhance the detection in paraglacial environments but also at lower altitudes, where vegetation limits the S1 suitability. Additionally, the use of S2 data could be extended to larger scales than used in this study. We present S2 offset tracking results for the four example landslides, however, coverage of the entire study sites would provide information on fast-moving slope instabilities to be included in the inventories. Furthermore, the combined use of DInSAR and offset tracking methods of various sensors can be applied to estimate local hazards. This is particularly crucial in the context of retreating glaciers that lead to glacier lake formation in high alpine terrain. Failures impeding such lakes can lead to lake damming as well as glacier lake outburst floods and are an increasing hazard in mountainous regions. The detection and inventorying of surface motion near glacier lakes using DInSAR and offset tracking methods can provide a systematic monitoring and estimation of the hazard. In such cases, and for rapid, hazardous landslides in general, in-situ data is an additional asset and source of validation (e.g., the Moosfluh sliding complex in Switzerland) [9].