The potential of point clouds for the analysis of rock kinematics in large slope instabilities: examples from the Swiss Alps: Brinzauls, Pizzo Cengalo and Spitze Stei

Lidar measurements and UAV photogrammetry provide high-resolution point clouds well suited for the investigation of slope deformations. Today, however, the information contained in these point clouds is rarely fully exploited. This study shows three examples of large-scale slope instabilities located in Switzerland, which are actively monitored for reasons of hazard prevention. We used point clouds acquired by terrestrial laser scanning to (1) identify differences in kinematic behaviour of individual rock compartments; (2) highlight active shear planes within the moving rock mass; (3) define the kinematic process driving the slope displacements; (4) model basal sliding planes based on the 3D surface movements of rock slides; (5) calculate exact displacement angles, (6) provide estimates on destabilised rock volumes. This information has significantly contributed to the process understanding and has thus supported decision-making in hazard management.


Introduction
Large slope instabilities pose a potential hazard for underlying infrastructure and settlements. Several established geodetic systems are used to monitor such instabilities. While spaceborne differential SAR interferometry (DInSAR) is frequently used to detect large slope deformations (Carlà et al. 2019), ground-based DInSAR (Gischig et al. 2009), terrestrial and airborne Lidar (Royán et al. 2014), and UAV-photogrammetry (Rossi et al. 2018) are often used to map the extent of instability more accurately. If the slope is accessible, single-point measurements are frequently carried out in addition to using a total station or in situ GNSS receivers (Limpach et al. 2016). DInSAR provides the highest accuracy of all monitoring techniques and good spatial coverage of the monitoring target. It is therefore ideally suited for early detection of small displacements in remote regions (Kenner et al. 2016). The main weakness of DIn-SAR is that only the 1D-deformation component in the radar line of sight is detected. Single point measurements are, therefore, often used to obtain information on the three-dimensional behaviour of the movement.
Lidar-or UAV-derived point clouds provide an increasingly powerful alternative. Compared to GNSS or total station data, they provide spatially distributed information on the slope instability. The accuracy of today's Lidar point clouds lies within the range of single-point methods under good measurement conditions (Friedli et al. 2019). Thus, point clouds are a very established dataset to investigate slope instabilities Rowlands et al. 2003;Teza et al. 2008). State of the art methods applied in numerous studies on landslide monitoring in the last 15 years have used point clouds to define slope deformations by by calculating one of the following: • elevation-difference models (Bitelli et al. 2004;Burns et al. 2010;Pellicani et al. 2019); • 2D horizontal displacement vector fields, based on rasterised surface structure patterns (Aryal et al. 2012;Carey et al. 2019;Kenner et al. 2014); • distances directly between the point clouds using established point-based methods such as C2C, C2M, or M3C2, which are described in detail in Method section (Fuad et al. 2018;Kromer et al. 2015).
A major disadvantage, described in more detail in Methods section, is that all these approaches capture only individual components of the 3D slope displacements. The raster-based methods capture either the horizontal component or, in the case of an elevation-difference model, the Eulerian version of the vertical component, i.e. the change in elevation at a given location in space. Point-based methods capture a mixture of vertical and horizontal components but underestimate the magnitude of the full 3D displacements clearly (Gojcic et al. 2020).
Thus, some studies deviate from the methodology to capture full 3D displacements at the expense of the spatial resolution of their results. Teza et al. (2007), for example calculated mean 3D displacement vectors for individual sub-areas of a landslide using the iterative closest point algorithm but did not take the rotational part of the movement into account. Other studies used point clouds to calculate 6-parameter transformation matrices for individual rigid body features on the landslide (Monserrat and Crosetto 2008;Oppikofer et al. 2009).
A recently developed processing method, allowing the calculation of highly resolved direct 3D deformation vectors between multitemporal point clouds (Gojcic et al. 2020), has provided new ways for more advanced usage of point clouds. This method has been successfully applied to capture the movement of a landslide (Gojcic et al. 2021), and its potential is further exploited in this study.
Absolute displacement vectors, especially in 3D, provide valuable information required for the assessment of a slope deformation. However, the key challenge for practitioners is to gain insight into the detailed compartmentalisation (i.e. boundaries of rock compartments moving homogeneously) and kinematics of rock slope instabilities, exceeding the absolute displacement rates of a point or an area. The type of movement is a central issue, as it is often not directly evident whether toppling and/or translational or rotational sliding predominates (Babiker et al. 2014;Glueer et al. 2019). Key questions concern the volume of the instability as well as the position and orientation of sliding and shear planes (Jaboyedoff et al. 2020). Furthermore, information on the structural integrity of the moving mass is important: what types of secondary deformation processes are taking place? Is the moving rock body deforming internally? Does it consist of one compartment, or can individual compartments be distinguished, and how do they move relative to one another? (Agliardi et al. 2001;Donati et al. 2019;Gischig et al. 2011;Stead and Wolter 2015).
Previous studies on point-cloud-based landslide analysis struggled to provide answers to these questions, mainly because smallscale deformation signals were covered by the much stronger largescale displacements or because the relevant movement components were not captured. In this study, we provide an overview of different data processing strategies which can significantly increase the information content of a point-cloud analysis. To this end, we present three examples of multi-million cubic meter rock slope deformations in the Swiss Alps with different characteristics. Terrestrial Lidar measurements have been carried out on all of them at regular intervals for more than one year. Tailored processing methods have been applied to each dataset, enabling the creation of novel products such as: I. spatially highly resolved maps of the displacement dip angle; II. 3D models of sliding planes derived from 3D surface displacements of rock slides; III. relative deformation maps show the relative reshaping of the rock mass during the absolute displacement process. The movement of individual rock compartments relative to each other gives valuable indications on the type of the overall moving process.
The results of our analyses are discussed from a geological perspective to highlight their relevance for each case.

Landslide Brienz/Brinzauls
The rock slope above the village Brienz/Brinzauls (Grisons) moves at a rate of several m/year towards the village. The village itself lies on a terrace that is currently moving at about 1.5 m/ year towards the valley incised by the river Albula. The instability above (referred to as "Landslide Berg") and below the village (referred to as "Landslide Dorf ") are not distinct landslides, but part of a large landslide complex that became unstable in prehistoric times and now involves more than 166 million m 3 of rock over an area of 2.4 km 2 . Sliding has been registered since the nineteenth century since a considerable rockslide detached from the eastern part above the village in 1878 and crept downslope by a few meters per day ("Igl Rutsch," document in Heim (1881)). For decades, the landslide moved at a rate of a few centimetres per year until acceleration began in 2000-2010, and increased rockfall activity triggered investigations of the slope above the village. Since 2011, a monitoring system comprising automatic total station measurements and periodic GPS measurements has served as an early warning system for potential rockslides threatening the village. The measurements showed dramatic accelerations that temporarily ceased between 2015 and 2017, following two particular dry winters, but have accelerated all the more since then. The velocity of the sliding mass beneath the village developed synchronously, although at a lower level, and associated ground deformation led to severe damage to houses and infrastructure. Since summer 2019, the monitoring system has been complemented with an automatic terrestrial DInSAR system and 11 automatic GNSS stations. In 2018, an exhaustive investigation was initiated to unravel the geological, hydrogeological, geomechanical, and kinematic conditions of both landslide sectors "Berg" and "Dorf," allowing to develop refined future failure scenarios and to assess possible mitigation measures (BTG Büro für Technische Geologie AG, 2022). The terrestrial laser scanning (TLS) data presented here are part of this ongoing monitoring program.
In the geological context, the landslide complex is located at the boundary of the North Penninic and the overlying austroalpine nappes. The landslide "Dorf " is formed predominantly by North Penninic turbidites ("Flysch"). In the landslide "Berg" nordpenninic Flysch is overlain by the Allgäu formation consisting of interbedded limestones and schists (Austroalpine), overlain by Raibler layers (yellowish band of rauwacke and dolomites), and Arlberg dolomite on top. The lithological contacts are not sedimentary but tectonic contacts that likely dip gently to moderately into the slope.

Rock slope failure Pizzo Cengalo
The Pizzo Cengalo (3369 m a.s.l.) is situated in the southern end of the Bondasca valley in the Bergell region in Southeastern Switzerland at the border to Italy. Together with the Pizzo Badile to its West, it represents the western end of the tertiary magmatic intrusion called the Bergell Pluton complex (Fig. 1).
In contrast to the compact Pizzo Badile, Pizzo Cengalo is characterised by several fault zones and more intense fracture cleavage. In the climbers' community, the northeast face of Pizzo Cengalo has been known to have a higher rockfall activity than the nearby Pizzo Badile (Miotti and Gogna 1985). The northeast slope of Pizzo Cengalo forms the western flank of a small and narrow, north-facing glacier cirque. The glacier at the foot of the Cengalo northeast slope is mainly fed by snow avalanches from the steep rocky slopes of this cirque.
The first large rockslide of this century occurred in July/August 2003 at Pizzo Cengalo, when a part of the "Kaspar Pillar" collapsed, which represents the northern boundary of the northeast slope ( Fig. 1). Since then, a series of rockfalls and bigger rockslides originated from the northeast slope (Table 1).
In December 2011, a major rockslide of about 1.5 million m 3 detached from the northeast flank of Pizzo Cengalo between approximately 3100 and 3000 m a.s.l. and developed into a rock avalanche that stopped in the valley bottom of Val Bondasca. Large parts of the detachment scar were covered with cleft ice, indicating the presence of permafrost. Prior to the event, increased rockfall activity and wide-open tension cracks had initiated comprehensive investigations by geologists. In 2013, the joint system was analysed in-depth, and a kinematic model was established (De Preux 2014). After the 2011 event, a comprehensive monitoring strategy was established, relying mainly on a portable radar interferometer (deployed once or twice a year), a terrestrial laser scanner (Terrestrial laser scanning section), and photo analysis (annual survey by helicopter). The monitoring indicated that large parts of the northeast slope were unstable. From 2012 to 2017, an increase in displacement velocity of the unstable rock mass was measured. Though it became clear that a collapse of the instability was imminent, an exact determination of the timing was not possible. On August 23, 2017, almost the entire unstable rock mass of about 3 million m 3 failed, causing a rock avalanche of 3.2 km in length, which led to human casualties on a hiking trail. The entrainment of glacier ice as well as entrainment and overloading of saturated, unconsolidated deposits of the former rock avalanche deposits led to the immediate formation of debris flows, which reached settled areas around Bondo and caused significant damage to infrastructure (Walter et al. 2020).

Rock slope instability Spitze Stei
The Spitze Stei rockslide above Lake Oeschinen in the Kandersteg region (Bernese Alps) comprises an area of ~ 0.5 km 2 , ranging in elevation between 2200 and 2900 m a.s.l. The slide involves approximately 20 million m 3 of rock material across two distinct rock formations within the Doldenhorn Nappe: the Öhrli formation, consisting of limestone and marl above approximately 2600 m a.s.l. and the Zementstein formation, consisting of marly and clayey schist below 2600 m a.s.l. The layers dip at an angle of approximately 30° towards NW. In the upper portion of the rockslide, bedrock is exposed at the surface; in the lower portion, the surface is covered by scree.
The geologic setting contributes to an unstable pre-disposition, which has led to large rock avalanches in the past. A rock avalanche with a volume of roughly 800 million m 3 occurred to the West of the current instability (Fisistock Bergsturz) approximately 3200 years BP (Singeisen et al. 2020;Tinner et al. 2005). The area below the current instability was subject to a large rock avalanche of approximately 2300 years BP (Oeschinensee Bergsturz (Köpfli et al. 2018)). This rock avalanche dammed the current Oeschinen lake.
The uppermost portions of the slide area (above 2800 m a.s.l.) lie within the permafrost, as suggested by modeling (Kenner et al. 2019) and confirmed by in situ electrical resistivity tomography surveys (Hilbich and Hauck 2019), borehole measurements (Geotest 2020a, b, c), and direct observations of freshly exposed permafrost ice (Geotest 2019). Below 2800 m a.s.l., there is a transition zone in which permafrost is already thawed in the upper 20 m, but still persists at larger depth. Some ice-rich permafrost exists, most notably in the Eastern portion of the slide area, where there is a small rock glacier. While the eastern part of the slide has moved for decades, according to orthoimage analysis (SLF, unpublished), other portions of the Spitze Stei instability began moving in the early 2000s, as shown by retrospective analyses of space-borne radar data via InSAR analyses (Caduff et al. 2021) and airborne orthoimages via feature tracking (Geotest 2020a, b, c). Acceleration over the past decade has been strong and nonlinear and accompanied by an expansion of the size of the unstable perimeter (Caduff et al. 2021;Geotest 2020a, b, c). An increasing number of growing fractures in the exposed bedrock, increasingly tilted blocks, and more frequent rockfall events from the Spitze Stei area have been clear signs of the progressing destabilisation. Today, displacement rates are on the order of several meters per year. Despite these high displacement rates, only minor rockfall events (volumes in the order of several 1000 m 3 ) have occurred so far.

Terrestrial laser scanning
Terrestrial laser scanning (TLS) was carried out at all three sites using a Riegl VZ6000 long-range laser scanner. Details regarding measurements dates, ranges, and resolutions are summarised in Table 2. Due to the long measurement ranges, the scans were carried out before sunrise at "Spitze Stei" to avoid measurement distortions caused by refraction (Friedli et al. 2019). The precision of the multitemporal scans ranges between 1 and 2 cm for Brienz, 2 and 3 cm for Pizzo Cengalo, and about 5 cm for Spitze Stei. These values refer to the ability of the method to reproduce unchanged terrain in multitemporal measurements. They, therefore, represent the accuracy or the level of detection for recorded changes. The values do not refer to the accuracy of a single point (i.e. small-scale measurement noise), instead, they represent the maximal size of large-scale distortions extending over square decametres, which corresponds to the scales of the slope deformations discussed.

Co-registration strategies
Co-registration of the multitemporal point-clouds was based on the internal geometry of the point clouds. Several commercial software implementations for this type of point cloud matching exist. Most of them represent variations of the iterative closest point (ICP) algorithm (Chen and Medioni 1991), which minimises the sum of point distances between two or more point clouds.
During the iterative process, areas showing particular high deviations between the point clouds, e.g. due to deformations of the scan object, are excluded prior to the following iteration. In this study, an algorithm called "multi-station adjustment," provided by the Riegl software RiSCAN Pro was used.
Depending on the geological question being investigated, ICP matching can be used in at least three different ways (Fig. 2). The most common approach is to define reference areas in stable terrain and to use a six-parameter transformation to reposition and reorient the scans. The algorithm will then match the repeat measurement to the reference measurement, causing stable areas to be congruent and deformation areas to show the absolute deformation rates between the point clouds (Co-registration A).
However, absolute large-scale displacements across a rock slope are often so strong that internal small-scale deformations within the moving rock mass cannot be detected. It is therefore advised to use the deforming area as a reference for the ICPmatching instead of stable terrain. Given that the analysed deformation area has roughly the same direction of movement, this approach causes the algorithm to compensate the mean absolute displacement of the rock mass, which then makes the internal deformation of the unstable rock mass visible (co-registration B). In simple words, co-registration method B pushes the rock instability back into its original position. It is then possible to compare the shape of the displaced rock mass with its original shape prior to the displacement. Thus, it becomes evident how the rock mass deformed internally during the process of displacement. This, in turn, provides indications of what type of displacement process took place, for example, a toppling or sliding process.
Depending on the type of movement, it may be useful to only compensate for the translation component of the movement. This will provide a visualisation of the rotational component and can, for example, help to distinguish tilting from sliding (co-registration C).

Deformation analyses
The most common and simplest way to compare two point clouds of the ground surface is to generate digital elevation models (DEM) as horizontal surface rasters and to subtract them from each other to obtain a digital elevation-difference model (DEDM). In addition to its simplicity, this approach easily allows the calculation of volumes, for example, of the detached mass after a rock slope failure.
However, the deformation analysis is clearly limited within a DEDM. Although it provides area-wide information, it only displays changes in elevation at certain points in space, delivering what might be described as a 1.5D product. Although horizontal displacements are not considered in a DEDM, they can cause strong elevation changes in the DEDM due to a shift of topography. Particularly in very steep terrain, this might cause unusable results, and the slope dependency of the elevation changes can lead to misinterpretations.
To avoid these problems, an alternative way to derive terrain deformation from point clouds is to compare the point-clouds directly, either by (1) computing the distances between closest points of two clouds (cloud to cloud, C2C), or (2) by considering a set of points to compute the local distance between two point-clouds (multiscale modelto-model cloud comparison, M3C2), or (3) by calculating the shortest distance between a point and the closest triangle of the reference mesh (cloud to mesh, C2M) (Cignoni et al. 1998). Compared to a DEDM, all the methods are less dependent on the slope angle and deformations are easier to interpret and to compare within the deformation model (Girardeau-Montaut et al. 2005). C2m and M3C2 are the most commonly applied methods. They generally provide similar results, however, dependent on the monitoring question and the surface roughness, C2M or M3C2 might be more appropriate (Barnhart and Crosby 2013). We used C2M in combination with co-registration method B (relative deformations) to analyse the change in the shape of a rock mass during its absolute displacement within the rock slope. However, when it comes to absolute displacements, none of the aforementioned approaches provide clear evidence on the amount and the direction of the movement. Instead, they show a mix of various effects such as: a) a slope and an aspect-dependent 1D component of the absolute movement; b) a small-scale deformation pattern caused by shifts of the smallscale surface structures (particularly if the movement direction is parallel to the surface); c) a potential internal deformation of the moving rock mass.
To determine realistic displacement rates, it is necessary to compare individual objects or structures in two point clouds. A first, relatively simple, and rapid approach is the calculation of horizontal 2D deformation vector fields. This approach is based on DEMs, which were processed with a high pass filter of a certain size. This filter eliminates the large-scale topography but preserves-depending on the filter size-the small-scale surface structure (Kenner et al. 2014). This surface structure provides a unique pattern across the entire DEM, which can be correlated among multitemporal DEMs unless the small-scale surface structure has changed too drastically since the reference measurement was carried out. The correlation window can be shifted over the entire high pass filtered DEM, providing 2D displacement vectors for each position. Here we use the particle imaging velocimetry correlation method introduced by Roesgen and Totaro (1995). The resulting displacement vector field (DVF) represents the absolute horizontal component of a slope movement. Faulty correlations are automatically filtered from the DVF based on the direction and magnitude of neighbouring vectors (Kenner et al. 2014). In the case of a very steep rock slope, the DEMs can be created in an elevation projection, causing the DVF to provide displacements in the XZ-or YZ-plane of a local coordinate system.
The horizontal 2D DVF has clear disadvantages in comparison to 3D DVF, especially if the vertical component of the displacement is large and/or varies within the slope. A large vertical component causes a strong underestimation of the actual displacements and variations of the latter to prevent the comparability of the displacement vectors within the 2D DVF. Gojcic et al. (2020) presented a new framework for point-cloudbased deformation analysis denoted as a feature-to-feature supervoxel-based spatial smoothing (F2S3), which avoids these problems by directly determining the 3D DVFs. F2S3 can be split into two steps, both inspired by recent advances in machine learning. In the first step, the initial 3D displacement vectors are established by computing the corresponding points across the epochs. These correspondences are established by comparing their local feature descriptors, i.e. high dimensional vectors describing the local geometric properties in the vicinity of each point (Gojcic et al. 2019). Due to the sampling pattern, motion, occlusions, and limited descriptiveness of the local feature descriptors, the initial correspondences include outliers. The second step of F2S3 is therefore devoted to filtering out the initial displacement vector field, based on the local rigidity assumption, i.e. local patches (supervoxels) are assumed to move as rigid bodies (Gojcic et al. 2020). The output of F2S3 is a dense displacement vector field comprising a 3D displacement vector for each point for which a reliable correspondence was established. For downstream processing, this dense displacement vector field can be summarised in the form of a regular grid of 3D displacement vectors, where each grid point is represented by a median 3D displacement vector of that cell.
Detailed deformation analysis of a rock slope generally requires a combination of at least two of the methods presented above. Each of these approaches has advantages and disadvantages but is most powerful to address different geological questions if used in combination.

Integrating a sliding plane from 3D deformation vector fields
The local surface movement direction of a pure rock slide (i.e. sliding on a single connected sliding plane as the only kinematic mode) is mainly controlled by the shape of the sliding plane. For example, a convex sliding plane will cause steeper 3D deformation vectors in the upper parts of the rock slide and increasingly flatter vectors in the lower parts. This relation was used by Carter and Bentley (1985) to propose a method to interpolate a profile of the sliding plane from 3D surface point measurements. Later on, this method was refined by Baum et al. (1998), taking into account the internal deformation of the sliding rock mass. While these methods worked with singlepoint measurements, the 3D DVFs presented in our study dramatically increase the information contained on the shape of sliding planes. For the first time, information on the orientation of the sliding plane exists area covering in a very high resolution, providing a high redundancy, consequently smaller influences of local disturbances due to shallow rock movements and thus higher robustness and accuracy.
To integrate a sliding plane from the 3D DVF, we calculated local planes, each defined by a single 3D deformation vector-representing the local dip vector-and the corresponding strike vector. These local planes can be represented by plane normals. The problem of shape reconstruction from local planes or surface normals arises in many computer vision applications, e.g. shape-from-shading and photometric stereo (Horn and Brooks 1986). There, the surface normals are turned into an estimate of the gradient of a 2D depth map encoding the surface geometry. The challenge is then to integrate this gradient field into the underlying depth map by adjusting the residuals caused by the high level of redundancy in this dataset.
A variety of numerical solutions have been proposed for this task (Quéau et al. 2018a). They mainly differ in their ability to handle: (a) noisy normal measurements; (b) a non-rectangular 2D integration domain; (c) a free boundary condition; (d) possible depth discontinuities. Variational methods offer a good compromise for this task: they minimise the discrepancy between the given gradient and that of the estimated depth. Hence, they are well suited for handling noisy 3D DVF measurements. Besides, they handle non-rectangular 2D grids and free boundaries, as is the case for the measured 3D DVF considered here. Since sparse vertical steps might occur in the sliding plane of the rock slide, the depth should be considered. For these reasons, we considered the total-variation-based integration method proposed in (Quéau et al. 2018b). It minimises the discrepancy between the measured gradient and that of the underlying depth map in the sense of the L1 norm (sum of absolute differences), which allows sparse depth discontinuities over the domain where data is available.
The output of the method is a 2.5D representation of the surface geometry in the form of a depth map. The resulting sliding plane is horizontally georeferenced, but z coordinates are in a local system since gradient integration is only possible up to an integration constant. The vertical offset of the sliding plane must therefore be corrected by aligning the modelled sliding plane to a known 3D point of the actual sliding plane, for example, to an outcrop area or a borehole lowered into the sliding plane. Intersecting the georeferenced sliding plane with the DEM ultimately reveals the thickness of the sliding mass.

Brienz-a sliding and toppling rock mass consisting of a complex set of interacting compartments
The unstable rock mass above the village Brienz (landslide "Berg") poses a substantial risk for the settlements and infrastructure of Brienz and Vazerol. The TLS perimeter focuses on the South oriented part of the Landslide "Berg" above Brienz (Fig. 1). We used coregistration method A in stable terrain to the East of the landslide and 3D displacement vector fields from F2S3 to quantify the absolute movement of the rock instability between August 5, 2020, and November 30, 2020 (Fig. 3). A distinct increase in 3D displacement Fig. 3 Magnitude of the absolute 3D slope deformation (mm/d) in the South oriented part of the rock slope instability above Brienz between August 5, 2020, and November 30, 2020. Yellow points mark boreholes drilled in the moving mass. Relief from airborne Lidar data acquired in 2019 rates is evident between the lower half of the measurement perimeter (~ 6 mm/day) and a belt located between approximately 1450 and 1700 m a.s.l. This rapid belt includes a zone directly below the "failure edge" (~ 13-23 mm/day), as well as a compartment towards the western "Catgeras" ridge. Within the rapid zone, the velocity increases from the base upwards. Velocities exceeding 17 mm/day are limited to two individual areas on the upper boundary of this zone. The high velocities, the strong disruption of the rock mass as well as the fast reaction of velocities to water infiltration (as seen from total station measurements) indicate that the rapid zone represents shallow movements of the rock mass superimposed on the slower and deeper-seated instability. At the "failure edge" and the "plateau" above, deformation rates abruptly decrease again to about 8 mm/day.
The absolute movement rates alone do not fully reveal the detailed rock kinematic processes taking place because the relatively high velocities conceal internal spatial velocity variations, i.e. deformations. Therefore, we used co-registration method B to compensate for the main movement of the landslide. The internal deformation pattern of the rock mass during its downslope movement thus became visible. For the visualisation, we used C2M and a 2D displacement vector field.
C2M revealed five distinct, parallel, and horizontal lineaments within the lower half of the moving mass (Fig. 4a, black arrows). These lineaments coincide with daylighting large-scale discontinuities (i.e. faults or fracture zones) that dip moderately into the slope.
The relative deformations at these lineaments indicate surfaceorthogonal shearing along the discontinuities of a few to several decimetres per year, causing relative heave of the surface below the rupture line and relative subsidence of the surface above the rupture line. The upper half of the sliding mass, which is delimited downslope by lineament 5, is concealed by relatively strong smallerscale deformations in the surroundings, which are hard to discern using C2M.
However, the 2D displacement vector field provides clearer evidence on the relative deformations of the sliding mass in the upper half and additional information on the lower half (Fig. 4b). The upper half can be subdivided into at least 9 individual compartments, which move significantly relative to each other. The boundaries of these compartments are defined by distinct changes in magnitude and/or direction of the relative deformation vectors over small areas. In most cases, they coincide with lineaments created from geomorphological and geological mapping (white lines in Fig. 4b), especially in the rapid zone they give additional evidence.
In the lower half of the unstable rock mass, at least one additional lineament became evident in the 2D deformation vector field (no. 4b). In contrast to the five lines found in the C2M visualisation (Fig. 4a), rupture line number 4b is characterised by apparent surface parallel shearing and compression and was therefore not discernible in the C2M visualisation. Geomorphological evidence confirms the existence of the rupture lines: Fig. 5 shows distinct changes in grain size and slope angle below and above the rupture Although the elevation angle of the surface movement (a) widely corresponds to the slope of the modelled sliding plane (b), this plane crops out at the black line in part b. However, the true sliding plane likely lies at more than 100 m depth here, as indicated by borehole deformation measurements. The black arrows mark ripples in the modelled plane, which indicate tensions during the process of normal integration and might indicate the location of plane discontinuities or deformation process changes. Relief from airborne Lidar data acquired in 2019 lines in talus-covered areas. In bedrock areas, the rupture lines correspond to continuous rock bands forming flat ledges in the rock slope. At this point, the interpretation of these rupture lines was ambiguous, they could either be a sign of compressive sliding on a concave surface or of toppling. In fact, a borehole east of these ripples (KB11, see Fig. 3) confirmed toppling at least on the "Caltgeras" ridge ( Figure S1).
To generate more evidence for the moving areas in the East of the Caltgeras ridge, we tested the hypothesis that sliding predominates here. To do this, we integrated a siding plane from the 3D vectors for this area. For vertical georeferencing, the sliding plane was attached to borehole KB10 (Fig. 3 and Figure S1), in which the sliding plane was found at around 48 m depth. The spatial pattern of the slope of the sliding plane and of the elevation angle of the 3D surface movement correspond well, suggesting a successful integration of the normal vectors during the calculation of the sliding plane (Fig. 6). Nevertheless, the modelled sliding plane crops out at 1250 m a.s.l., which is about 100 m above the village (Fig. 6). However, according to borehole deformation data, the actual sliding plane is at around 150 m depth in borehole KB1 next to the village (Figs. 3 and S1). Thus, it is not possible to derive a plausible sliding plane from the surface movements, implying there must be either a different deformation process taking place (e.g. toppling instead of sliding) or that the sliding plane is not continuous but interrupted by vertical steps or superimposed by a secondary sliding plane with strongly diverging orientation. Small ripples at the eastern and western borders of the modelled sliding plane might indicate the locations of such discontinuities or areas where predominant sliding gives way to a different kinematic mode such as toppling. As the prominent rupture lines described above might be a sign of toppling as well, we consider this as the most likely movement process here. Our methods allowed a detailed partition of the Brienz rock slope instability into individually moving rock compartments, which were not discernible from absolute deformation data. Their distinction and the analysis of their translation relative to each other allowed conclusions on the type of deformation process. By integrating a potential sliding plane from the surface displacements, we could exclude sliding on a continuous plane as the sole deformation process. In fact, both sliding and toppling are likely to occur in concert within the unstable rock mass.

Pizzo Cengalo-toppling of five granite compartments
The absolute displacements of the Pizzo Cengalo NE-face between July 12, 2013, and September 13, 2016, were quantified using the coregistration methods A and F2S3 (Fig. 7). This analysis revealed a 3D deformation area of about 135,000 m 2 . The displacements increased gradually from zero at the base of the deformation area to about 60 cm at the top. A useful point of orientation in the following figures is an about 170 m high rock pillar called the "Kasper" pillar in the Northwest of the instability (Fig. 1).
The displacement magnitude alone provided little information on the detailed deformation process. The 3D DVF from F2S3 allows, however, the calculation of any movement-and direction components such as the elevation angle (dip) or azimuth (dip direction) of the displacements. The elevation angle mapped in Fig. 9 shows a striking step-like pattern, where the steps occur exactly at the outcrop lines of the joint set 3. Below these outcrop lines, i.e., at the upper edge of a joint-bordered rock body, the elevation angle is smaller and increases sharply above the outcrop lines, i.e. at the lower edge of a joint-bordered rock body. We interpret this as a toppling movement along with set 3. Set 1 delimits the toppling instability laterally, set 2 towards the back. The inset in Fig. 9 explains the decrease of the elevation angle of the movement from bottom to top of a toppling rock compartment.
We additionally applied co-registration methods B and C in combination with a C2M visualisation. Co-registration method B confirmed the results suggested by the analysis of the elevation Fig. 9 Elevation angle of the slope movement in the Pizzo Cengalo northeast face (from F2S3). The red, blue, and black lines mark the three joint sets, corresponding to Fig. 8. Caused by their different position relative to the centre of rotation, the elevation angle of the movement is larger at the lower edges of a toppling rock compared to the upper edges (see inset). The figure is not true to scale, the height of the "Kasper" rock pillar is given for reference Fig. 10 Magnitude of the absolute 3D slope deformation of the rock slope instability Spitze Stei between July 19, 2020, and July 17, 2020. The instability is divided into 6 sectors, showing differing deformation rates and directions. The yellow points mark the location of three boreholes, which TB4 ( Figure S3) was used to validate our model results. Orthophoto © Geotest AG angle. Clear signs of shearing appeared along the outcrop lines of Set 3 ( Figure S1). The shearing became apparent by offsets of a few cm along the joints in form of a relative heave of the rock compartment below the daylighting discontinuity and relative subsidence of the rock compartment above, similar to the pattern observed at the Brienz rock slide ( Figure S2).
The mean deformation magnitude of the instability (irrespective of arithmetic sign) from C2M was 12.4 cm for co-registration A (absolute values), 10.3 cm for co-registration B (compensation of translation and rotation), and 12.1 cm for co-registration C (compensation of translation, no compensation of rotation). This means co-registration B reduced the mean deformation magnitude by 21 mm, while co-registration C caused a reduction of 3 mm only. However, the compensation of translations should strongly reduce any displacements caused by sliding processes. Since no such reduction occurred, we do not expect any sliding component, and the entire movement of the rock slope can be explained by toppling, i.e. rotation of rock compartments separated by set 3 discontinuities. Our results correspond well with the structural analysis carried out prior to the rock slope failure (De Preux 2014).

Spitze Stei-where are the sliding planes?
At least 6 sharply defined sectors with clearly differing displacement rates can be distinguished at the rock and debris slide Spitze Stei (Fig. 10). Over a period of 60 days, between July 19, 2020, and September 17, 2020, sector 5, which forms the uppermost portion of the slide, showed 3D displacements of about 50 to 80 cm. Sector 6, a subsection of the rocky summit portion moved slightly faster. The debris-covered sector 4 at the western base of sector 5 moved up to 2.5 m over the same time period. A sharp increase in movement is visible between sectors 5 and 4, causing the development of an extension zone between the foot of the rock portion and the debris compartment. The same applies to the border between sectors 1 and 4. The higher displacement rates in sector 4 are likely caused by shear and sliding processes occurring in multiple layers superimposed on a deeper-lying basal sliding zone. Sector 2 showed displacement rates of around 25 cm during this period. Sector 3 represents a transitional zone between sectors 2 and 4, with displacement rates increasing towards sector 4.
Below the border between sectors 5 and 3, a stable rock patch exists in sector 3, sharply separated from movement above in sector 5. We inferred from analysis of the deformation field, point-cloud structure, and field photos that the sliding plane labelled A daylights here (Figs. 11 and 12). This sliding plane A dips northwestwards and constitutes the basal sliding plane of sectors 1, 4, 5, and 6 ( Fig. 11). Although the rock patch below the outcrop line (coloured formation in Fig. 12) is stable (confirmed by measurements of two surveying prisms installed in 2020), the talus slopes around and below are moving. The corresponding movement is initiated by the motion of sectors 2 and 3. Sector 2 must be sliding on an independent, deeper-lying plane B. The spatially homogeneous deformation rates point towards a deep-seated rock slide. Outcrop locations of the siding plane B of sector 2 are marked in Fig. 12. Sector 3 constitutes a transition area between the rock masses sliding on planes A and B.
Based on the analysis of historical orthophotos, sectors 2 and 3 have already been active since the 1970s, while the movement of sectors 1 and 4 started during the 1990s, and the deformation onset of sectors 5 and 6 was after 2010.
Sliding planes were integrated from the 3D surface displacements for sectors 2 and 5. In sector 1, large data gaps prevented the integration of a sliding plane. Indicated by strong changes in movement direction and magnitude, the basal sliding plane in sectors 3 and 4 are superimposed by the motion of loose rock material, whose magnitude and direction changes during the course of the year, indicating the existence of multiple shear layers with varying activity. As our method cannot distinguish between multiple movement processes, it cannot be applied here. Sector 6 is assumed to show toppling in addition to the gliding process.
Plausible results were obtained for the integration of sliding plane B in sector 2, which is characterised by spatially homogeneous surface movements with a relatively small elevation angle compared to the slope angle of the terrain. Sliding on a continuous plane may best explain the motion pattern observed. The slope angles of the modelled sliding plane correspond well and area-wide with the elevation angles of the 3D surface movements. This indicates that a single sliding plane can explain the surface movement well. Only small tensions occurred in the modelled sliding plane during the process of normal integration. We corrected the vertical offset of the integrated sliding plane by attaching it to the outcrop lines at the base of the slope, indicated by bulgy structures and the deformation vector field. Though an overall validation of our result is not possible, the 41.5 m deep borehole TB4 (Figs. 10 and S3), drilled in the centre of sector 2, does not contradict the modelled sliding plane: Like the other two boreholes shown in Fig. 10, TB4, unfortunately, did not reach the sliding plane. This agrees with the 51 m depth of the modelled plane at the location of the borehole. Moreover, the modelled outcrop line at the base of the slide fits well with the course of the observed outcrop line, indicating a realistic model result.
In the lower third of sector 2, the terrain has a slope of approximately 35°. Here, the thickness of the sliding mass increases steadily from bottom to top (Fig. 13). Above, the terrain flattens to about 30°, and the thickness of the sliding mass remains constant at about 45 to 50 m. A rocky ridge forms the uppermost part of the rock slide. The great thickness of the sliding mass at the upper boundaries of the instability indicates an almost vertical backward detachment surface of the rock slide. This observation corresponds to former rock slope failures in the close vicinity, which showed sliding planes of 30 to 35° steepness and subvertical backward detachment surfaces (Fig. 1). The thickness of the sliding mass under the ridge, exceeding 70 m locally in Fig. 13, might be overestimated because the sliding plane algorithm was not able to model the shape and exact location of the subvertical backward detachment planes correctly. Considering this uncertainty, the total volume of this sector is estimated to lie between 2.75 and 2.95 million m 3 .
The motion of sector 5 is difficult to describe with sliding alone. The slope angles of the integrated sliding plane were not coherent with the elevation angles of the surface movements, nor did the modelled and observed outcrop lines correspond very well. This is likely caused by an almost elevation-parallel discontinuity at about 2660 m a.s.l., which divides sector 5 into a lower part (5a) and an upper part (5b). In the eastern part of sector 5, this discontinuity causes a distinct rift at the surface (III in Fig. 14b). While sector 5a slides at approximately 30° on plane A according to the elevation angle of the movement (Fig. 14a), the trapezoidal or wedge-profiled sector 5b subsides into an extension zone formed between the slide movement of 5a and a steep backward detachment surface (Fig. 14). This results in a sharp increase of the elevation angle of the movement to more than 50° in 5b (Fig. 14a). In addition to settling, the eastern parts show a toppling towards the East, deviating the total displacement direction from NW to NE. These observations imply that the rock mass does not slide on a continuous rotational sliding surface but on a biplanar surface consisting of a sliding surface and a steep detachment zone.
There is a remarkable bulge at the base of sector 5 (Fig. 11a). The gain in mass along this elevational belt is visible in the C2M visualisation for the absolute displacements in Fig. 11a. In the eastern part of sector 5, the strong bulging is obviously caused by the sliding process of sector 5 along sliding plane A which daylights here (Fig. 11). The bulging belt continues in the same size and magnitude westwards in the debris material, suggesting that sliding plane A continues in the same direction along the border of sectors 5 and 4 (Fig. 11a) and crops out under a layer of talus here. This observation excludes the possibility of another rock compartment sliding on plane A below this outcrop line in sector 4. Indeed, the elevation angles of the movement in sector 4 imply that the thickness of sector 4 might not exceed 20 m to 25 m in its eastern part (4a) but are much larger in its western part (4b). This in turn indicates that plane A might not be overlain by another rock compartment in sector 4a but only be covered by a layer of rock debris here (Figs. 11a and 14). It would furthermore imply that the eastern part of sector 5 is-aside from the debris burden-kinematically free at its base.
The deformation processes and the interactions of the individual rock compartments of Spitze Stei are complex. Nevertheless, our results provide useful indications that are currently integrated into a broader geological model: Outcrop locations of the basal sliding plane A are directly evident between sectors 5 and 3 and indicated by a bulging of the debris slope between sectors 5 and 4. The eastern part of Sector 4 seems to consist only of loose material which slides on plane A and is pushed over the terrain edge at the lower boundary of sector 4. Sector 5 thus gains increasing kinematical freedom at its base. In sector 2, a pure sliding process, affecting a volume of almost 3 million m 3 , appears most likely. We see clear indications for steep backward detachment surfaces of both rock slides in sector 2 and 5. The steep detachment surface causes subsidence of the uppermost portions of sector 5 into an extension zone.

Three-dimensional deformation vector fields and angles of movement
We used 3D DVF for the analysis of all three slope instabilities. The accuracy of these vector fields was analysed by Gojcic et al. (2021), using the Brienz landslide as an example. The comparison with other measurement methods generally reveals small deviations which are not due to measurement or method errors but are caused by the fact that the methods record other movement signals. Nevertheless, the displacement vectors of F2S3 showed deviations to reference point measurements of a few cm only.
The displacement angles, especially the elevation-or dip-angles were of particular value at all sites. At Pizzo Cengalo they provided evidence of a toppling process in accordance with an independently prepared kinematic analysis. At Spitze Stei they allowed for manual expert profiling of the instability in locations where automatic sliding surface modelling was not possible. Although not shown here, the elevation angles confirmed as well our findings at the Brienz landslide. As the interpretation of these angles is qualitative, there is no validation in terms of independent measurements. However, pronounced changes in 3D movement direction occurred at all three sites exactly along prominent geomorphological lineaments. We consider this as proof of the reliability of this analysis.

Sliding plane models
Detailed information on the shape and location of sliding planes provides insights undoubtedly very helpful for estimating, for example the volume and criticality of slope instability. The idea of sliding plane modelling from surface displacements is over 35 years old (Carter and Bentley 1985). The only obstacle to its widespread application was the lack of area covering 3D displacement information. This information is now provided by F2S3. Due to the highly redundant dataset, the method became much more complex in comparison to the first approach by Carter and Bentley (1985). The basic principle, however, remained the same: it is based on the hypothesis that the differential directions of the surface displacements of a pure rock slide reflect the orientation of the sliding plane at the corresponding location. Considering this hypothesis as true, our method would do a perfect job. In those parts of the rock instabilities in which sliding can be assumed as a clearly dominant process, the dip angle of the surface displacements corresponds very well with the slope angle of the modelled sliding plane. In the case of sector 2 of the Spitze Stei instability, for example, the RMSE between both angles was less than 3°.
There are uncertainties that cause deviations between the modelled and the true sliding plane. Often, there are secondary movement processes such as subsidence, toppling or sliding plane parallel rotations that superimpose the movement signal of the sliding process. Although vertical rotations, induced by the shape of the sliding plane, are captured by our method, such rotational slides might cause internal deformations of the sliding rock mass, which influence the surface displacements and cause additional deviations to the sliding plane model (Baum et al. 1998).   Fig. 11b. The rock formations coloured in green, magenta, and blue are stable and located directly below the basal gliding plane A (sectors 1;4:5). The dashed red line marks the approximate outcrop location of gliding plane B of sector 2 in debriscovered terrain (deduced from deformation data). Red arrows mark locations along the outcrop where bedrock is exposed. The dashed orange line delimits the main deformation area towards the east. The dashed white line shows the deformation boundary between the two basal gliding planes A and B. The yellow points mark surveying prisms installed in the slope and numbered to facilitate orientation. The distance between prism 4 and 5 is about 160 m Fig. 13 Thickness of the sliding rock mass in sector 2. The red, white, and orange lines correspond to those in Fig. 12. The black dashed line marks the intersection of the modelled sliding plane with the DEM. Borehole "TB4" is marked with a yellow point (compare Fig. 10). Orthophoto © Geotest AG Since all these distortions cause an increase of the before mentioned RMSE between the slope angle of the sliding plane model and the dip angle of the 3D displacements, the resulting RMSE value may be considered a measure of accuracy for the modelled sliding planes. For example, in Spitze Stei sector 5, where we expect more complex displacement processes, this RMSE value increases to 6.3°.
Direct validation of our results would be most elegant. For the case presented in this study, two indicators point towards the reliability of our results: the perfect reproduction of the observed outcrop lineament of Spitze Stei sliding plane B and agreement with the observations in borehole TB4 at the same site. A broader validation using borehole information is intended with upcoming datasets.

Relative registration methods
The rather simple approach to compensate the main displacement signal and to analyse the relative deformation of the instability instead of the absolute deformation of the slope increased the information content of the analysis significantly. At the Brienz landslide, we could identify at least 14 individual rock compartments with differing kinematical properties using coregistration method B. Systematic relative movement signals between the compartments in the lower half of the investigated perimeter indicated a toppling process here. At Pizzo Cengalo, we could show that the sum of deformations did not change when applying co-registration method C instead of co-registration method A. This was proof that the slope movement was entirely caused by rotation, i.e. a toppling process.
As the internal geometry of the point clouds does not change when applying a relative registration method, the accuracy of such relative deformation signals is equivalent to the accuracy of the absolute deformation signals. As for the angles of movement, the interpretation of these relatively registered datasets is a qualitative one. The reliability of this analysis is thus mainly shown by the fact that pronounced deformation signals occurred exactly along prominent geomorphological lineaments.

Conclusions
Point clouds of rock slope surfaces contain information reaching far beyond absolute 2D displacement data. We presented different processing methods and showed how they could contribute to a detailed kinematical model of different rock slope instabilities. The compensation of the main deformation signal in point clouds revealed unseen secondary processes in the moving rock masses. These can be visualised using 2D deformation vector fields or simple C2M algorithms. Using various strategies of relative point-cloud registration allows the distinction of toppling and sliding processes. Direct 3D deformation data prevent an underestimation of the absolute displacement rates. Mapping the individual direction components of a 3D movement such as elevation angle or azimuth can provide valuable insights into the displacement processes. Furthermore, the 3D deformation vector fields provide the basis for the integration of a 3D sliding plane model. In the case of continuous sliding planes, such models allow for the calculation of potential rockfall volumes. In case they fail to infer a consistent sliding surface, the method helps indicate kinematic complexity, such as a combination between sliding and toppling, superposition of sliding surfaces at different levels or the existence of discontinuities in the sliding plane. The kinematic observations are key in assessing possible failure scenarios and estimating the hazard potential of rock slope instabilities. Set into a broader geological and geotechnical frame, they contribute to the understanding of causes, progress and criticality of deep-seated rock slope instabilities. The data and funding of this work were provided to a large extent by the Offices for Forest and Natural Hazards of the Swiss Cantons Bern, and Grisons and the communities of Albula, Kandersteg, and the Bregaglia, to whom we extend our sincere thanks for the constructive cooperation. The Arge Alp project "Influence of permafrost on landslides and rockfalls" financed part of the Lidar measurements on Pizzo Cengalo. Marcia Phillips contributed valuable inputs to this study.

Author contribution
Robert Kenner: main author, idea, measurements, analysis, and writing; Valentin Gischig: geological input for all three sites and writing; Zan Gojcic: calculation of 3D deformation vector fields and writing; Yvain Quéau: normal integration for modeling sliding planes and writing; Christian Kienholz: Geological input for Spitze Stei, field measurements, and writing; Daniel Figi: geological input for Brienz, field measurements, and writing; Reto Thöny: geological input for Brienz, field measurements, and writing; Yves Bonanomi: geological input for Pizzo Cengalo and writing.

Funding
Open Access funding provided by Lib4RI -Library for the Research Institutes within the ETH Domain: Eawag, Empa, PSI & WSL. Indirect funding came from the Offices for Forest and Natural Hazards of the Swiss Cantons Bern and Grisons and the communities of Albula, Kandersteg, and the Bregaglia.

Availability of data and material
The point-cloud data is subject to restrictions.

Code availability
The non-commercial code for normal integration is available online https:// github. com/ yqueau/ normal_ integ ration.

Declarations
Conflict of interest The authors declare no competing interests. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.