Multi-Annual Kinematics of an Active Rock Glacier Quantified from Very High-Resolution DEMs : An Application-Case in the French Alps

Rock glaciers result from the long-term creeping of ice-rich permafrost along mountain slopes. Under warming conditions, deformation is expected to increase, and potential destabilization of those landforms may lead to hazardous phenomena. Monitoring the kinematics of rock glaciers at fine spatial resolution is required to better understand at which rate, where and how they deform. We present here the results of several years of in situ surveys carried out between 2005 and 2015 on the Laurichard rock glacier, an active rock glacier located in the French Alps. Repeated terrestrial laser-scanning (TLS) together with aerial laser-scanning (ALS) and structure-from-motion-multi-view-stereophotogrammetry (SFM-MVS) were used to accurately quantify surface displacement of the Laurichard rock glacier at interannual and pluri-annual scales. Six very high-resolution digital elevation models (DEMs, pixel size <50 cm) of the rock glacier surface were generated, and their respective quality was assessed. The relative horizontal position accuracy (XY) of the individual DEMs is in general less than 2 cm with a co-registration error on stable areas ranging from 20–50 cm. The vertical accuracy is around 20 cm. The direction and amplitude of surface displacements computed between DEMs are very consistent with independent geodetic field measurements (e.g., DGPS). Using these datasets, local patterns of the Laurichard rock glacier kinematics were quantified, pointing out specific internal (rheological) and external (bed topography) controls. The evolution of the surface velocity shows few changes on the rock glacier’s snout for the first years of the observed period, followed by a major acceleration between 2012 and 2015 affecting the upper part of the tongue and the snout.


Introduction
Because of their abundance in many poorly-glacierized mountain regions of the world [1], rock glaciers are key features to understand the response of the high altitude cryosphere to climate change [2].Rock glaciers typically flow downslope at relatively low and steady rates, exhibiting surface velocity typically ranging from a few cm to 2 m per year [1].In the Alps, synchronous inter-annual fluctuations of surface velocities have been observed on many rock glaciers, suggesting a control of the ground thermal state on the deformation rate of the permafrost [3,4].Some cases of rock glaciers' destabilization were also observed in various mountain chains, such as in the Alps [5,6] or in the Andes [7], becoming a potential source of hazards [8].These temporal changes in rock glaciers' deformation patterns pointed out various mechanical processes involved in their behavior [9].In this context, developing a sound knowledge of the spatial and temporal variabilities of the surface velocity of rock glaciers, as a proxy of their deformation and an indicator of their climatic controls, is needed.It is indeed a prerequisite for future efforts aiming at analyzing their physico-mechanical behavior and the driving factors, among which the climate probably plays a substantial role [3,[10][11][12].
Because investigating the internal deformation of rock glaciers is very costly and limited to a single one-dimensional vertical profile along boreholes, the surface displacement is the most easily "accessible" variable for characterizing rock glacier dynamics and deformation patterns.The kinematics of the topographical surface integrates the overall vertical profile of deformation, including creeping and shearing, and their changes in space and in time are valuable information for understanding the factors that rule the deformation.Surface kinematics of rock glaciers are generally measured and monitored either remotely, with terrestrial, aerial or spatial vectors (see, e.g., [13,14]), or with in situ measurements, by geodetic methods (total station, differential GPS; see, e.g., [15]).The first set of approaches has the advantage of a wide spatial coverage, using imagery with a pixel resolution of a few decimeters to a few meters (ground sample distance (GSD)).It is nevertheless limited by the temporal availability of datasets, typically several years to decades.In the second set of methods, the surface topography can be acquired more frequently (as often as an operator can go out in the field to carry out surveys), but it results in poorly representative data in spatial terms, because it consists of punctual information (generally a marked block at the surface that can be tracked from year-to-year).Other methods for continuous measurements and monitoring of ground surface changes, like continuous L1-GPS [16] or fixed camera [12], are allowed by very recent technological developments, but remain up to now costly and not so easy to implement.
Terrestrial and airborne laser scanning (TLS, ALS) and terrestrial photogrammetry (or structurefrom-motion/multi-views stereophotogrammetry (SFM-MVS)) constitute intermediary approaches.Because of their relative simplicity in handling, these approaches easily allow one to generate (very) high-resolution digital elevation models (DEM).In addition, their repetitivity can be annual or even at smaller time intervals, mostly depending on the terrain characteristics and the range of surface velocities.So far, these methods have rarely been used for rock glacier studies [17][18][19][20], and no investigations have yet addressed the potential of those methods for long-term monitoring of the surface kinematics.

Aims and Objectives
In order to investigate the multi-temporal kinematics of rock glaciers and, ultimately, to better understand their climatic controls, we took advantage of an easily-accessible site with one of the longest time-series of ground-truth data (since the early 1980s; see [21]).On this active rock glacier located in the French Alps, we carried out several TLS, ALS and SFM-MVS surveys.The objectives of the present study are to describe and evaluate very high resolution DEMs of the studied rock glacier, to determine its surface velocity fields at various temporal scales, to evaluate the spatial characteristics of the rock glacier kinematics and to discuss the geomorphological significance of the results.
After presenting the study site (Section 3) and the methods used (Section 4), we describe the datasets (Section 5) and their analysis (Section 6).These last two sections include: (i) an evaluation of the data quality by estimating the vertical and horizontal errors; (ii) the determination of the surface displacement fields at various temporal scales; and (iii) the analysis of the surface displacement changes through time and space.Finally, the significance of those results in terms of possible causes of the spatio-temporal changes of rock glacier surface displacements is discussed (Section 7).

Main Characteristics
The Laurichard active rock glacier is located in a transitional climate zone between the northern French Alps, affected by an oceanic regime, and the southern French Alps, mostly subjected to the Mediterranean influence.Using regionalization results from the SAFRANmodel provided by [22], the mean zero-isotherm altitude (ZIA) in the Laurichard area over the 1958-2002 period can be estimated to be around 2450 m asl; whereas the mean annual sum of precipitation is around 1000 mm at 1800 m asl, and the mean snow cover duration is about 220 days per year [23].
The debris that constitutes the rock glacier is predominantly of granitic lithology, falling down from a densely-fractured 600 m-high rock wall.This gravitational activity feeds the rock glacier with relatively large boulders, which constitute an open-work debris accumulation that largely covers the slopes around.
The Laurichard rock glacier has a rather simple morphology, with an 800 m-long and 100-200 m-wide tongue, flowing down between 2650 m a.s.l. at the rooting zone (foot of the rock wall) and 2450 m a.s.l. at the snout.The activity of the rock glacier is evidenced by: (i) the presence of transversal ridges and furrows; (ii) a steep and unstable latero-frontal talus; and (iii) locally, by the looseness and the instability of the open-work active layer.

Previous Studies and Monitoring Activity
The Combe de Laurichard was chosen by Francou [24] for achieving different geomorphological studies of scree slopes dynamics.Among other activities, Francou originally marked more than 40 blocks, among which 15 have been regularly measured since 1984 up to present (yearly since 1999, every 1-3 years before).Initially, blocks were positioned using a tacheometer [21], then a total station was used from 1992-2013, and since 2012, a differential dual-frequency GPS has been used in real-time kinematic (RTK) mode.This still ongoing long-term survey has allowed quantifying the annual surface velocity of the Laurichard rock glacier along its main flow line for more than three decades [25,26] with an overall centimetric accuracy.
Comparison between the 2005 TLS datasets and a digital elevation model (DEM) representing the topography in 1975 (from digitizing the 10-m contours of the 1:25,000 IGN topographic map) shows that a large part of the rock glacier movement occurs on its right-hand side [20].To account for this specific pattern of the rock glacier flow, new marked blocks were therefore added in 2008 to the original network of points (Figure 1).

Datasets Collection
Since 2005, we conducted several measurement campaigns on the rock glacier aiming at generating DEMs that fulfill the following requirements: (i) a wide spatial coverage, at least the lower half of the rock glacier; (ii) a spatial resolution finer than the coarser blocks at the surface, the size of which is roughly estimated to be greater than 50 cm; and (iii) a quality that will allow at least an annual inter-comparison of the DEMs.
For that purpose, we performed several TLS and SFM-MVS campaigns.We also benefited from an ALS survey performed in 2012.After post-processing the initial datasets (see the following section), six DEMs have been generated at a resolution ranging from 5 cm-50 cm.
All DEMs were referenced in the RGF93 grid coordinate system based on the Lambert 93 conformal conic projection and the RAF09 geoid model for elevation reference (standard coordinate reference system of the French IGN, Institut National de l'Information Géographique et Forestière).

TLS Datasets
For the TLS surveys, we used an Optech R Ilris 3D scanner (Teledyne Optech,Vaughan, ON, Canada) (in 2005, 2006 and 2011) and an Optech R Ilris 3D-LR (long range) scanner (Teledyne Optech, Vaughan, ON, Canada) (in 2013).Details on the measurement capacities of the devices are described in Table 1.
For each of the four TLS surveys, between one and three stations were used to minimize the areas not viewed by the laser.
Between 5 and 12 artificial (35 cm-diameter spheres fixed on stable and moving rocks) and/or natural targets were measured during each campaign, using dual-frequency GNSS receivers, in differential mode (DGPS), allowing a relative accuracy of a few centimeters.The flight was performed by Sintégra R on 16 and 17 August 2012, with a Riegl LMS Q680i laser scanner on board a Cessna 206, flying at a nominal elevation above the ground of 610 m.The sky was clear with limited wind (<10 kts).The point density after filtering was 7 pts/m 2 .The embedded dual-frequency GNSS and inertial motion unit, coupled with a GNSS base station allow a positional accuracy of 5 cm and an angular accuracy of 0.01 • of the laser scanner position at each time during the flight.
Sintégra performed the post-processing (point cloud modeling, filtering, DEM generation and orientation), as well as the evaluation of the dataset by comparing to in situ GNSS measurements (45 points on a flat artificial surface, yielding a vertical accuracy (RMSE) of the 3D point cloud within 2 cm) [27].

SFM-MVS Datasets
On 16 September 2015, 27 pictures were taken from the Pyramide de Laurichard, i.e., facing the rock glacier (cf. Figure 1), with a Nikon D800 reflex camera, equipped with a 35-mm fixed-focal lens.The weather was cloudy, with mostly diffuse solar radiation, allowing a rather constant lighting of the terrain and homogeneous white balance distribution during the survey.
The initial raw (.nef Nikon file format) images were converted to .tiffformat using XN Converter freeware.The images were processed in the commercial Photoscan Agisoft R software, with the following procedure: (1) selection of images using a quality index (7 images removed); (2) camera alignment and markers retrieval; (3) dense cloud processing; (4) export DEM at a 10-cm pixel resolution.Georeferencing of the 3D model is based on DGPS measurements of the artificial targets (spheres) and additional natural targets (blocks) visually detected on the 2012 ALS-derived DEM.

Co-Registration of the DEMs
Once each DEM had been individually georeferenced, they were all co-registered on stable areas with the 2005 DEM used as a reference.This was performed using a standard best-fit alignment procedure in InnovMetric R Polyworks software (InnovMetric Software Inc., Quebec, QC, Canada), which minimizes by a least-square method the distance between the two 3D models.The accuracy of the co-registration is assessed by the statistics of the error distribution (Table 2).In order to compute 3D displacements of the rock glacier between each pair of DEMs (2006-2005, 2011-2006, 2012-2011 and 2015-2012; note that the 2013 DEM appeared as too noisy and unfavorable for image correlation), the elevation grids were processed in the IMCORRmodule within SAGA-GIS [28,29].The algorithm retrieves pixels pattern between two georeferenced images and produces shapefiles (points and lines) containing the 2D and 3D displacements.Three parameters can be set: the search chip size, the pattern size and the grid sampling interval.
All the original DEMs were resampled at a 10-cm resolution and set in a common 2201 × 2951 pixels grid system.Several parameter combinations were tested, the best scores being achieved for the following parameter values: search chip: 64 or 128 pixels; pattern: 16 or 32 pixels; grid interval: 3 or 4.
The surface displacement datasets were manually cleaned due to border effects between grids with different local spatial coverage and artifacts (e.g., related to random local similarity of the coarse blocky surface) and then exported into a point-shapefile.This latter dataset was then interpolated (inverse distance weighting method, with a search radius of 10 m and a Gaussian weighting) into a 2-m resolution grid.

Comparison of Computed Surface Displacements with Ground-Based Measurements
As mentioned previously, 35 blocks located along three lines crossing the rock glacier (Figure 1) have been annually measured by the Parc national des Ecrins since 1999 [25].Interannual displacements are therefore available for those points and can be compared, in the three dimensions, to the surface displacements (distance and direction) computed from images' correlation.

Evaluation of the Vertical Accuracy
To assess the vertical accuracy of the datasets, we performed a pixel-to-pixel difference of the successive pairs of DEMs on stable areas (Table 3).For the 600-800 thousands points compared, the absolute maximal elevation difference is below 2.5 m, with a standard-deviation (SD) and a root-mean-squared error (RMS) within 15 and 27 cm.As expected, the distribution of the DEMs' difference is centered (the average equals zero), except for the 2015-2013 DEMs comparison, for which the 2015 model is locally suffering from broad-scale disturbances probably due to an uncorrected radial distortion inherent to an insufficient calibration in the SFM-MVS process [30].

Evaluation of the Displacement Quantification
A preliminary visual inspection of the shaded-relief DTMs made from the 10-cm (2005,2006,2011,2013,2015) and 50-cm (2012) DEMs (cf.Supplementary Material) clearly shows the coherent patterns of the rock glacier deformation: stable areas outside the rock glacier, main central flow line along the maximal slope with reduced displacements towards the sides, higher velocity in the steepest upper part of the rock glacier and a bend of the movement towards its right orographic side.
Mapping the surface displacements derived from image correlation and from ground measurements shows a good agreement between the two approaches (Figure 2).The direction of the vectors is satisfactorily reproduced, as well as the their amplitude with an RMS error of 0.16 m/yr (Figure 3).

Spatio-Temporal Variation of the Rock Glacier Velocity
Though varying in amplitude, the spatial distribution of the movement remained rather homogeneous in space, with a very fast zone (>1 dm/yr) in the upper surveyed area and a very slow area (a few m/yr or less) on the left flank of the tongue (Figure 4).The orographic-right bend is also visible at every period, as well as marked border effects with a clear velocity lowering from the center flow line to the lateral sides.A localized lack of displacement vectors on the 2006-2011 datasets is due to a shallow dry debris-flow that took place over the steepest part of the rock glacier.This drastic morphological change, visible on the shaded-relief DTMs, impeded pattern retrieval by the image correlation algorithm.The mean ground-measured velocity (bold dashed line in Figure 5) and the velocity fields computed from image correlation over the whole rock glacier surface (red bold line in Figure 5) both display a close mean value, suggesting that the design of the annual topographic survey is well representative of the average rock glacier movement.Deviations of ±5-15% are also observable between both curves.Though this is within the error range of the velocity derived from image correlation, it can nevertheless be assumed that this could come from the heterogeneous velocity patterns: temporary acceleration or deceleration of specific areas of the rock glacier would be captured by the glacier-wide measurements derived from DEM correlation, but not by the network of blocks surveyed.The spread (standard deviation) of the image correlation-derived values also shows the overall large spatial variability of the velocity, especially changing between the rapid central flow line and the slower sides.
The longitudinal velocity profile (Figure 6

Deformation Patterns
The deformation pattern is here quantified by the relative spatial variation of velocity along the longitudinal profile (Figure 7), which mostly allows one to identify the compressive and extensive sectors along the profile.A compression area is visible near the front (horizontal distance > 200 m), due to the decrease in velocity just uphill of the front and the consequent buttressing effect.Higher up along the profile, a lateral shift of the compressive-extensive zone seems also observable between the two periods considered in our measurements, amounting to approximately 10 m, in coherence with the longitudinal translation of the ridges and furrows.Some significant changes, between the two periods, of the longitudinal velocity variation can be seen around 115 and 125 m (−20%) and around 130 and 140 m (+15%).

Representativeness of the Estimated Surface Velocity
When statistically comparing the variability of the velocity measured from topographic surveys of marked blocks and the one of the velocity derived from image correlation of very high-resolution DEMs, we notice that: (i) the mean values of spatial variability from both methods are very close (0.21 m/yr for velocity derived from DEMs, 0.23 m/yr for velocity derived from topographic surveys), which means that the network of blocks is well representative of the whole behavior of the rock glacier; (ii) the temporal variability is also quite similar; therefore, even if the time interval between DEMs can be longer than the annual time span between two topographic surveys, a similar temporal signal is captured by both methods; (iii) the non-stationarity of the two methods is well identified: the mean velocity derived from DEMs increases, as well as its variability.

Topographic Control and Rheological Properties
No clear control of slope (derived from the 2012 DEM resampled to 2 m) on surface velocity can be detected when considering the entire profile (Figure 8).In order to account for the possible effect of surface roughness (many blocks are bigger than the 2-m DEM pixel size) and of the interannual variation of velocity, we also compared mean annual surface velocity measured with GPS between 2003 and 2016 with slopes derived from a 10-m DEM: similarly, no clear relationship appeared.Nevertheless, when considering only the 50-150-m section of the longitudinal profile (plain line section in Figure 7), a significant relationship is found for slope higher than 20 • (Figure 8).This specific pattern may be related to the specific rheological characteristics of the ice-debris mixture in a section of the profile that is mostly affected by creep, without acceleration and deceleration related to slope variation along the main flow line.

Limits and Perspectives of the Methods
Very high resolution DEMs, obtained from laser scanning or photogrammetry, have already been used for rock glacier studies.For instance, [31] used multi-temporal LiDAR DEMs to derive velocity fields on a very fast rock glacier (surface velocities 2-7-times higher than the Laurichard case), whereas [32] compared TLS-derived DEMs to airborne orthophoto on a variety of moving landforms, including rock glaciers.Nevertheless, our study is the first to assess the accuracy of LiDAR-and SFM-MVS-derived DEMs for quantifying multi-temporal surface displacements of rock glacier with 'normal' velocities (<2 m/yr).
We showed that interannual surface displacement fields of the Laurichard rock glacier can be derived in a rather straightforward way using high resolution DEMs.Though it is less precise than the geodetic survey of marked blocks, it is a complementary approach that allows describing precisely the spatial characteristics of the displacements.The combination of various sources of DEMs (terrestrial LiDAR, aerial LiDAR, SFM-MVS) is also possible, though it may complicate the selection of correct parameter values for the best retrieval of surface movement with the image correlation algorithm.One limit of this approach is the incomplete spatial coverage of terrestrial-based surveys, which leads to a lack of data in masked/hidden areas due to the limited view from the terrain.A good alternative could be, as proposed by [18], to use an unmanned aerial vehicle to produce a high resolution DEM that covers most of the study area.

Conclusions
By evaluating the very high resolution DEMs and surface displacement maps derived from those datasets, we can conclude that the various methods used to acquire the initial data (TSL, ALS, SFM-MVS) and the image correlation employed for further deriving displacement fields from DEMs are appropriate for investigating the multi-temporal kinematics of an active rock glacier.The accuracy we achieved for the DEMs, with horizontal and vertical errors around 0.15-0.3m (excluding the ALS survey) and for the velocity fields (error of 0.16 m/yr), would also even allow one to quantify changes over a time interval shorter than a year, for instance to quantify displacement during summer.If both methods show similar spatial and temporal variability and are thus adapted for monitoring the pluri-annual velocity changes, temporary and localized changes in velocity, not captured by the topographic surveys of the blocks, may also occur as suggested by the differences between mean values derived from both methods.Therefore, in order to investigate the seasonal behavior of rock glaciers and to detect localized changes, we suggest to perform, every 1-3 years, TLS, ALS or terrestrial/drone SFM-MVS surveys in order to produce 10-50-cm DEMs, in addition to the annual topographic survey.

Figure 1 .
Figure 1.Map of the study area with the location of: (1) the targets (spheres with a diameter of 30 cm, anchored on stable terrain) used for georeferencing; (2) the annually surveyed marked blocks (DGPS); and (3) the cameras for the 2015 survey (symbols are oriented with the optical axis of the lens).The lower-left inset shows the general location of the study area in the French Alps (source: Open Street Map).

Figure 2 .
Figure 2. Surface velocity (m/yr) on the Laurichard rock glacier for the period 2011-2012 quantified by image correlation on DEMs (red lines) and by ground topographic (total station and DGPS) measurements of marked blocks (black arrows).

Figure 3 .
Figure 3. Surface velocity derived from image correlation (computed from the values of the 10 closest points around the actual position of the marked blocks) versus 3D surface velocity measured on marked blocks.

Figure 4 .
Figure 4. Maps of the surface velocity of the Laurichard rock glacier at three different time periods, as derived from image correlation using HR DEMs.

Figure 5 .
Figure 5. Evolution of the Laurichard rock glacier mean velocity since the early 2000s, derived from ground topographic measurements (total station and DGPS; bold black dashed line) and from image (DEM) correlation (red lines).

Figure 6 .
Figure 6.Variation of velocity along the longitudinal axis of the Laurichard rock glacier between 2005 and 2015, derived from image correlation.The upper panel shows the velocity profile for each period, and the lower panel depicts the variations of velocity between the different periods.

Figure 7 .
Figure 7. Relative spatial variation of velocity along the longitudinal axis of the Laurichard rock glacier for the slowest period (2005-2006, in black) and the fastest period (2012-2015, in red).Positive values indicate an increasing velocity towards downslope (extensive deformation), whereas negative values indicate a slowing down (compressive deformation).The plain line is used in the central section of the profile, which is referred to later on (see Figure8).

Figure 8 .
Figure 8. Velocity against slope along the longitudinal profile (circle), with values for the central section (horizontal distance between 50 and 150 m; black dot).The red line is the regression fitted to the values for the central section for a slope angle above 20 • .

Table 1 .
Technical overview of the terrestrial laser scanners used in the study (from the Optech R factsheets).

Table 2 .
Accuracy of the datasets, with absolute values based on the manufacturer factsheet (TLS) or software processing results (SFM-MVS with Agisoft Photoscan), and relative values based on the co-registration error of each DEM on 2005 DEM.* Data according to the manufacturer factsheets.

Table 3 .
Vertical differences (in meters) between successive pairs of DEMs in stable areas around the rock glacier (sub-vertical rock face, flat Quaternary moraine and gently inclined fossil rock glacier).